ION SOLVATION D Y N A M I C S IN B I N A R Y S O L V E N T S By Tyler James Frederick Day B. Sc. (Chemistry) TJ .B.C. -O.U.C. , 1994 A T H E S I S S U B M I T T E D IN T H E R E Q U I R E M E N T S D O C T O R O F P A R T I A L F U L F I L L M E N T O F F O R T H E D E G R E E O F P H I L O S O P H Y in T H E F A C U L T Y O F G R A D U A T E S T U D I E S (Department of Chemistry) We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F B R I T I S H C O L U M B I A June 1999 © Tyler James Frederick Day, 1999 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department The University of British Columbia Vancouver, Canada Date XL ir/rn DE-6 (2/88) Abstract The dynamics of selective ion solvation in binary solvents is investigated using molec-ular dynamics simulations. Initial studies of Stockmayer solvents focus on the dependence of the usual solvation response function, S(t), on solvent composition and on the relative polarity of the solvent species. Particle solvation response functions, P(t), are also introduced which describe the compositional relaxation of the first solvation shell. It is shown that the selective solvation process can be well described by a simple phenomenological model based on the ideas of elementary chemical kinetics. This model is useful and helps in the identification of two distinct timescales present in the selective solvation process. These are associated with a rapid electrostriction step during which the total number of particles in the first shell increases to its equilibrium value, and a slower spatial redistribution process during which the composition of the first shell achieves equilibrium. The redistribution phase depends on the rate of mutual ion-solvent diffusion and also on the rate of particle exchange between the first and second shells. A detailed analysis of the exchange process indicates that exchanges occur on virtually a one-to-one basis with the insertion of a stronger dipole into the first shell being mirrored by an almost immediate ejection of a weaker one. An extension to more realistic systems consisting of both water-methanol and water-dimethylsulfoxide (DMSO) mixtures was also carried out. It is found that, despite the presence of hydrogen-bonding, the dynamical behaviour in water-methanol systems is very similar to that observed for simple Stockmayer solvents. For the water-DMSO mixtures, however, it is found that the relative sizes and geometries of the solvent species ii can have a substantial influence on the preferential solvation process. For positively charged solutes in water-DMSO the physical picture does not differ greatly from the water-methanol case. However, for negative solutes the DMSO component shows only a weak response to the charge, and the solvation process consists largely of water molecules moving slowly to the solute through an essentially static DMSO medium. The results also illustrate that the usual solvation response function, which depends on the total solute-solvent energy, is not a very sensitive probe of selective solvation dynamics. This insensitivity has since been noticed in recent experimental studies, and an alternative function that appears to be more sensitive to the solvent composition near the solute has been suggested. A comparison of simulation results with time dependent density functional theory is also carried out. It is shown that the nonlinear version of the theory is necessary in order to obtain a good quantitative description of selective solvation dynamics. The linear theory is only of qualitative value. Also, attention is drawn to a previously little appreciated problem which arises when one attempts to compare time dependent den-sity functional theory with computer simulation or experimental results. The difficulty involves matching the theoretical and absolute timescales. Finally, an extension of the theory to Stockmayer systems reveals that the theory in its present form is insufficient to capture the behaviour of these systems. i i i Table of Contents Abstract ii List of Figures vii List of Tables xii List of Symbols xiii Acknowledgement xv 1 Introduction 1 2 Simulation Methodology 7 2.1 Statistical Mechanics 9 2.2 Intermolecular Potentials 14 2.2.1 The Multipole Expansion . . 17 2.2.2 Interaction Site Models 19 2.3 Periodic Boundary Conditions 21 2.4 Molecular Dynamics 25 2.4.1 Molecular Orientation 26 2.4.2 Time Integration 28 2.4.3 The Equations of Motion • • • • 3 1 3 Solvation Dynamics in Simple Solvents 34 3.1 Model and Computational Details 34 iv 3.2 Results and Discussion 37 3.3 A Simple Phenomenological Model 48 3.4 Further Results and Remarks on Mechanism 56 3.5 Summary 63 4 Solvation Dynamics in Realistic Solvents 66 4.1 Model and Computational Details 67 4.2 Results and Discussion 70 4.2.1 Water-Methanol 72 4.2.2 Water-Dimethylsulfoxide 81 4.3 Summary 93 5 A Comparison with Time Dependent Density Functional Theory 95 5.1 Theoretical Overview 96 5.1.1 The Self-Diffusion Constant 99 5.2 Results and Discussion 104 5.3 Summary 117 6 Summary and Conclusions 121 Bibliography 125 Appendices 129 A The Potential Due to a Point Dipole 129 B The Ewald Sum 132 C The Reaction Field Method 138 v D The Generalized Langevin Equation 140 vi Lis t of Figures 1.1 The structure of Coumarin 153 3 2.1 Form of the standard Lennard-Jones potential for the interaction energy of two spherical particles 16 2.2 The effect of the Euler angle rotations on a water-like molecule 27 3.1 Solvation response functions for systems of different size. The solid curve is System 4 (256 particles) and the dotted curve is System 7 (500 particles). 38 3.2 Solvation response functions for different boundary conditions. The solid curve is System 5 (Ewald method) and the dotted curve is System 6 (re-action field method) 39 3.3 Solvation response functions for System 1 (dotted) and System 2 (dashed) along with that of System 0 (solid) for comparison 40 3.4 The number of solvent particles of each type (labelled W and S) together with the total number present in the first solvation shell for System 1. The corresponding solvation response function is also shown 42 3.5 Solvation response functions for Systems 3, 4, and 5 along with that of System 0 for reference 45 3.6 The number of solvent particles of each type (labelled W and S) together with the total number present in the first solvation shell for System 5. The corresponding solvation response function is also shown 46 3.7 Pair distribution functions for System 3 with unionized (a) and ionized (b) solute. The curves are for Solvent W (solid) and Solvent S (dotted). . 47 vii 3.8 Particle solvation response functions obtained from simulation (dotted) and from fitting of Eqs. (3.3.8) (solid) for System 2 53 3.9 Particle solvation response functions obtained from simulation (dotted) and from fitting of Eqs. (3.3.8) (solid) for System 3 54 3.10 Particle solvation response functions obtained from simulation (dotted) and from fitting of Eqs. (3.3.8) (solid) for System 5 55 3.11 A comparison of ramping effects on the composition of the first solvation shell as a function of time. The dotted and solid curves are for Systems 2 and 5, respectively. The dashed curve was obtained using Eq. (3.4.1) and the long-dashed curve using Eq. (3.4.2) 59 3.12 Modified shell counts as described in the text for System 3 (a) and System 5 (b). The curves are (Nw(t)) (W), (Ns(t)) (S) and (N(t)) (Total). . . . 61 3.13 A comparison of the influence of varying solvent dipole moments. The curves are for System 5 (solid), System 2 (dotted) and System 8 (dashed). 62 4.1 A schematic of the 4 models employed in this chapter. The relative site sizes are drawn to scale, while molecular geometry is approximated. The hydrogen sites possess no Lennard-Jones component, and are represented by dashed circles 71 4.2 Solute-solvent radial distribution functions for the systems consisting of 50% water and 50% methanol 74 4.3 The solvation response functions, Sj(t), particle solvation response func-tions, Pj(t), and first shell coordination numbers, (Nj(t)), for System WM50—. The labels correspond to water (W), methanol (M), and the total (T) : . . . . 75 viii 4.4 The solvation response functions, Sj(t), and particle solvation response functions, Pj(t), for System WM50+. The labels are as in Fig. 4.3. . . . 78 4.5 Solute-solvent radial distribution functions for the systems consisting of 10% water and 90% methanol 79 4.6 The solvation response functions, Sj(t), and particle solvation response functions, Pj(t), for System W M 1 0 - . The labels are as in Fig. 4.3. . . . 80 4.7 Solute-solvent radial distribution functions for the systems consisting of 50% water and 50% DMSO 83 4.8 Solute-solvent site-site pair distribution functions for DMSO in Systems WD50N, W D 5 0 - and WD50+ 85 4.9 A snapshot showing all (a) DMSO and (b) water molecules present within a distance of 7.0 A of the negatively charged ion. The colours correspond to methyl (black), sulfur (orange), oxygen (blue), hydrogen (red), as well as the ion (green) 86 4.10 A snapshot showing all (a) DMSO and (b) water molecules present within a distance of 7.0 A of the positively charged ion. The colours are as in Fig. 4.9 87 4.11 Solute-solvent energy components, (Uj(t)), and first shell coordination numbers, (Nj(t)), for System WD50—. The labels correspond to water (W), DMSO (D) and the total (T) 88 4.12 The solvation response functions, Sj(t), and particle solvation response functions, Pj(t), for System WD50—. The labels are as in Fig. 4.11. . . . 91 4.13 The solvation response functions, Sj(t), and particle solvation'response functions, Pj(t), for System WD50+. The labels are as in Fig. 4.11. . . . 92 ix 5.1 The equilibrium pair distribution functions for (a) the initial state and (b) the final state, as obtained from simulation with xs = 0.25 and €Xs/€xs = 6. The solid curves are for Solvent W and the dashed curves for Solvent S. 107 5.2 A comparison of particle solvation response functions obtained from simu-lation and linearized theory for System 1. In the upper, middle and lower panels the D/Db ratio is 1, the value given by the theory of Araki and Munakata, and the value obtained by fitting, respectively. In each panel, from top-to-bottom, the curves are Pw(t), Ps{t) and Pr(t) 108 5.3 A comparison of particle solvation response functions obtained from simu-lation and linearized theory for System 2. Curves and panels are as in Fig. . 5.2, although the Pr(t) curve has been omitted, as the statistical noise prevented any useful consideration 109 5.4 A comparison of particle solvation response functions obtained from sim-ulation and linearized theory for System 3. Curves and panels are as in Fig. 5.2 110 5.5 A comparison of particle solvation response functions obtained from simu-lation and linearized theory for System 4. Curves and panels are as in Fig. 5.2, although the Pr(t) curve has been omitted, as the statistical noise prevented any useful consideration I l l 5.6 A comparison of particle solvation response functions obtained from sim-ulation and nonlinear theory for System 1. Curves and panels are as in Fig. 5.2 112 5.7 A comparison of particle solvation response functions obtained from simu-lation and nonlinear theory for System 2. Curves and panels are as in Fig. 5.2, although the Pr(t) curve has been omitted, as the statistical noise prevented any useful consideration 113 x 5.8 A comparison of particle solvation response functions obtained from sim-ulation and nonlinear theory for System 3. Curves and panels are as in Fig. 5.2 114 5.9 A comparison of particle solvation response functions obtained from simu-lation and nonlinear theory for System 4. Curves and panels are as in Fig. 5.2, although the Pr{t) curve has been omitted, as the statistical noise prevented any useful consideration 115 5.10 A comparison of particle and solvation response functions obtained from theory and simulation for System 2 of Chapter 3. The solid and dashed curves are simulation and theoretical results, respectively. From top to bottom the curves are Pw(t) [Ss(t)], Ps(t) [S(t)] and P(t) [Sw(t)}. The timescale is plotted in units of the rotation diffusion constant for the weakly polar species, 118 A . l Diagram for the derivation of the potential at a point P due to a point dipole moment 131 xi List of Tables 3.1 A summary of the systems simulated. The reduced dipole moments and mole fractions are given. Note that p* — 0.8 and T* = 1.35 for all systems studied 41 3.2 The composition of the first solvation shell at equilibrium. N!j and Nj (j = W, S) are the average numbers of Solvent W and Solvent S particles present within a radius of 1.3er for both ionized (i) and unionized (u) solute particles. The error limits represent one standard deviation 44 3.3 Relaxation times for selected systems obtained from fitting Eqs. (3.3.8) to simulation results 52 4.1 Model parameters for the species considered. Distances are in Angstroms, energies in kJ/mole, charges in e, and angles in degrees 69 4.2 A summary of information for the systems considered. The system names are as described in the text. r c is the cut-off distance and JVw at and iVmeth/DMS0 are the average number of particles of a particular type found within rc of the solute at equilibrium. The relaxation times, r e i and r r e were determined by fitting the phenomenological equations 72 5.1 Mole fractions, xs,w, and interaction parameters for the systems considered. 106 xii Lis t of Symbols Symbol Description c(r) Direct correlation function. D The self diffusion constant. Db The bare diffusion constant. E Electric field. g(r) The radial distribution function. G(r,t) The van Hove function. H The Hamiltonian. J(r,t) Particle flux. kB The Boltzmann constant. I Moment of inertia. k Reciprocal lattice vector. K The kinetic energy of the system. m Mass. n The number of trajectories. n Lattice vector. N The number of particles. P Generalized particle momenta. Pa{t) The particle solvation response function for species a. Q Point charge. q Generalized particle coordinates. xiii Qo,Qi,Q2,Q3 Quaternion parameters. r Intermolecular separation vector. r Intermolecular distance. Sa(t) The solvation response function for species a. t Time. T Temperature. U The potential energy of the system. V Volume. X Fractional composition. a Ewald convergence parameter. e Lennard-Jones energy parameter. e' Dielectric constant. /j, Dipole moment vector. V The gradient operator. v Fluorescence emission frequency. cj) Electrostatic potential. p Density. a Lennard-Jones length parameter. Tei,rei Electrostriction and redistribution timescales. xiv Acknowledgement I wish to thank my supervisor, Dr. Gren Patey, for his excellent guidance and support during the past five years. I would also like to thank all the postdocs, grad students, undergrads, and visitors who have come and gone over the years. They've helped make this a truly enjoyable place to work, both academically and socially. Most of all, I would like to thank my family for their constant support and encouragement. And of course, I would like to thank the Natural Sciences and Engineering Council of Canada for their financial support. xv Chapter 1 Introduction The study of condensed phase chemistry has seen a dramatic increase in recent decades. Developments in the theoretical treatment of the liquid phase, coupled with the refinement of computer simulation techniques, has led to a dramatic increase in the study of a wide variety of topics, including phase behaviour, interfacial phenomena and liquid crystals. Indeed, the current generation of computing facilities even allows for the study of biologically important macromolecules, a subject which has been beyond the scope of any simulation treatment until now. With the seemingly endless increase in computing power each year, new areas of study are constantly being opened up. One area of condensed phase study that has received significant attention in recent years is that of solvation dynamics [1-22]. Solvation dynamics may be most generally defined as the response of a solvent to a change in the characteristics of a solute parti-cle. The nature of this change may vary considerably, the more typical examples being simple ionization or electronic excitation. The excited states of many molecules possess substantially different electrostatic characteristics than do the ground states. This is typ-ically seen as a change in the magnitude, and possibly direction, of the dipole moment, although more subtle changes are also possible. Other examples of changes may include an alteration of the size or geometry of a solute molecule, leading to a local restructur-ing of the solvent. Studies such as these are essential for establishing a groundwork for the understanding of more complicated liquid phase phenomena. For instance, it is well known that the rates of most chemical reactions in solution vary considerably with the 1 Chapter 1. Introduction 2 particular choice of solvent. The reason for this is typically attributed to the solvation state of the reactants, which thereby determines how accessible they are to other reactant species. In addition, the ability of the solvent to react to the formation of an intermediate or final product may play a significant role in determining the stability of that product, thereby affecting the overall rate of reaction. The vast majority of existing theoretical and computer simulation work in this partic-ular area has been concerned with the dynamics of ion solvation in pure, polar solvents. A n initially neutral solute particle is ionized at time t = 0 , and the subsequent response of the solvent, typically monitored via the solute-solvent interaction energy, is then ob-served. Simulation studies on simple models of spherical molecules with embedded dipole moments have shown that the solvent relaxation is comprised of two phases. The initial phase consists of a rapid reorientation of all the solvent particles. Solute molecules in the presence of a strong electric field, such as that due to the newly introduced charge, will obviously have a preferred orientation which differs considerably from the average orientation in the absence of a field. Thus, following ionization, the strong orientation-dependent ion-dipole interaction results in the observed reorientation of the solvent. This is followed by a somewhat slower relaxation phase by which the solvent density about the ionized solute increases. This, again, is due to the strong ion-dipole interaction which acts to "compress" the solvent in the immediate vicinity of the solute particle, and is often referred to as electrostriction. Various theoretical approaches have also been undertaken in the study of these systems, and have yielded similar results. As with all studies of a theoretical nature, it is always desirable to have analogous experimental results with which to compare. While a direct experimental measurement of the solute-solvent energy (the quantity most often used in theoretical and computer simulation studies) is obviously not obtainable, appropriate experimental techniques have Chapter 1. Introduction 3 been developed which allow for the study of this particular phenomenon. These experi-mental studies [23-29] are typically carried out using solute molecules such as Coumarin 153, the structure of which is shown in Fig. 1.1. While species such as this are obviously more complex than the small, spherically symmetric ions typically employed in theo-C F 3 Figure 1.1: The structure of Coumarin 153. retical and computer simulation treatments, they possess two important characteristics which make them particularly useful for solvation studies. First, they are capable of being photoexcited to states which possess substantially different electrostatic characteristics than are typical of their ground state. This may simply be a substantial increase in dipole moment, as in the case of Coumarin 153 (6.55 Debyes —> 14.2 Debyes), or may correspond to complete ionization. These excited states, by necessity, must possess lifetimes which are considerably longer than the typical time required for solvent relaxation for them to be of practical use. The second important characteristic of these species is that they are fluorescent. Upon absorbing light from an external source, typically a laser, they then re-emit the energy at a characteristic frequency. This characteristic frequency is directly related to the internal electronic energy states of the molecule. Upon excitation, the solute molecule subsequently loses energy through collisions with the surrounding sol-vent, thereby decaying to the lowest vibrational mode of the excited state. Spontaneous Chapter 1. Introduction 4 emission then returns the molecule to its ground electronic state, emitting radiation in the process. The wavelength of this radiation is determined by the electronic energy level gap within the molecule. The application of an external electric field, such as that due to the environment in a polar fluid, will induce a change in this gap, and thereby a change in the emission frequency, referred to as a Stokes' shift. Thus, as the microscopic envi-ronment of the solute molecule changes, such as occurs during solvent relaxation, there is a corresponding change in the Stokes' shift. This technique allows for the convenient monitoring of the relaxation of the electric field at the solute, which is directly related to the solute-solvent energy, as measured in theoretical and computer simulation studies [4]. Indeed, experimental results obtained in pure solvents have yielded a similar picture to that from theoretical treatments, despite the often more complicated nature of the solute. While the amount of study directed towards solvation in pure systems has been con-siderable, the related topic of solvation in binary solvents has received comparatively little attention. While at first this may seem to be a trivial and unnecessary compli-cation of the study of pure systems, there are in fact considerations which make these systems significantly different. While the solvent relaxation observed in pure systems consists of rapid reorientation followed by electrostriction, one might expect that a third phenomenon will become important in binary systems. This additional factor can be seen by simply considering the equilibrium state for both a neutral and ionized solute particle immersed in a binary solvent. In all likelihood, the composition of the immediate environment of the solute particle will not be representative of the system as a whole. It is expected that a slight preference for one solvent over the other will be present. Factors such as the solute-solvent interaction will not be identical for both species, and consid-erations of molecular size and geometry may lead to packing constraints which would also influence the solvation state of the solute. This phenomenon is typically referred Chapter 1. Introduction 5 to as preferential or selective solvation. In all probability, the nature of the preferen-tial solvation (that is, the particular equilibrium solvent composition in the immediate vicinity of the solute) will be different for both the neutral and ionized case. Thus, in order to achieve equilibrium following ionization, a substantial change in composition of the solvation shells of the solute must occur, in addition to the relaxation behaviour observed in pure systems. This redistribution phase can only be accomplished through the translational displacement of solvent particles, and is expected to take place on a substantially slower timescale than the other two relaxation modes. Solvation studies in such binary systems would have important chemical applications in several areas. Many chemical reactions require a completely anhydrous environment to proceed as desired. Even trace amounts of water will either result in unwanted side reactions or even prevent reaction entirely. Typically, addition of species such as chloride salts is made to prevent such interference. Chloride ions, when added to such systems, will scavenge any water molecules present, thus preventing them from taking part in these side reactions. The degree to which this scavenging takes place, as well as its rate, is of obvious importance. Considerations of preferential solvation also play an obvious role in systems where one of the two solvents is, in fact, a reactant. The presence of the other, non-reactive solvent species will result in competitive solvation. Thus, the rate at which the reactive species can penetrate to the solute will greatly influence the rate of the overall reaction. The remainder of this thesis is divided into 5 chapters. Chapter 2 presents an overview of the statistical mechanical and computer simulation techniques required for the study of the liquid phase. Chapter 3 presents simulation results for relatively simple solvent models. This serves as a useful initial point before considering more sophisticated sys-tems. Also presented in this chapter is a simple phenomenological model which can be applied not only to the simulation results presented here, but may also prove useful in Chapter 1. Introduction 6 the analysis of experimental results. In Chapter 4, the applicability of the observations and conclusions of Chapter 3 to more realistic systems is investigated. This entails the study of ion solvation in water/methanol and water/DMSO systems, including the ap-plication of the phenomenological theory. Chapter 5 consists of a comparison between the results of computer simulation and theoretical calculations performed using time-dependent density functional theory (TDDFT) , done in collaboration with Dr. Akira Yoshimori of Kyushu University, Japan. Finally, Chapter 6 presents a summary and overview of the conclusions reached in this thesis. Chapter 2 Simulation Methodology In the study of a physical phenomenon such as that of interest here, there are many tools and methods available. One such tool which has seen an increase in both power and practicality in recent years is computer simulation [30]. With the unprecedented increase in available computing power witnessed over the past few decades, the study of formerly intractable physical problems has now become commonplace. In order to tackle these problems two complementary simulation techniques have been developed, Monte Carlo (MC) and molecular dynamics (MD). Each of these approaches provides an alternate means of obtaining equilibrium and, in the case of M D , dynamical properties of a system. M D simulations involve the straightforward integration of the equations of motion for the entire system. Equilibrium quantities can then be calculated as time averages as the system evolves. In contrast, the Monte Carlo approach provides a statistical means of probing the equilibrium state of a system. Beginning with an initial configuration, small perturbations, typically in the orientation and position of constituent particles, are made. Acceptance of each such a "Monte Carlo move" is made based on the energy of the new configuration relative to the previous one from which it was obtained. If the move yields a lower energy, the new configuration is accepted automatically. If it yields a higher energy, the move is accepted with a probability determined by a Boltzmann factor. Thus, after many such moves have been performed, the entire collection of configurations which have been generated form an ensemble corresponding to a given temperature and pressure. Equilibrium properties can then be determined by averaging over the ensemble. 7 Chapter 2. Simulation Methodology 8 The usefulness of either of these two techniques depends on the particular problem being investigated. The strength of the Monte Carlo approach lies in the versatility of the permutation algorithm. The canonical (constant NVT) ensemble can be easily sampled [30], and extension to the grand canonical ensemble (constant pVT) is easily facilitated by the allowance for creation and destruction of particles in addition to the standard moves. Conversely, the advantage of the molecular dynamics approach is a direct result of the integration of the equations of motion, thus reproducing the actual time evolution of the system. Consequently, M D simulations are able to provide information on time dependence. This can either be the simple estimation of transport properties such as diffusion constants, or the more sophisticated study of nonequilibrium phenomena. How-ever, the integration of the equations of motion which is inherent in the M D approach can also be problematic. If one makes the obvious choice of Newtonian equations of motion, the resulting system will be in the microcanonical (constant NVE) ensemble. While tech-niques for constant N V T and N P T studies have been developed, the resulting equations of motion which must be employed are by definition non-Newtonian. Consequently, care must be taken to ensure that these modified equations are not having unphysical effects on the system. This is of particular importance in the case of nonequilibrium systems. If, for instance, the temperature is being fixed, one must ensure that the method for adding and removing energy from the system is not interfering with the observed physical be-haviour. Many techniques have been developed for this, and will be discussed later in this chapter. Due to the obvious importance of time dependence in the study of solvation dynamics, this thesis will rely exclusively on the molecular dynamics approach. Chapter 2. Simulation Methodology 9 2.1 Statistical Mechanics A system of N classical particles can be fully described by a set of generalized coordinates, q = (qi, q%, ...qjv) along with their conjugate momenta, p = (pi, P2, •••P/v) [31]. In the most general case, a given molecule requires 6 coordinates and 6 momenta to describe it completely, and thus a system of N such particles requires 12N parameters to be fully characterized. The 127V dimensional space which contains all such possible states is referred to as the "phase space" of the system. For particles containing symmetry, the dimensionality of this space is decreased somewhat. For instance, a system composed of axially symmetric particles only requires ION variables, while one composed of spherically symmetric ones will require only 6N. A system evolving in time will generate a trajectory in this multidimensional space determined by the initial state point and the equations of motion. The Hamiltonian, H, which governs the system is given by [31] H(q,p) = K(p) + U(q), (2.1.1) where K(p) is the kinetic energy and U(q) is the total potential energy. If the particular equations of motion being employed are "Hamiltonian", as is the case for Newtonian mechanics, then the evolution of the trajectory is determined by Hamilton's equations q = — (2.1.2) dp p = - — . (2.1.3) aq In such a case, the trajectory evolves over a constant energy surface in phase space. That is, it passes through points for which the total energy of the system is fixed. Alternatively, non-Newtonian forms for the equations of motion are also possible. For instance, if one instead desires to maintain a constant temperature, a modified form of the equations of motion may be chosen, as will be seen later. Chapter 2. Simulation Methodology 10 Given the generalized coordinates, q(i), and their conjugate momenta, p(t), as a function of time for a given equilibrium trajectory, one can calculate the average value of a given observable quantity, A, as a time average over the trajectory [30] Of course, in terms of computer simulation it is not possible to extend to infinite time, but rather a sufficiently long trajectory is generated so as to produce a reasonable estimate of the quantity A, assuming that the system is at equilibrium. While it is conceivable that one may wish to calculate such an average in a system which is not at equilibrium, this is not normally the case. Typically, the system is evolved in time until equilibrium has been achieved, and an estimate oi A is then determined over a subsequent trajectory. The progress of a system to equilibrium is in fact monitored via Eq. (2.1.4). Various quantities such as energies and structural factors are calculated over a series of short time intervals. When no discernible drift in these quantities is present, it is assumed that the system has achieved equilibrium. The possibility of "bottlenecks" in phase-space, corresponding to metastable states, is always of concern [30]. However, for most condensed phase systems a suitably long equilibration period is sufficient to overcome this problem. Then, once equilibrium has been achieved, the desired quantities may be estimated from Eq. (2.1.4). One of the advantages of computer simulation lies in the near limitless number of properties which may be calculated. With the coordinates and momenta of all particles available as a function of time, quantities which are impossible to measure experimentally are often trivial to determine from simulation data. In the present study of preferential solvation, for instance, the equilibrium structure of the solvent is of particular interest. While overall structural information is experimentally available from techniques such as neutron diffraction [32], the contribution attributable to one particular component in a multi-component system is not readily discernible. While some techniques such as (2.1.4) Chapter 2. Simulation Methodology 11 N M R have been used in an attempt to obtain such information [33,34], the reliability of such techniques is still not well established. Additionally, these techniques provide information on only equilibrium structure, and are not capable of the high degree of time resolution required for non-equilibrium studies. However, such data is easily obtained from computer simulation. The most informative means of describing liquid structure is through the pair distribution function, g(r,Q) [30,35]. This function, in its most general form, measures the density of particles at position r with orientation $7 relative to the bulk density. Thus, regions of greater than normal density will correspond to values of g(r,Q) greater than 1, while regions of low density will correspond to values less than 1. In the present work a somewhat more simplified form, the radial distribution function g(r), is also of interest. This function is simply a measure of the density of particles at a distance r from a reference particle, ignoring orientational information. For the study of preferential solvation of interest here, the reference particle will typically be the solute particle. Thus, the solute-solvent radial distribution function provides detailed information about the distribution of the solvent species relative to the solute particle. In addition to equilibrium structural information, the availability of dynamical data is of critical importance to the study of solvation dynamics. Following a change in the properties of a solute particle (eg. ionization) the nature of the subsequent solvent re-laxation is of primary interest. A time averaged quantity such as that represented by Eq. (2.1.4) is of limited usefulness. The value obtained would be an average over the entire relaxation process and would contain no detailed information about the dynamics, or about the state of the system at a particular instant in time. Instead, the instanta-neous value of the desired observable quantity, A(q(t),p(£)), must be used. However, such instantaneous quantities are susceptible to a considerable degree of statistical noise, particularly in small systems such as those of practical use in computer simulations. To overcome this difficulty, the concept of a time average is replaced with that of a trajectory Chapter 2. Simulation Methodology 12 average [30]. Let (ri(0), r 2(0), ... r„(0)) be a series of "initial" points in phase space representing a state of the system which is far from equilibrium. Each of these points will then generate a nonequilibrium trajectory as the system relaxes to the equilibrium state. The trajectory-averaged time-dependent value of a particular quantity, < A(t) >, is then determined as an average over the trajectories themselves [30] <A{t)>=-J2An(t) (2.1.5) n i where An(t) is the instantaneous value of A at time t in trajectory n. Provided that a sufficient number of trajectories are available over which to perform the averaging, a reasonable estimate of < A(t) > can be obtained. In terms of the study of solvation dynamics, the initial points Tj are selected from the equilibrium state of the initial neutral system. The desired change in the solute-solvent interaction is then introduced (eg. ionization), and the time evolution of each of these state points to the new equilibrium state is monitored. The quantity most typically employed to monitor trajectory-averaged solvent relax-ation is the solvation response function, Sa(t), defined as _ (U{t) - U(oo)) ba[t)~ (u(o)-u(oo)y [ 2 - L b ) where U(t) is the solute-solvent energy at time t, U(0) is the value at time t = 0, and [/(oo) is the final equilibrium value. The solvation response function is defined in such a way that it has an initial value of 1 at t = 0, and a final value of 0 as t —> oo. Thus, it becomes immediately obvious for a particular S(t) curve how far the system has progressed towards equilibrium. Additionally, the dimensionless nature of Eq. (2.1.6) allows one to easily compare the results for two different systems without having to account for differences in absolute energy scales which may arise due to differences in model parameters. Chapter 2. Simulation Methodology 13 As mentioned previously, the solute-solvent energy cannot be obtained experimentally. Instead, the time-dependent Stokes' shift, u(t), is measured. It is generally assumed [4] that this quantity behaves similarly to the solute-solvent energy, and thus the experi-mental analog of the solvation response function, the spectral response function S„(t), i/(0) — v(oo) where u(t) is the peak in the fluorescence spectrum at time t, is determined. Thus, Eqs. (2.1.6) and (2.1.7) allow for the direct comparison of theoretical and experimental results. Although the interaction energies have proven sufficient in the study of solvation in one-component solvents, the extension to binary solvents means that additional informa-tion may be of interest. While the energetic relaxation may capture the nature of the overall relaxation, it contains no explicit information about the particular solvation state of the solute. To obtain such information, a structural quantity is required. While the time-dependent radial distribution function may seem to be the ideal candidate, it turns out to be impractical for nonequilibrium studies. To obtain a reliable time-averaged g(r) in an equilibrium system requires that a relatively lengthy trajectory be generated. The number of trajectories which would be required to obtain a meaningful trajectory-averaged radial distribution function for a nonequilibrium system is computationally prohibitive. Instead, a more coarsely grained quantity which is less sensitive to statistical noise, but still contains the relevant information is required. The most obvious choice is the compo-sition of the first solvation shell of the solute as a function of time. More accurately, it is the number of solvent a particles present within the first solvation shell of the solute at time t, Na(t). For the purposes of the current study, the first solvation shell is defined to be that region within a distance r c of the solute, where rc is dependent on the particular system being studied and is typically taken to be the first minimum in the solute-solvent pair distribution function. This quantity allows for the monitoring of the solvation state Chapter 2. Simulation Methodology 14 of the solute as it evolves in time. A function analogous to the solvation response function can then be defined, related to the composition of the first solvation shell. This particle solvation response function is given by _ (Na(t) - Na(oo)) ~~ / AT (n\ AT7^T\' (2.1.8) {Na(0) - Na(oo)) and has the same benefits of the solvation response function, in that it is a dimensionless quantity which contains no system specific information such as absolute coordination numbers, etc. 2.2 Intermolecular Potentials For a system of N interacting particles with positions TN, the general form for the total configurational potential energy, U, is given by U(rN) = J2 + E + E ri> r*) + - (2- 2- 1) i ij ijk where the first term is dependent on individual particle positions (as might be important in the presence of external fields, etc. ), the second term represents a pairwise sum, the third a triplet sum, and so forth. While ideally it would be desirable to include as many such terms as possible, the increase in computational cost for even the three-body term is considerable. Consequently, although three-body and higher interactions may constitute a substantial contribution to the total configurational energy in dense systems, most simulation studies, including those presented in this thesis, are restricted to pairwise additive terms only. Instead, most models employed in these studies are devised in such as way as to implicitly contain important many-body contributions such as polarization in the potentials. This is accomplished by assuming a typical liquid density and adjusting the parameters in the pairwise additive potentials so as to reproduce the effect of including three-body and higher order terms. Chapter 2. Simulation Methodology 15 The form for the pair potential itself may take many forms. For simple, nonpolar spherical particles (such as argon), the potential may be well represented by both a short-range repulsive, and long-range attractive term. The attraction, resulting from dispersion forces between particles, is known experimentally to vary as r - 6 , where r is the intermolecular distance. The repulsive contribution due to the overlap of electron orbitals as two particles are brought into close proximity results in a steeply repulsive barrier at short separations. There is no unique functional form for this repulsion, but it is most commonly represented by either an exponential or r~ 1 2 term. For the purposes of this thesis, the latter choice will be employed, leading to the standard Lennard-Jones (LJ) potential for the interaction of two spherical, nonpolar particles [30] where e defines the "well depth" of the potential, and o determines the effective size of the particle. A schematic of a typical L J potential is shown in Fig. 2.1. The cross-interaction between species possessing different e and o values is typically determined using the Lorentz-Berthelot rules, o\i — (crn + a22)/2 and t\i = v/ene22-While the Lennard-Jones potential is very useful and has been employed in many theoretical and simulation studies, it is obviously rather limited. Should one desire an orientation dependent potential, or one which includes long-ranged electrostatic forces, one must devise a more robust model. While several approaches have been developed, only two are of interest here. The first of these, the multipole expansion, attempts to reproduce the potential due to an arbitrary charge distribution by defining electrostatic "multipoles" located at the center of the distribution. The second, somewhat more versatile approach is the "interaction site model." This approach reproduces the overall potential by distributing various sites within the molecule, each of which is defined by its own set of potential parameters. (2.2.2) Chapter 2. Simulation Methodology 16 Distance Figure 2.1: Form of the standard Lennard-Jones potential for the interaction energy of two spherical particles. Chapter 2. Simulation Methodology 17 2.2.1 The Multipole Expansion Consider a system containing N point charges, qt, with positions r;. The potential at a point P, with position rp , due to this charge distribution is given by [36] ^ ) = E ^ S . (2-2.3) Expansion in a Taylor series yields m = E r L - • V ^ - + \ E qlrlrl : V 2 ^ - - (2.2.4) i rP . r P 2 . r P where rp = | rp | and V is the gradient operator defined as V = ° - + ^ + ^- (2.2.5) o x a y a2; with the differentiations being made at P. This expression is valid so long as the point P lies "outside" the distribution of charges. That is, r; must be significantly greater than the distance between any two charges. The first term in Eq. (2.2.4) is the contribution due to the net charge on the system. The second and third terms are contributions from the dipole moment, p, and quadrupole moment, Q, respectively, where these quantities are given by /i = Eftri (2.2.6) and Q = ^EWi- (2-2-7) These are typically referred to as "point" multipoles, due to the assumption that the field is being sampled at a distant much greater than the separations between the charges, and thus the charge distribution itself may be represented as a point. While the above derivation is the most general form for the field due to an arbitrary charge distribution, in the current work we are only concerned with the dipole term. The potential at point Chapter 2. Simulation Methodology 18 P due to a point dipole can be obtained by combining Eq. (2.2.4) with the definition Eq. (2.2.6) to yield 4>{P) = -p-V- (2.2.8) r or, by applying the gradient operator explicitly, 4>{P) = (2.2.9) Additionally, since the electric field, E, is related to the potential as E(P) = - V 0 ( P ) , (2.2.10) one can write the field due to the point dipole as E(P) = - £ . (2.2.11) This allows one to calculate the interaction energy of two point dipoles,' UOD, given by uDD = p • E = 3 ( ^ • ^ • r ) - (2.2.12) Similarly, the interaction energy of a point charge q and a point dipole, uID, is given by UiD = q<t> = q ^ . (2.2.13) It should be noted that the above derivation assumes that the charge distribution con-sists of a number of point charges. For the case of an actual molecule, these moments are calculated by replacing the sums over charges with integrals over continuous charge distributions. While the above derivation is useful when higher order multipoles are involved, it is not very instructive as to the physical meaning of the moments themselves. A n alternate description of a dipole moment proves somewhat more intuitive. Consider two point Chapter 2. Simulation Methodology 19 charges, +q and —q, separated by a vector, 1. The dipole moment of this configuration is given simply by p = ql. (2.2.14) The corresponding idealized dipole can be obtained by replacing 1 with si, q with q/s, and taking the limit as s —> 0. This has the effect of preserving the overall dipole moment, p, while reducing the spatial volume occupied by the dipole itself to a single point. The potential for a point dipole derived in this manner, as is outlined in Appendix A, is identical to that resulting from the multipole expansion. The combination of a L J particle with an embedded point dipole is typically referred to as a "Stockmayer" particle, and will be of considerable importance to this thesis. The multipole expansion method is useful for investigating the effect of charge, dipole moment, etc. on the behaviour of the system being studied, and is computationally cheap to implement. However, as mentioned above, the expansion is only really valid so long as the point at which the field is being sampled is sufficiently far from the charge distribution itself. Thus, if one intends to model a realistic molecule, such as water, in the condensed phase, the usefulness of the multipole expansion is somewhat limited. Consequently, an alternative approach must be taken if one wishes to move from simple solvent models to more realistic ones. 2.2.2 Interaction Site Models The interaction site approach [30] produces a more robust intermolecular potential at the cost of increased computational complexity. Instead of representing a molecule as a single site characterized by a set of potentials, it distributes multiple sites within a single molecule. Each site is characterized by its own potential, as well as its location, s, relative to some origin (typically the molecular center of mass). While these potentials Chapter 2. Simulation Methodology 20 may be of any form, including a multipole expansion, they more typically consist of a single point charge, a LJ-type potential, or both. While variants of the interaction site model have been developed which allow for movement of these sites relative to the origin of the molecule so as to allow for bond vibrations, rotations, bending, etc. , this thesis will be concerned with rigid models only. While this may appear to limit the choice of molecules to exceedingly simple species, certain simplifications are possible. For instance, in the modelling of methanol, it would appear that internal rotations about the C-0 bond must be accounted for. However, if one uses a single interaction site to represent the entire methyl group, as opposed to one for carbon and one for each hydrogen, the problem of internal rotations is avoided. This "unified atom approach" assumes that the internal rotations about the C-0 bond have a negligible effect on the intermolecular potential of the methanol molecule itself. As the typical timescales of these bond rotations is substantially faster than that of molecular rotation and translation, the use of a "rotationally averaged" single site potential for the methyl group is sufficient. Thus, a molecule such as methanol could be represented by a simple three site rigid model. The total configurational potential energy of a system constructed from interaction site models is simply the sum over all possible pairs of sites treated as if they were independent particles. However, one exception must be made. The interaction of two sites within the same molecule is omitted. The inclusion of this intramolecular potential would be an unphysical, albeit constant, contribution to the total energy. Forces acting on sites are also calculated pairwise, as though the sites were independent particles. The total force, F, acting on a molecule is then the vector sum of all the forces on its constituent sites, sites F = £ fi, (2.2.15) Chapter 2. Simulation Methodology 21 where f; is the force on site i. Similarly, the torque on a molecule can be calculated as sites r = E f , x s « , (2.2.16) i where Sj is the location of site i relative to the center of mass. It should be noted that the multipole expansion and interaction site approaches are by no means mutually exclusive. Models are possible where multipole expansions are conducted about several locations within a given molecule. Thus, the result will be sites which carry not only point charges, but may also carry point dipoles, quadrupoles, etc. However, such treatments are rare, and are not of relevance to this thesis. 2.3 Periodic Boundary Conditions With the recent advances in available computing power, the range of problems which are computationally accessible is increasing at a remarkable rate. Despite this, however, most computer simulation studies are limited to systems typically containing several hun-dred to a few thousand particles. While simulations of much larger systems (i. e. several million) have been carried out, these are limited to relatively simple models and usually require substantial supercomputing facilities. However, even these relatively large sys-tems are exceedingly small when compared to real-world phenomena. A simulation of one million water molecules, while an enormous task computationally, is still far from the size of even a single drop of water. The limitation of simulation studies to micro-scopic systems leads to problems associated with surface effects. Simulations of systems consisting of only several hundred molecules means that most lie in the surface region. Thus, the straightforward simulation of the bulk liquid phase is not possible. To overcome this limitation, periodic boundary conditions are employed [30,37]. This technique involves surrounding the primary simulation cell with periodic replicas of itself in all three dimensions. Obviously the chosen geometry for the primary cell must be of a Chapter 2. Simulation Methodology 22 space-filling form, and while several choices are possible (such as rhombic and truncated octahedral) the simplest choice, and that which will be employed in this thesis, is cubic. Under the periodic boundary convention, as the status of each particle in the primary cell is updated, so are all of its periodic images. Thus, if one particle leaves the primary cell through one face, its periodic image enters through the opposite face, and the total number density of the system remains constant. The fact that the system now possesses no effective surface means that all particles are in "bulk" surroundings, and the problem of surface effects has been removed. The interactions are typically calculated subject to a "minimum image" convention, by which a given molecule interacts with the closest replica of each of the other molecules. Care must be taken, however, that the periodicity of the system is not affecting the physical behaviour of the system itself. For instance, many phenomena such as phase transitions exhibit density fluctuations on characteristic lengthscales. Due to the periodicity of the simulation, any such fluctuations which are of lengthscales longer than the primary simulation cell itself will be suppressed. However, the work presented here is conducted far from any such phase transitions, and as such should not depend greatly on the particular size of the simulation cell, provided a reasonable choice is made. A further problem arises if the form of the potential being used in the simulation possesses long-range contributions, as might be the case for molecules containing dipole moments or point charges. Under the minimum image convention, each molecule will only interact with images of other molecules which lie within a distance less than or equal to the box length being used. Any contributions from the bulk liquid beyond this distance are ignored. While this is a negligible contribution when considering short-range potentials such as the L J form, it can be very significant for the longer ranged electrostatic interactions. Consequently, an alternative to the minimum image convention, which will account for interactions over many periodic images, is required. Chapter 2. Simulation Methodology 23 The simplest form for such a method would be to sum over all replicas of the simula-tion cell out to a given distance. If we assume point charges for simplicity, this leads to the following form for the total configurational potential energy [30] tf = ££ ] r^b> (2-3.1) where n is the lattice vector for a given replica of the primary cell. The term i = j is omitted when n = 0, so as to eliminate the interaction of a particle with itself. Unfortunately, Eq. (2.3.1) is only conditionally convergent for long-range potentials such as those arising from charges and dipoles. That is, the result of the summation depends on the order in which the sum over n is performed. The most obvious choice is to systematically perform the summation over cells within a certain distance of the primary cell. Thus, more and more "layers" of cells are added on, similar to the layers of an onion. The standard technique for performing this summation is the Ewald method. A mathematical derivation of the Ewald expression is outlined in Appendix B, while only a physical interpretation will be presented here. Suppose one wishes to apply the Ewald summation method to a particular configura-tion of point charges. Imagine that each point charge, qi: is surrounded by a spherically symmetric charge distribution, p(s), of equal magnitude, but opposite sign. While any reasonable form for this distribution is possible, it is most often taken to be Gaussian p(s) = qia3exp-a2s2/^2, (2.3.2) where a determines the width of the distribution and s is the distance from the original point charge. The result is a "screening" of the interaction between the original point charges, leading to a much shorter-ranged potential. The parameter a, often referred to as the "convergence parameter", can be used to change the effective range of the potential by varying the width of the charge distribution. An increase in a results in a much sharper distribution, and thus a shorter-ranged interaction. Chapter 2. Simulation Methodology 24 In order to regain the potential due to only the original point charges themselves, a second, cancelling distribution must then be included. This distribution is also of the Gaussian form, Eq. (2.3.2), but with the sign of the original point charge. The contribution to the potential due to this cancelling distribution is calculated in Fourier space. The result of this process yields the Ewald expression for a periodic lattice of point charges [30] Z i=l j=l |n|=0 l r * j ^ "I + - 7 3 J2mTrexp(-k2/4a2)cos(k-rij) - < 2- 3- 3> 2 = 1 where erfc is the complementary error function (see Appendix B), k = 2irn/L2 are the reciprocal vectors in the lattice summation, and L is the cell length. The first term in Eq. (2.3.3) is the interaction of the screened point charges. The number of lattice vectors, n, which must be included depends on the width of the distribution used. Typically, a value of a is chosen so that the minimum image convention may once again be used. That is, so that the sum over n need only consider the contribution from n = 0. The second term in Eq. (2.3.3) is due to the Fourier space summation over the cancelling charge distributions. In this case, the sum must be performed over enough reciprocal vectors, k, to ensure that the proper, converged total energy is obtained. The third term in Eq. (2.3.3) is a correction term. In the above derivation, the calculation of the potential at a particular lattice point has included a contribution from the charge distributed about that point. This contribution must thus be removed via the third term in Eq. (2.3.3). A n alternative means of incorporating long-range potentials into systems with peri-odic boundary conditions is the reaction field method [30,36]. In this approach, pairwise interactions are only explicitly calculated when the intermolecular separation is less than Chapter 2. Simulation Methodology 25 a prescribed cut-off distance. Beyond this distance the liquid is approximated as a polar-izable dielectric continuum. This has considerable advantages in terms of computational cost as well as simplicity. The reaction field method is only of moderate interest to the work presented here, and a brief overview is given in Appendix C. 2.4 Molecular Dynamics As mentioned above, the work presented in this thesis relies exclusively on the molecular dynamics approach. This method involves the numerical solution of the equations of motion for a system of N particles, consisting of one or more sites which interact via a pairwise additive potential. If Uijfaj) is the potential between sites i and j , separated by the vector r^, then the total force on site i, f;, due to all other sites is given by fi = E~V uO-( r«)- (2-4-1) As mentioned previously, the total force on a molecule j, F,-, can then be determined, and is merely the vector sum of the forces on its constituent sites, sites Fj = E (2-4.2) i Similarly the total torque on the molecule, Tj, is given by sites Tj = E ft x si5 (2.4.3) i where Si is the position of site i relative to the molecular center of mass. This expression assumes that the individual sites themselves experience no torque. While this is true of sites carrying only point charges, it will not be true of sites possessing higher order moments. For the purposes of the work in this thesis, however, the only molecules possessing such moments will be the Stockmayer particles, which are single site models. Thus the force and torque on the particle is simply that of the single interaction site. Chapter 2. Simulation Methodology 26 2.4.1 Molecu la r Orienta t ion When dealing with nonspherically symmetric potentials such as those of interest here, a systematic method of describing the orientation of a given particle must be found. Historically, the most common approach is through the use of the Euler angles, (</>, 6, ip). These angles describe three successive rotations of an object in three-dimensional space. Let (xo,yo,zo) represent an orthogonal body-fixed axis system. The first rotation is by an angle, </>, about the z0 axis to the new set of axes, (xi, yi, z0). The second rotation is about the x\ axis by an angle 9 to the new axes (x\, 2/2, £/)• The final rotation is about the Zf axis by an angle ip to the final set of axes, (xf, yj, Zf). A graphical representation of these rotations is shown in Fig. 2.2. The rotation matrix, A, which transforms from space-fixed to body-fixed coordinates subject to the three Euler rotations is given by [30] ^ cos 4> cos ip — sin r/> cos 9 sin ip sin <p cos ip + cos (p cos 9 sin ip sin 9 sin ip ^ — cos (p sin ip — sin <p cos 9 cos ip — sin 0 sin ip + cos (p cos 9 cos ip sin 9 cos ip ^ sin 4> sin 9 — cos (p sin 9 cos 9 J (2.4.4) While the Euler representation of the orientation of a three-dimensional body appears to be sufficient, there is a problem inherent in this formulation. Consider the case where 9 — 0. This results in the other two rotations becoming degenerate, as they are now both about the original z axis. This leads to singularities in the orientational equations of motion which, while obviously not representing any real physical effect in the system, make the numerical solution unstable. To overcome this problem, an alternative choice for the representation of molecular orientations is desirable. One such approach due to Evans [38] has proven to be very successful. Due to the singularities inherent in the equations of motion based on three variables, a system of four "quaternion" parameters is employed instead. The relationship between these parameters Chapter 2. Simulation Methodology Figure 2 . 2 : The effect of the Euler angle rotations on a water-like molecule. Chapter 2. Simulation Methodology 28 and the more familiar Euler angles is given by (2.4.5) g0 = cos \9 cos \((/> 4- ip) qi = sin |r?cos\{(f) - ip) q2 — sin ^ 6s'm^((f) - ip) q3 — cos |#sin | (0 + ip). The use of four parameters to describe a system with three degrees of freedom means that one of the variables is redundant, and the representation is thus not unique. Uniqueness is obtained by imposing the further condition q2o + Ql + Ql + Ql = i- (2.4.6) The rotation matrix, A, in terms of these quaternion parameters now becomes ( ol + q\-q\-q\ 2(qiq2 + q0q3) 2(gig3 - q0q2) ^ A= 2(q1q2-q0q3) ql - q\ + q\ - q\ 2(q2q3 + q0Qi) • (2.4.7) ^ 2(q1q3 + q0q2) 2(q2qz - q0qx) ql - q2 + q\ + q\ j This approach provides a singularity-free method of representing molecular orientations. In fact, the quaternion method has become the preferred method of dealing with orien-tation in molecular dynamics simulations. 2.4.2 Time Integration The classical equations of motion for a system of N particles consists of a series of ordinary differential equations which can be solved by the standard finite difference approach. If the positions, orientations, and velocities at time t are known, one can attempt to accurately determine these quantities at time t + 6t. Countless methods are available [30] for the numerical solution of these equations, but only two are of interest here. Chapter 2. Simulation Methodology 29 For systems with continuous potentials, one can reliably predict the positions, orien-tations, velocities, etc. at time t + 8t by employing a Taylor expansion about time t. For instance, the predicted position of particle i at time t + 5t, rf(t + 5t), would be given by r?(i + St) = n(t) + 6tti{t) + \sevi + (2.4.8) where the dots indicate differentiation with respect to t. Dynamical information, such as the acceleration of the particle, af, can also be predicted a?(t + St) = ai(t) + Stai(t) + ^St2ai + (2.4.9) The predicted quantities at time t+St are based solely on their values and time derivatives at time t. No information about intermolecular forces or the equations of motion has been explicitly considered. The new position is predicted based solely on the "history" of the particle's current trajectory. So long as the particles interact via continuous potentials, this yields a moderately accurate prediction. A correction step can then be applied which explicitly takes into account the intermolecular forces. This involves determining the particle's acceleration due to these forces, and comparing with the predicted acceleration. This yields an estimate of the required correction Aa;(i + St) = alit + St) - af(t + St). (2.4.10) This correction can then be applied to the remaining quantities r?(t + (Jt)=r?(< + <Jt)-l-ciAai(t + <5t) (2.4.11) aci(t + St) = api(t + St) + c 2Aa;(£ + St), (2.4.12) where c\ and c 2 are constants. The optimal choice for these constants depends on the number of time derivatives included in the prediction step. This algorithm was developed by Gear [30], and is generally known as the Gear predictor-corrector. In theory, the Chapter 2. Simulation Methodology 30 prediction and correction step should be iterated until complete convergence has been achieved. However, this is extremely costly in terms of computational time, as it must be carried out for each particle at each timestep, and typically only one prediction and correction cycle is used. Provided that the initial prediction is sufficiently accurate, and the particles interact through well-behaved potentials, this has proven to be a very robust method, particularly for models such as Stockmayer particles. An alternative to the Gear predictor-corrector is the Verlet algorithm [30]. Once the acceleration of a particle at time t, a;(£), has been obtained from the equations of motion, the updated velocity at time t+^5t, V j ( £ + |6"£), can be determined via the simple relation V i ( * + = Vi(t - ^6t) + Sta^t). (2.4.13) This updated value of the velocity can then be used to update the particle's position at time t + 5t n(t + St) = n(t) + 5tVi(t + l-5t). (2.4.14) As is evident from Eq. (2.4.13), the velocities are stored at the mid-steps, t — \5t and t + \8t instead of at t and t + 5t. The algorithm, originally developed by Verlet, is typically referred to as the "leap-frog" approach for this reason, and has proven to be somewhat more stable for interaction site models than has the predictor-corrector. Since any calculation of the temperature or Hamiltonian of the system requires that potential energy and kinetic energy both be determined at the same instant of time, the velocity at time t, when required, can be calculated as the simple average of the values at t — ^5t and t + \5t Vi(t) = \{vi{t + \st) + Vi(t - l-5t)). (2.4.15) Chapter 2. Simulation Methodology 31 2.4.3 The Equations of Motion The preceding discussion of time integration has implicitly assumed that some means was available of obtaining the particle accelerations, ^(t). While the force, Fj(t), may be obtained directly from the intermolecular potentials, the relationship between the force and acceleration can take several forms. The most obvious choice would be that corresponding to Newtonian mechanics F, = m ^ , (2.4.16) where mi is the particle mass. As Newtonian mechanics conserve energy, a system evolv-ing in time subject to these particular equations of motion would do so at constant E. As Newton's equations are what govern the trajectories of real objects within the classical limit, the use of anything but the Newtonian form in a classical simulation may seem somewhat unjustified. Nonetheless, there are many occasions for which other forms for the equations of motion are desirable. The most common such instance, and that of interest here, occurs in the study of constant temperature, T, systems. Suppose, for example, one wishes to study a system under the influence of a temporally varying external field. In this case, although the particle trajectories themselves may be evolving subject to Newtonian mechanics, the presence of the external field acts as a continuous source of energy for the system. This will result in a constant increase in the total energy of the system as the applied field does work upon it. Thus, the temperature of the system will increase as well, as the periodic boundary conditions which were required to remove surface effects have also removed any surface boundaries through which the excess energy would normally escape. In fact, this proves to be a problem even in equilibrium systems with no such external source of energy. Due to the slight numerical inaccuracies inherent in any time integration technique, the energy of the system is not rigorously conserved. As a result, even the temperature of a closed system may drift Chapter 2. Simulation Methodology 32 substantially over time if no effort is made to constrain it. Thus, Newtonian mechanics is insufficient, and a slight modification of the equations of motion is required. The simplest approach to maintaining a constant temperature is that of velocity rescaling. The instantaneous' temperature, T, of a system can be calculated from the individual particle velocities as 2 N 1 T = -i—T^-m^il2. (2.4.17) 3Nkb V 2 1 The system can then be rescaled to the desired temperature, T", simply by multiplying all particle velocities by the factor yjT'/T. This technique allows one to periodically correct for any drift in the overall temperature due to numerical inaccuracies, while maintaining Newtonian trajectories between rescalings. While this approach is sufficient to compensate for such drifts in temperature, a somewhat more elegant approach may be desired. An alternative method involves the formulation of equations of motion which conserve temperature, rather than energy. In this case, Eq. (2.4.16) is replace by [30] m i a i = F i-e(r ,v)v i , (2.4.18) where £(r, v) "adjusts" the acceleration so as to maintain a constant temperature. The choice for £ which minimizes the difference between these new equations of motion and the original Newtonian trajectory is given by [30] £( r ,v) = ^7'*\ (2.4.19) Thus, the use of Eqs. (2.4.18) and (2.4.19) allows for the sampling of the canonical ensem-ble while still generating trajectories which differ only marginally from the Newtonian ones. However, it is not clear what effect such non-Newtonian equations would have on nonequilibrium behaviour, and thus the constant temperature formulation is best Chapter 2. Simulation Methodology 33 restricted to equilibrium calculations only. Similar techniques are also available for con-ducting constant pressure simulations, but are not of significance to the work presented in this thesis. Chapter 3 Solvation Dynamics in Simple Solvents The majority of computer simulation studies conducted to date have involved the study of relatively simple models. The reason for conducting research on what might seem to be somewhat unrealistic systems is two-fold. In the study of liquid phase phenomena, as in all studies, it proves highly instructive to begin with a relatively simple model. This allows one to freely investigate the effects of various properties such as molecular mass or polarity without the unnecessary complications which would arise should one choose to model a more realistic system. Once the baseline behaviour for these simple models has been investigated, one may then move on to more sophisticated studies, the goal of which is to explore how the observed behaviour applies to "real-world" systems. The other motivation for studying these simple solvent models is to better facilitate a comparison with theoretical treatments which tend to be limited to comparatively simple systems. Results generated via computer simulation serve not only as a means of studying the robustness of the theoretical treatment, but comparison with theory may also shed light on the simulation results. The solvent models presented in this chapter have been chosen with the above two motivations in mind. 3.1 Model and Computational Details Our system of interest consists of a solute particle, s, with position r s , immersed in a binary mixture of nonpolarizable dipolar solvent particles with positions r*. The two solvent species will be labelled Solvent W for the more weakly polar solvent, and Solvent 3 4 Chapter 3. Solvation Dynamics in Simple Solvents 35 S for the more strongly polar. These solvent particles will be characterized by the dipole moments fiw and p,$- The solute particle, while initially neutral, carries a point charge q following ionization. The total configurational energy, U, of the system is given by U=W ug*(ry) + £ ufsR(vis) + \ Y: t i g f l ( r y , to) + £ < s ° ^ s , /O , (3-1.1) ij i ij i where = r; — r,- and ris = — r s . In all calculations both the short-range (SR) solvent-solvent energy, ufR, and the short-range solvent-solute energy, ufR, were taken to be of the usual Lennard-Jones potential given by Eq. (2.2.2) The remaining terms, ufjD and are the dipole-dipole (DD) and ion-dipole (ID) interactions, respectively. Except where otherwise stated, the long-range dipole-dipole and ion-dipole interac-tions were accounted for using the Ewald summation method as described in Section 2.3 and Appendix B. The solute particle, when ionized, is considered to be at infinite dilu-tion and thus is only present in the primary simulation cell. Consequently, the ion-dipole energy is identical to the energy of insertion of a single point charge into an irregular dipole lattice and can be expressed in the form [37] * = - | t T « < M . • - rf* • V ) 2*£Lpd) , (3.L2) where V is the volume of the simulation cell, a is the convergence parameter which partitions the energy between the real-space and A;-space components, and erfc(x) is the complementary error function. This result differs from the case of finite ion concentration only by a constant term. However, it should be pointed out that in the case of finite concentration the system as a whole would carry an infinite charge in the absence of counterions and thus the Ewald sum would be divergent. For our purposes though, the equation is perfectly valid as the system carries only a net charge of q. In the present calculations a value of a = 6.4 was used and all lattice vectors, k, satisfying k2 < 62 were included in the sum. This has been shown to yield satisfactory energy calculations [37]. Chapter 3. Solvation Dynamics in Simple Solvents 36 In all calculations, the Lennard-Jones parameters o and e as well as the mass, m, were taken to be the same for all species. The system is characterized by specifying the reduced dipole moments, p* = /^/(ea 3) 1/ 2 , the reduced temperature, T* = kbT/e, the reduced moment of inertia, I* = I/(mo2), the total reduced density, p* = po3, and the solvent mole fractions, Xj. The total density, p, includes all three species. All systems were studied at p* = 0.8 and T* = 1.35. The moment of inertia for both solvent species was fixed at I* = 0.025 and the reduced ionic charge on the solute was chosen to be q* = (g 2/ae) 1/ 2 = 12.0. These reduced parameters lie in the physically realistic range. The time integration was carried out via a fourth order predictor-corrector algorithm as described in Section 2.4.2 using a reduced timestep of t* = (e/mo2)ll2t = 0.0025 with rotational motion being accomplished through the use of quaternions [30]. The systems were initially equilibrated at constant temperature for 150000 timesteps, with averages accumulated for a further 50000 timesteps. A total of 20-40 Configurations were extracted from the equilibrium trajectory with the uncharged solute particle, and were subsequently used as initial conditions for the nonequilibrium simulations. All MD calculations following ionization were carried out at constant energy. The solute is charged at t — 0 and progress to equilibrium can be monitored via the solvation response function S(t) (3.1.3) (Uss(0) - CL(oo)) ' and particle solvation response function (3.1.4) discussed previously. Chapter 3. Solvation Dynamics in Simple Solvents 37 3.2 Results and Discussion Molecular dynamics calculations were performed on nine primary systems labelled 0 through 8 which are summarized in Table 3.1. Systems 4 and 7 have identical solvent characteristics and mole fractions, but System 4 contains a total of 256 particles while System 7 contains 500. The solvation response functions for each are shown in Fig. 3.1, along with error estimates for three points on each curve. These error estimates represent one standard deviation as calculated from 20 independent trajectories. It can be seen that the discrepancies between the two simulations are small, and well within the limits of uncertainty. This indicates that 256 particles are sufficient for the dynamics under consideration here. Consequently, all remaining systems in Table 3.1 were simulated using a total of 256 particles. The possible influence of boundary conditions is also of concern. System 6 is identical to System 5 (see Table 3.1) with the exception that reaction field boundary conditions (see Appendix C) were employed. The solvation response functions for both systems are shown in Fig. 3.2, and again the discrepancies which occur are well within the uncertainty limits. This strongly suggests that our choice of the Ewald summation method in accounting for the long-range forces is not having any serious unphysical effects on the results. Systems 1 and 2 differ only in composition and their solvation response functions are shown in Fig. 3.3, along with that of System 0 for comparison purposes. System 0 provides results typical of pure solvents. The curves, as with all that will be presented in this chapter, have been plotted on a timescale determined by choosing the reduced parameters to correspond to those of argon [30] (e = 119.8 K and o = 0.341 nm), resulting in a timestep of 5.4 fs. Somewhat surprisingly, the solvation response functions for Systems 1 and 2 differ only little from that typical of pure solvents. The response functions are still dominated by the initial reorientational motion which accounts for the Figure 3.1: Solvation response functions for systems of different size. The solid curve is System 4 (256 particles) and the dotted curve is System 7 (500 particles). Figure 3.2: Solvation response functions for different boundary conditions. The solid curve is System 5 (Ewald method) and the dotted curve is System 6 (reaction field method). Chapter 3. Solvation Dynamics in Simple Solvents 40 0.8 h 0.6 \-S(t) 0.4 "i i i i i i i i i i i i i r i i i l i i i i l i 0.2 0 - 0 . 2 • 1 1 1 1 1 » 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 o 5 10 15 20 25 Time (ps) Figure 3.3: Solvation response functions for System 1 (dotted) and System 2 (dashed) along with that of System 0 (solid) for comparison. Chapter 3. Solvation Dynamics in Simple Solvents 41 Table 3.1: A summary of the systems simulated. The reduced dipole moments and mole fractions are given. Note that p* = 0.8 and T* = 1.35 for all systems studied. System Pvv Xyy Xs 0 2.00 - 1.00 0.00 1 1.00 2.00 0.75 0.25 2 1.00 2.00 0.95 0.05 3 0.50 2.00 0.75 0.25 4 0.50 2.00 0.90 0.10 5 0.50 2.00 0.95 0.05 6 0.50 2.00 0.95 0.05 7 0.50 2.00 0.90 0.10 8 1.00 3.00 0.95 0.05 vast majority of the relaxation energy. However, the subsequent decay to equilibrium following reorientation now occurs on a somewhat longer timescale than in pure solvents. This can be attributed to the spatial redistribution of the solvent species which must occur to achieve equilibrium. Evidence of this behaviour can be clearly seen by examining the time evolution of the number of solvent particles of each type within the first solvation shell, as shown for System 1 in Fig. 3.4. Although the solvation response function serves as a measure of the relaxation of the total system, it has been shown previously [10] that its behaviour is dominated by effects in the immediate environment of the solute particle, and thus we focus our attention there. As is evident in Fig. 3.4, the immediate environment of the solute particle is indeed becoming progressively richer in Solvent S as time progresses, with the concentration of Solvent W decreasing in a similar fashion. Interestingly, the obvious correspondence between the two curves is a result of the fact that the total num-ber of solvent particles in the first solvation shell is effectively constant throughout most of the equilibration process. The implications of this and the nature of the redistribution will be discussed later in this chapter. Chapter 3. Solvation Dynamics in Simple Solvents 42 <Nj(t)> 10 8 6 4 2 Total W 0.8 0.6 S(t) 0.4 0.2 0 -0.2 i i i i i i i i i i i i i i i i i i i i i i i i i i • ' i » i i i i ' i ' i i i i i i i i i i i i i I d 0 5 10 15 20 25 Time (ps) Figure 3.4: The number of solvent particles of each type (labelled W and S) together with the total number present in the first solvation shell for System 1. The corresponding solvation response function is also shown. Chapter 3. Solvation Dynamics in Simple Solvents 43 While the presence of the predicted solvent redistribution is clearly evident from Fig. 3.4, it appears to have relatively little effect on S(t). This can most likely be attributed to the dominance of the energy relaxation by the initial reorientational motion. The ion-solvent energy resulting from the subsequent redistribution of the solvent species ac-counts for a comparatively small fraction of the total. As a result, the solvation response functions appear remarkably similar to those obtained for pure solvents. With this in view, it should be possible to choose the various system and solvent parameters such that the relative importance of the reorientational motion is substantially reduced. The most obvious way to accomplish this is to reduce the dipole moment of the less polar species. The result should be a decrease in the total energy obtained from the reorien-tational motion of the solvent, and a corresponding increase in the relative importance of the solvent redistribution of interest. This change in relative polarities should also have the less obvious effect of increasing the discrepancy between the initial and final solvated states. As the two solvent species become more dissimilar, the observed degree of preferential solvation for both the ionized and unionized solute particles should be correspondingly more prominent. These two factors should also serve to increase the dependence of the rate of relaxation on concentration. It was with these possible effects in mind that Systems 3, 4, and 5 were considered. The solvation response functions for Systems 3, 4, and 5 are shown in Fig. 3.5. These systems are similar to 1 and 2 except that the dipole moment of Solvent W is reduced to = 0.50. In contrast to the results discussed above, S(t) now shows a substantial deviation from the behaviour typical of pure solvents. As expected, the significance of the reorientational motion of the solvent, although still important, has been substantially reduced. The remainder of the solvation energy then relaxes through the longer timescale behaviour which was only partially evident in the solvation response functions for Systems 1 and 2. Again, this can be associated with the redistribution of the solvent species with Chapter 3. Solvation Dynamics in Simple Solvents 44 Table 3.2: The composition of the first solvation shell at equilibrium. TV- and A/" (j = W, S) are the average numbers of Solvent W and Solvent S particles present within a radius of 1.3<r for both ionized (i) and unionized (u) solute particles. The error limits represent one standard deviation. System iVfr 7V| Nls 1 6.3 ± 0.3 1.4 ± 0.3 3.5 ± 0.3 4.7 ± 0.3 2 7.6 ± 0.2 0.3 ± 0.1 6.0 ± 0.3 2.6 ± 0.4 3 6.1 ± 0.1 1.6 ± 0.2 2.6 ± 0.2 5.7 ± 0.2 4 7.5 ± 0.1 0.4 ± 0.1 3.8 ± 0.1 4.6 ± 0.2 5 7.5 ± 0.1 0.4 ± 0.1 4.1 ± 0.1 4.2 ± 0.2 respect to the solute particle that is best illustrated by considering (Nj(t)), as shown for System 5 in Fig. 3.6. As in the earlier example, the increase in abundance of Solvent S particles in the first solvation shell is matched by a corresponding decrease in the abundance of Solvent W, resulting in a relatively constant total number. As shown in Fig. 3.5, the relaxation behaviour also shows a significant dependence on system composition. This is to be expected, as a system with a lower concentration of the more polar species will presumably obtain even less energy from solvent reorientation, thus increasing the importance of the solvent redistribution evident in Fig. 3.6. The degree of preferential solvation which does occur can be seen by examining the equilibrium compositions of the first solvation shell for both the ionized and unionized solute, and some results are summarized in Table 3.2. The change in composition of the first shell when the solute is ionized reflects the increased affinity for Solvent S. The preferential solvation can also be clearly seen by examining the solute-solvent pair distribution functions before and after ionization. Fig. 3.7 shows the initial and final pair distribution functions for System 3. As expected, the neutral solute maintains solvation shells of approximately the same composition as the bulk liquid, with a slight preference for the less polar species. Indeed, this slight preferential solvation of a neutral solute Chapter 3. Solvation Dynamics in Simple Solvents 45 Figure 3.5: Solvation response functions for Systems 3, 4, and 5 along with that of System 0 for reference. Chapter 3. Solvation Dynamics in Simple Solvents 46 0 5 10 15 20 25 Time (ps) Figure 3.6: The number of solvent particles of each type (labelled W and S) together with the total number present in the first solvation shell for System 5. The corresponding solvation response function is also shown. Chapter 3. Solvation Dynamics in Simple Solvents 47 Figure 3.7: Pair distribution functions for System 3 with unionized (a) and ionized (b) solute. The curves are for Solvent W (solid) and Solvent S (dotted). Chapter 3. Solvation Dynamics in Simple Solvents 48 is to be expected when the solvent-solvent interactions are considered. The more polar particles are more likely to be found in the bulk where they can be completely solvated by other polar species, thus leaving the less polar particles to solvate the nonpolar solute. In contrast, the ionized solute exhibits a very high degree of preferential solvation by the more polar species, as indicated by the very high initial peak in the appropriate distribution function. It is this relative difference in solvation states that determines what portion of the relaxation energy will be obtained by the initial reorientation, and what portion will result from the subsequent solvent redistribution. 3.3 A Simple Phenomenological Model To better understand the underlying features of the solvent relaxation, it is useful to con-sider a simple phenomenological model. For instance, if one considers the first solvation shell of the solute particle, one can regard the individual solvent particles as occupying small regions, or "sites", within this volume. These sites would obviously have no phys-ical interpretation and simply provide a useful way of representing the solvation state of the solute. Such a site could be occupied by either a Solvent W particle, in which case it will be labelled it IW, or a Solvent S particle, in which case it will be labelled it IS. Additionally, the inclusion of "unoccupied sites", / , would allow for changes in the total coordination number of the solute. If one makes the assumption that the probability of occupation of a given site is independent of the occupancy of the remaining sites, one can formulate the simple phenomenological equations I + S - I S (3.3.1a) / + W - IW , (3.3.16) Chapter 3. Solvation Dynamics in Simple Solvents 49 where the k{ represent rate constants for the various processes. This particular model is somewhat analogous to certain chemical processes, such as ligand binding to a surface on which there are numerous potential binding sites. In such a case, the presence of two or more potential ligands in solution will yield competitive binding processes with the equilibrium surface coverage dependent on the relative affinities of the two species for the binding sites. In the case of ion solvation of interest here, one can regard the surface as having been replaced by the solute particle, with the sites now representing regions within the first solvation shell that can be occupied by competing solvent species. One can then define [IS] and [IW] as the fractions of these sites occupied by each species, and [I] as the fraction of unoccupied sites such that [/] + [IS] + [IW] = 1 . (3.3.2) The rate equations for [IW] and [IS] can then be written as d[IW] dt d[IS] = k'3[I] - k4[IW] , (3.3.3a) = k[[I] - k2[IS] , (3.3.36) dt where it has been assumed that the bulk concentrations of Solvent W, [W], and Solvent S, [S], are constant, and the pseudo rate constants k[ = ki[S] and k'3 = k3[W] have been defined. Under these conditions, the rate equations given in Eqs. (3.3.3) have the solutions [IW] - [iW]oo = ae~t/Tl + be~t/T2 , (3.3.4a) [IS] - [/5]oo = ce-t/Tl + de~t/T2 , (3.3.46) where a, b, c, and d are constants, [JWjoo and [IS]^ are the final equilibrium values, and Tip are given by = i [ ( ^ + fc2 + ^ -f- fc4) ± [(^ + /c2 + ^ + A;4)2 - 4(A;'1 + A;2)(A;^ + fc4) + 4A;2A;4]^ ] • (3.3.5) Li Chapter 3. Solvation Dynamics in Simple Solvents 5 0 Also of interest is the total fraction of occupied sites, [O], in the first solvation shell and this is given by [O] - [O]oo = (a + c)e-*/Tl + [b + d)e~t^ . (3.3.6) A minor problem arises in that the simulation results describe the first shell in terms of the number of particles present, while Eqs. (3.3.4) and (3.3.6) pertain to the fraction of occupied sites. Thus, to compare directly with simulation results it would be necessary to specify the total number of available sites. This problem can be circumvented by consid-ering particle solvation response functions as defined earlier. From the phenomenological theory one obtains PW{T)=[BO-BL=CWE~T,TI+(L - CW)E~T/T2• (3-3-7A) Ps(t) = Cse^ + (1 - C 5 ) e - t / T 2 , (3.3.76) P(t) = Ce~t/Tl + (1 - G > - ' / T 2 , (3.3.7c) where [IW]0 is the initial value of [IW], Cw = a/(a + b), Cs = c/(c + d) and C — (a + c)/(a + b + c + d). Also, the subscripts W and S refer to Solvents W and S, and P(t) is the particle solvation response function for the total number of particles. Equation (3.3.7c) is obtained using [O] as given by Eq. (3.3.6). It should be noted that if only one timescale occurs in the relaxation (i.e., TX = r 2) then all three expressions reduce to a single exponential. In principle, it is possible to fit Eqs. (3.3.7) to the simulation results using the five parameters Cw, Cs, C, T\ and r 2 . However, a satisfactory fitting requires that all three curves be considered simultaneously. While this can be attempted and is briefly discussed below, it is more instructive to further simplify the equations on physical grounds. From the M D results one can make several observations. Initially, there is a small rapid increase in the total number of particles, (N(i)), in the first solvation shell. This Chapter 3. Solvation Dynamics in Simple Solvents 51 is essentially an electrostriction effect arising from the strong charge-dipole interactions when the ion is "switched on". In other words, when the ion is initially charged solvent particles of both types in its immediate vicinity are quickly "pulled" closer to the ion. Following this initial period, (N(t)) remains effectively constant while Solvent S and Sol-vent W redistribute much more slowly until equilibrium is reached. The redistribution is symmetric with the increase in (Ns(t)) with time being mirrored by a corresponding de-crease in (Nw(t)). These observations are particularly true for systems dilute in Solvent S. Also, they are very similar to the steady-state behaviour observed in some sequen-tial chemical reactions where the concentration of an intermediate remains effectively constant following a brief initialization period. Such behaviour is observed in systems where the intermediate, in our case an unoccupied site I, is formed slowly but rapidly consumed. Taking account of these observations, it is assumed that selective ion solvation occurs in two distinct steps, electrostriction and redistribution, and that two distinct timescales are involved. These are characterized by an electrostriction (el) relaxation time, r e i , which describes the rate of change of the total number of particles, (N(t)), in the first solvation shell, and a redistribution (re) relaxation time, r r e , which describes the rate of compositional change in the first shell. Further, it is assumed that r r e 3> r e i , such that that (N(t)) is constant after a short initial period. This is a steady-state approximation and implies that the redistribution of solvent particles is completely symmetric. Return-ing to Eqs. (3.3.7), we now identify T\ as re\ and r 2 as r r e (note that this identification is arbitrary since Eqs. (3.3.7) are symmetric in T\ and r 2 ) . Our steady-state assumption implies that b — —d and hence C = 1. Eqs. (3.3.7) then simplify to the form Pw(t) = Cwe-t/T»l + (1 - Cw)e-t/T™ , (3.3.8a) Ps{t) = Cse'tlTA + (1 - Cs)e't^e , (3.3.86) Chapter 3. Solvation Dynamics in Simple Solvents 52 Table 3.3: Relaxation times for selected systems obtained from fitting Eqs. (3.3.8) to simulation results. System 2 0.774 ps 77.1 ps 3 5.22 ps 32.8 ps 5 1.04 ps 29.4 ps P(t) = e~t/Tei . (3.3.8c) Assuming that Eqs. (3.3.8) apply, re\ can be obtained by fitting P(t). The remaining three parameters, r r e , Cyv and Cs, are then found by fitting Pw(t) and Ps(t) simultane-ously. The fitted curves as well as the M D results for Systems 2, 3 and 5 are shown in Figs. 3.8-3.10 and the relaxation times obtained are given in Table 3.3. The simulation results used in the fits included 40 trajectories to reduce the noise. We see from Figs. 3.8-3.10 that Eqs. (3.3.8) fit all three systems very well. It is particularly gratifying to note that for System 2 [Fig. 3.8] and System 5 [Fig. 3.10] the initial increase in Pw(t) due to electrostriction is completely captured by the simple model. From Table 3.3, we see that for Systems 2 and 5 which are dilute in Solvent S (Xs = 0.05), the relaxation times obtained are well separated and r r e S> r e i in accord with our assumptions. However, for System 3 where Xs = 0.25, r r e and r ei.are not as well separated, but excellent fits [Fig. 3.10] are still obtained for all three particle solvation response functions. A possible explanation is as follows. In discussing Eqs. (3.3.7) above, it was noted that if there is only one timescale in the process, then Pw(t), Ps(t) and P(t) become identical and decay as single exponentials. This limit remains satisfied by Eqs. (3.3.8) if r r e = re\, and likely this partially accounts for the good fits found for System 3 despite the fact that the distinct timescales assumption made in obtaining Eqs. (3.3.8) is less well satisfied. Five-parameter simultaneous fits of all three response functions to Eqs. (3.3.7) were Figure 3.8: Particle solvation response functions obtained from simulation (dotted) and from fitting of Eqs. (3.3.8) (solid) for System 2. Chapter 3. Solvation Dynamics in Simple Solvents 54 t (ps) Figure 3.9: Particle solvation response functions obtained from simulation (dotted) and from fitting of Eqs. (3.3.8) (solid) for System 3. Chapter 3. Solvation Dynamics in Simple Solvents 55 Figure 3.10: Particle solvation response functions obtained from simulation (dotted) and from fitting of Eqs. (3.3.8) (solid) for System 5. Chapter 3. Solvation Dynamics in Simple Solvents 56 also attempted. This proved to be a somewhat tedious undertaking as the results ex-hibited a slight sensitivity to the choice of initial parameter values. Nonetheless, the relaxation times obtained from these fits were in reasonably good agreement with those given by Eqs. (3.3.8), providing further support for the underlying physical picture dis-cussed above. 3.4 Further Results and Remarks on Mechanism The previous discussion regarding the phenomenological model was restricted to an ex-amination of the first solvation shell of the solute as a function of time, and resolved the behaviour into electrostriction and redistribution timescales. While the nature of the electrostriction is relatively obvious, it is not clear what determines the rate of the solvent redistribution. At sufficiently low concentrations one would expect that the number of Solvent S particles present in the vicinity of the solute at the moment of ionization to be extremely small. Consequently, the rate of solvation would be closely related to the rate at which this species could diffuse through the bulk liquid. More accurately, it would de-pend on the ion-solvent mutual diffusion rate subject to the long-range ion-dipole force. In contrast, a system with a high concentration of Solvent S would possess an abundance of the more polar species throughout the solvation shells of the solute particle. In fact, at sufficiently high concentrations it might be possible for the first shell to reach near equilibrium solely via exchanges with the second shell. A further redistribution of sol-vent species throughout the remaining shells would be required to reach true equilibrium, but the overall rate of relaxation of the first shell would be largely determined by the exchange process. Finally, at more moderate concentrations, one might expect to find systems where both processes contribute to the rate of solvent redistribution. While distinguishing between these cases proves quite difficult, it is possible to devise Chapter 3. Solvation Dynamics in Simple Solvents 57 simulations which are instructive. For instance, if we examine the relaxation times for Systems 2 and 5 as listed in Table 3.3 we see that the rate of solvent redistribution decreases sharply when the dipole moment of the weakly polar species is increased from 0.5 to 1.0. This decreased rate can be attributed to two possible factors. As the dipole moment of Solvent W is increased, the dielectric constant of the bulk liquid would likewise increase. The result would be a reduction in the long-range solvent-averaged ion-dipole force, thus decreasing the rate at which any field dependent diffusion might occur. Also, the increased dipole-dipole interaction between Solvent S and Solvent W could slow the diffusion rate. If diffusion is the dominating factor in the relaxation mechanism, these two effects could account for the observed decrease in rate. On the other hand, the increased dipole moment will also serve to increase the ion-solvent attractive force for Solvent W, thus making it less likely to be displaced from the first solvation shell. If the exchanges are limiting the rate of redistribution, then one would also expect the observed dependence on dipole moment. To determine which of these scenarios is more likely the true cause, it is instructive to construct systems which share important characteristics of both Systems 2 and 5. For instance, we can vary the dipole moment of the less polar Solvent W particles with distance from the solute particle such that Mn) = \ + \e~^l°-^ , (3.4.1) where Vj is the distance between the solute particle and the ith Solvent W particle. This choice will yield dipole moments of approximately 1.0 when in the vicinity of the solute particle, but decaying rapidly to 0.5 further into the bulk. Such a system would possess diffusive properties similar to System 5, but with first shell effects more closely tied to those of System 2. The opposite result can be obtained by choosing /4(r0 = 1 - I e-2<"/*-i> a . (3.4.2) Chapter 3. Solvation Dynamics in Simple Solvents 58 This choice should yield a system with a rate of diffusion more closely related to that of System 2, but with the exchange process presumably closer to that of System 5. The (Ns(t)) plots obtained with both rampings are compared with those of Systems 2 and 5 in Fig. 3.11. It can be seen that both choices of ramping yield systems with first shell relaxation rates faster than that of System 2, but slower than that of System 5. The fact that an intermediate result is obtained in both cases is a strong indication that both of the processes discussed above play significant roles at the concentration considered. It is also interesting to look at the exchange process more closely. The actual mi-croscopic relationship between (Nw(t)) and (Ns(t)) is not clear from the results plotted above. While the two species appear to undergo a simultaneous one-to-one exchange process, yielding a relatively constant (N(t)) curve at longer times, the results shown have been averaged over many trajectories. Therefore, any specific relationship between the insertion of a particle of one species and the ejection of a particle of the other would be obscured by the averaging. Consequently, it is instructive to examine the composition of the shell immediately prior to and following the addition of a Solvent S particle within specific trajectories. To accomplish this, each individual trajectory was examined and the points in time at which insertion of a Solvent S particle occurred were identified. For each such occurrence the trajectory was shifted in time such that the insertion took place at t = 0, and the results of all trajectories were then averaged. It should be noted that some trajectories will be included multiple times in the averaging process, provided multiple insertions occurred. In such a case the trajectory, though present multiple times, would be shifted differently for each occurrence. The results of this analysis for Systems 3 and 5 are shown in Fig. 3.12. The nature of this calculation is such that the curve corresponding to (Ns(t)) will exhibit a jump discontinuity of unit magnitude at t = 0. It is the behaviour of (Nw(t)) with respect to this time origin that is of interest here. Examining the curves in Fig. 3.12, it is apparent Chapter 3. Solvation Dynamics in Simple Solvents 59 5 i i i i i i i i i i i i i i i i i i i i i i i i i i i 4 \ -Figure 3.11: A comparison of ramping effects on the composition of the first solvation shell as a function of time. The dotted and solid curves are for Systems 2 and 5, respectively. The dashed curve was obtained using Eq. (3.4.1) and the long-dashed curve using Eq. (3.4.2). Chapter 3. Solvation Dynamics in Simple Solvents 60 that the insertion of a Solvent S particle is .followed almost immediately by an ejection of a Solvent W particle already present in the first solvation shell. The result is a sharp peak in the total number of solvent particles, (N(i)), at t = 0. While these observations are not entirely surprising, they are significant. They indicate that the exchange process is not dependent on fluctuations which momentarily eject a particle from the first shell, thus creating a "hole" into which an insertion can occur. Rather, the attractive force between the ion and the more strongly polar solvent "pulls" a Solvent S particle into the first solvation shell. The "overcrowding" then results in an almost immediate ejection of one of the Solvent W particles, returning the total number to its value prior to the insertion. Thus, one has an almost simultaneous particle exchange with only an extremely shortlived intermediate state in which the total number deviates from its average value. Finally, it should be noted that the results discussed above focussed on the influence of the weakly polar solvent on the solvation rate. As it is the species present in the higher concentration, it is expected to strongly influence processes such as exchange and bulk diffusion. Nonetheless, one would expect the dipole moment of the more polar species to also play an important role due to the direct interaction with the ion. Therefore, an additional system, System 8, was considered in order to investigate the effect of increasing the dipole moment of Solvent S. Simulation results for this system are shown in Fig. 3.13 together with those of two other systems for reference. Somewhat surprisingly, the resulting exchange rate appears to be unchanged in moving from a system in which fi*s — 2.0 (System 2) to one in which p*s = 3.0 (System 8). Increasing the attractive interaction between the ion and Solvent S appears to have little effect. This is possibly explained by the fact that with the larger dipole moment Solvent S is more strongly solvated by Solvent W, thus decreasing its mobility. Apparently, this effect cancels the increased attraction between the ion and Solvent S. Chapter 3. Solvation Dynamics in Simple Solvents 61 -20 -10 0 10 20 t (ps) Figure 3.12: Modified shell counts as described in the text for System 3 (a) and System 5 (b). The curves are (Nw(t)) (W), (Ns{t)) (S) and (N(t)) (Total). Chapter 3. Solvation Dynamics in Simple Solvents 62 Figure 3.13: A comparison of the influence of varying solvent dipole moments. The curves are for System 5 (solid), System 2 (dotted) and System 8 (dashed). Chapter 3. Solvation Dynamics in Simple Solvents 63 3.5 Summary In this chapter, M D simulations have been used to investigate the dynamics of selective ion solvation in model binary solvents. It was found that even for systems with relatively large differences in solvent polarity (eg., a factor of two), the usual solvation response function, S(t), is not very dependent on solvent composition and remains quite similar to the results obtained for pure solvents. This can be true despite a substantial change in the composition of the first solvation shell when the solute is charged. The reason for this is that a large fraction of the total ion-solvent energy is still obtained by simple reorientation of solvent particles and only a much smaller fraction is related to the sub-sequent redistribution of solvent species. This means that for many systems S(t) will be rather insensitive to the redistribution which occurs in selective solvation. Indeed, recent experimental work by Nishiyama and Okada [39] have confirmed this relative insensitiv-ity of the solvation response function. Instead, the authors suggest that monitoring the relaxation of the spectral width, instead of the spectral maximum, appears to be a more useful approach when studying solvation dynamics in binary solvents. Of course, if the difference in solvent polarity is large enough, S(t) does show signif-icant deviations from the pure solvent result and the influence of solvent redistribution becomes evident. Under these conditions, S(t) also exhibits a strong dependence on the composition of the binary solvent. However, the relaxation rate does not have a simple dependence on the relative magnitudes of the solvent dipole moments. For example, decreasing the dipole moment of Solvent W increases the relaxation rate, whereas in-creasing that of Solvent S does not have a parallel effect. When the dipole moment of Solvent S is increased, the increase in the ion-solvent force appears to be largely offset by the slower diffusion rate due to the stronger dipole-dipole interactions between the solvent species, and hence there is little change in the relaxation rate. Chapter 3. Solvation Dynamics in Simple Solvents 64 Partially because of the insensitivity of S(t) discussed above, and partially to acquire more physical insight into the nature of the selective solvation process, particle solva-tion response functions, Pj(t), have been introduced which describe the relaxation of the composition of the first solvation shell. It has also been shown that a simple phenomeno-logical model based on ideas of chemical kinetics gives expressions for these functions which yield excellent fits to the M D results. This is very useful and provides a relatively simple physical picture of the selective solvation process. Essentially, the process can be separated into two parts, a rapid electrostriction phase during which the number of solvent particles (of both species) in the first shell increases, and a much slower redis-tribution phase during which the composition of the first shell comes to equilibrium. If one focuses on the first solvation shell, which is clearly the most important region, then the entire selective solvation process can be well described by two relaxation times: r e i , associated with the electrostriction, and r r e , associated with the redistribution of solvent species. Simulations and analyses aimed at gaining a detailed understanding of the microscopic mechanism by which species redistribution occurs in the first solvation shell have also been carried out. These calculations indicate that two processes are of near equal importance in the systems considered. The overall rate of redistribution depends on the rate of ion-solvent diffusion in the bulk and also on the rate of exchange at the first shell. Further, it was found that the exchanges occur on nearly a one-to-one basis, with only a very brief period between the insertion of a particle (Solvent S) and the ejection of another (Solvent W). In summary, a pleasingly simple physical picture of the dynamics of selective ion solvation emerges from the present calculations. It is by no means certain to what extent this picture will hold for more realistic models. For example, recent work [8,40] has shown that for ion solvation in pure methanol the formation of solvent-solute hydrogen Chapter 3. Solvation Dynamics in Simple Solvents 65 bonds involves relatively slow reorientational motion, and such effects may well appear as additional slow modes in mixed solvents where hydrogen bonding is important. The next chapter is aimed at answering this, as well as other, as yet unanswered questions. Chapter 4 Solvation Dynamics in Realistic Solvents The results of the previous chapter establish a foundation upon which a more de-tailed understanding of solvation dynamics can be built. As mentioned previously, the motivation for employing relatively simple models is to provide insight into the funda-mental underlying principles of the phenomenon. Without the unnecessary complication inherent in more sophisticated solvent models, the dominant factors which are important in solvation in these systems became apparent. Now, with the basic groundwork having been established, more complicated, and thus more realistic, models may be considered. While the results for the Stockmayer systems provided considerable insight into the nature of the solvent relaxation as well as the effect of polarity, it still left many questions unanswered. Relatively simple effects such as differences in molecular size were not investigated. Such differences might have a considerable effect on equilibrium preferential solvation. If one solvent species is considerably smaller than the other, it should be able to solvate an ion in greater number, thus increasing the observed affinity for that species. Other effects such as molecular geometry may also be important. Highly directional potentials, such as those involved in hydrogen-bonding species, have been shown by Landyi and Skaf [41] to have a considerable effect on solvent relaxation in other systems. Also of interest is asymmetry with respect to positive and negative charges. Many species are known to have vastly different solvating properties with respect to positively and negatively charged solutes. Thus, in a system which is otherwise identical, a simple change in charge sign may yield different results. 66 Chapter 4. Solvation Dynamics in Realistic Solvents 67 It was with these questions in mind that the systems considered in this chapter were chosen. The first set of such systems consists of a water-methanol solvent, and is useful as an initial investigation into realistic systems. These are relatively simple species, but are subject to very strong hydrogen-bonding interactions. While these interactions are known to be important on shorter timescale relaxation, it is not clear what effect they will have on the somewhat longer timescales of interest here. The second set of systems presented employs water and dimethylsulfoxide (DMSO) as solvents, and serve as a means of investigating the effect of molecular size and geometry, as well as charge asymmetry. 4.1 Model and Computational Details The systems considered consist of a single solute particle, s, with position r s , immersed in a binary solvent. The species involved are modelled using the standard site-site ap-proach discussed in Section 2.2.2. Two sites, i and j, with positions r; and r,-, interact through a potential function, Uij, consisting of both short-range (SR) and long-range (LR) components. The total configurational energy of the system is given by U = \ £ uij(rij) + YI uis(ris) , (4-1.1) ij i where the solute-solvent energy, uis, has been included as a separate term, = r-j — r,-, and ris = Tj — r s . Terms in the first summation corresponding to sites within the same molecule are omitted. The short-range potentials are of the usual Lennard-Jones (LJ) form. Each site may also have a point charge, (expressed in units of the elementary charge, e), resulting in long-range coulombic interactions, which are once again accounted for using the Ewald summation method (see Appendix B). As was the case in the pre-ceding chapter, we consider a single solute particle which is ionized only in the primary simulation cell. Thus, the long-range solute-solvent interaction, M^r, is the energy of Chapter 4. Solvation Dynamics in Realistic Solvents 68 insertion of a single point charge into an irregular lattice of point charges given by [37] A value of a — 5.4 was used in the current work, with all lattice vectors, k, satisfying k2 < 27 included in the summation. Several models have been developed for the species we wish to consider here. It should be noted that most model potentials used in condensed phases simulation studies are not designed to accurately reproduce the isolated, two-particle potential. Instead, they are developed so as to reproduce the known observed bulk properties of the liquid. For instance, it is experimentally known that the dipole moment of water is 1.85 Debyes in vacuo. However, in the liquid phase, the dipole moment increases to approximately 2.6 Debyes due to polarization effects from neighbouring molecules [42]. Thus, a non-polarizable model which is meant to be used for liquid phase studies would possess a higher dipole moment than would be required to reproduce the simple water-water pair potential. The other parameters in the model are typically adjusted so as to reproduce other quantities, such as density, self-diffusion constants, etc. A n additional requirement for any model which will be used in the current study is that they be transferable. That is, despite the fact that these models are typically developed with the pure liquid in mind, the potential must yield acceptable results when combined with other species. Much work has been done on transferable potentials for many species, including those of interest here. The widely used S P C / E model [43] will be adopted for water, along with the HI model of Haughney et al. [44] for methanol. The choice of models available for DMSO is somewhat smaller, and that of Vaisman and Berkowitz [45], which is based on neutron diffraction data as well as on ab initio calculations will be used here. The models for both methanol and DMSO employ the unified atom approach with respect to their methyl Chapter 4. Solvation Dynamics in Realistic Solvents 69 Table 4.1: Model parameters for the species considered. Distances are in Angstroms, energies in kJ/mole, charges in e, and angles in degrees. DMSO Methanol Water o 2.940 o 3.083 a 3.166 0 e 0.276 O t 0.731 O e 0.650 q -0.459 q -0.728 q -0.848 a 3.560 o - o -s e 0.844 H e - H e -Q 0.139 q 0.431 q 0.424 a 3.600 o 3.861 C H 3 e q 0.669 0.160 C H 3 e q 0.758 0.297 IOSC 107.2 LCOE 108.5 IHOH 104.9 LCSC 99.2 r(CO) 1.425 r(OH) 1.000 r(SO) 1.496 r{OH) 0.945 r(SC) 1.800 Solute o 4.417 e 0.493 q ±1.0 groups, thus removing any necessity to account for any bond rotations which would otherwise be required. Our solute model is simply that of the chloride ion used by Chandrasekhar et al. [46]. For comparison purposes, the same L J parameters were also employed for the positive and neutral solutes. While a positively charge chloride ion is obviously unrealistic, this model will allow for the investigation of the effects of charge sign without introducing other contributions which would arise from using a solute of different size. A summary of the parameters for each model is given in Table 4.1, and a schematic of each is shown in Fig. 4.1. A l l cross-species L J interactions were calculated using the Lorentz-Berthelot rules [30], aX2 = (an + cr22)/2 and e i 2 = v /ene 22-The simulation cell consisted of 1 solute and 255 solvent molecules, selected so as to correspond to the desired composition. In all cases the temperature was 298 K and the cell dimensions were chosen to be consistent with the experimental density at this temperature. The time integration was carried out using a standard Verlet leap-frog Chapter 4. Solvation Dynamics in Realistic Solvents 70 algorithm, with a timestep of 4 fs for equilibration purposes, and 1 fs for the trajectories themselves. Quaternions were employed to constrain molecular geometry. A l l systems were equilibrated at constant temperature for 1 ns, with equilibrium averages accumu-lated over a further 160 ps. For the neutral systems, a total of 10 configurations were extracted during this period. These were employed as initial points for ten nonequilib-rium trajectories, which were carried out at constant energy. As in the preceding chapter, we will consider both the solvation response function, S(t), as well as the particle sol-vation response function, P(t). Also of use are the equations obtained from the simple phenomenological model, both in their full Pj(t) = Cje-t/Tl + (1 - Cj)e-^n , (4.1.3a) PT(t) = CTe~t/Tl + (1 - CT)e-t/T2 , (4.1.36) and somewhat simplified Pj(t) = Cje-t/Tel + (1 - Cj)e-t/T" , (4.1.4a) PT(t) = e~t/T* , (4.1.46) forms. The relaxation times, re\ and r r e , can be obtained by fitting Eqs. (4.1.3) or (4.1.4), depending on the validity of the steady-state assumption for the system being considered. 4.2 Results and Discussion A summary of the systems considered is given in Table 4.2. The system names are meaningful. The first two letters specify which two solvent species are present (W -water, M - methanol, D - DMSO), while the number indicates the percent composition of the first solvent species (note that this is water for all systems discussed here). The + or — indicates the sign of the charge on the solute particle s (\qs \ = l.Oe), and N denotes Chapter 4. Solvation Dynamics in Realistic Solvents 71 Figure 4.1: A schematic of the 4 models employed in this chapter. The relative site sizes are drawn to scale, while molecular geometry is approximated. The hydrogen sites possess no Lennard-Jones component, and are represented by dashed circles. Chapter 4. Solvation Dynamics in Realistic Solvents 72 Table 4.2: A summary of information for the systems considered. The system names are as described in the text. rc is the cut-off distance and A / w a t and i V m e t h / D M S o a r e the average number of particles of a particular type found within rc of the solute at equilibrium. The relaxation times, re\ and r r e were determined by fitting the phenomenological equations. System r c (A) iVwat -^meth/DMSO Tel (PS) Tre (ps) WM50N 4.2 0.7 0.5 - -W M 5 0 - 4.2 5.1 2.0 3.3 27.8 WM50+ 4.2 6.2 2.2 0.5 24.1 WM10N 4.2 0.0 0.6 - -W M 1 0 - 4.2 1.3 4.5 6.6 104.5 WD50N 6.2 5.2 6.7 - -W D 5 0 - 6.2 9.3 6.6 0.4 130.0 WD50N 5.3 0.8 0.2 - -WD50+ 5.3 6.2 3.7 1.1 60.9 the neutral solute. Thus, for example, System WM10— consists of 10% water and 90% methanol with a negatively charged solute. 4.2.1 Water-Methanol The first five systems listed in Table 4.2 employ water-methanol solvents. These systems are of particular interest due to their potential for strong hydrogen bonding. The breaking and reforming of hydrogen bonds, which must occur to facilitate the restructuring of the solvation shells of the solute [41], might have a significant effect on the overall relaxation mechanism. Additionally, while S P C / E water is only marginally more polar than HI methanol (i.e., a dipole moment of 2.35 Debyes (D) versus 2.33 D), it is the smaller of the two solvents. The C H 3 group of methanol is quite bulky, which not only prevents strong electrostatic interaction of the methyl group with other species, but also makes the methanol molecule considerably more bulky as a whole. This, combined with methanol's slightly lower dipole moment, should result in water being the preferred solvent for both positively and negatively charged solutes. Chapter 4. Solvation Dynamics in Realistic Solvents 73 The equilibrium solute-solvent radial distribution functions for the first three related systems, consisting of 50% water and 50% methanol, are given in Fig. 4.2. As ex-pected, the neutral solute molecule exhibits a moderate preference for methanol; the more strongly interacting water molecules prefer the bulk solution where they can be completely surrounded by other polar species. Conversely, the positively and negatively charged solutes both exhibit a pronounced preference for water. It should be noted that the solvation shells about the neutral solute are considerably more relaxed than for the ionized cases. This contrasts with the earlier observations for Stockmayer solvents, where the difference in peak positions for the neutral and ionized solutes was smaller, resulting in less electrostriction following ionization. The somewhat different nature of the more realistic solvent models is not surprising given that the structure of bulk water itself is relatively loose compared to simple Stockmayer liquids. The significant discrepancy in the effective radius of the solvation shells of the neutral and charged solutes can also be seen in Table 4.2. Included in Table 4.2 are the compositions of the first solvation shell of the ionized solutes. Note that the radii r c defining the first shell are the same for positive and negative solutes in water-methanol mixtures. It can be seen from Table 4.2 that for the water-methanol systems there are on average few solvent particles of either species within rc of the neutral solute. Thus, it is clear that electrostriction is very significant for these systems. The importance of the electrostriction process can be seen clearly by examining the composition of the first solvation shell as a function of time. (Nj(t)) results for System WM50— are shown in Fig. 4.3 along with the solvation response functions, Sj(t), and the particle solvation response functions, Pj(t). As is evident from Table 4.2, the solvent density within the cut-off distance of 4.2 A at the moment of ionization is relatively small. However, immediately following ionization, this density increases rapidly as the solvent constricts under the newly created electric field. This is followed by the somewhat slower Chapter 4. Solvation Dynamics in Realistic Solvents 74 (r) 8.0 6.0 4.0 2.0 0.0 6.0 4.0 2.0 0.0 6.0 4.0 -2.0 -0.0 WM50N W a t e r • • • - M e t h a n o l WM50+ J I I I I IZ_J I I L J I L 0.0 2.0 4.0 6.0 8.0 10.0 r (A) Figure 4.2: Solute-solvent radial distribution functions for the systems consisting of 50% water and 50% methanol. Chapter 4. Solvation Dynamics in Realistic Solvents 75 1.0 Q Q I I I I I I I I I I I ' 0 10 20 30 40 50 t (ps) Figure 4.3: The solvation response functions, Sj(t), particle solvation response functions, Pj(t), and first shell coordination numbers, (Nj(t)), for System W M 5 0 - . The labels correspond to water (W), methanol (M), and the total (T). Chapter 4. Solvation Dynamics in Realistic Solvents 76 redistribution phase during which the total coordination number of the solute remains effectively constant, but the composition of the first solvation shell changes as methanol is partially replaced by water molecules. This behaviour is also apparent in the Pj(t) functions included in Fig. 4.3. The relaxation times obtained by fitting Eqs. (4.1.4) are given in Table 4.2. We note that, apart from the increased importance of electrostriction, the results for this system are in qualitative accord with the earlier observations for simple Stockmayer solvents. Thus, while hydrogen bonding might well influence fast reorientational relaxation, it does not appear to be of qualitative importance for the slower processes that govern selective solvation dynamics. Also of interest are the solvation response functions, Sj(i), included in Fig. 4.3. Sr(t) depends on the total ion-solvent energy and approaches equilibrium on the rapid elec-trostriction timescale with only a small contribution attributable to the slower redis-tribution phase. Although the two individual solvent species are energetically far from equilibrium [see Sw(t) and SM(t)], the contribution to the total ion-solvent energy from the redistribution process is quite small compared to the contributions from electrostric-tion and fast solvent reorientation. Nonetheless, the redistribution contribution is clearly not negligible since without it no preferential solvation would occur. These observations are again of particular importance for experimental studies which rely exclusively on the total solvation response function, ST(t), as a means of monitoring the relaxation. A l -though the Sr CO curve suggests that the relaxation is essentially complete after ~ 10 ps, it is obvious from the solvent-specific curves that a significant amount of solvent redistribution is still occurring. Thus, it is necessary to measure different properties if redistribution is to be studied experimentally, such as the spectral width relaxation discussed earlier [39]. The response functions for the related system, WM50+, are shown in Fig. 4.4. These Chapter 4. Solvation Dynamics in Realistic Solvents 77 results are qualitatively similar to those for the negatively charged ion, with one signif-icant difference. The rate of electrostriction is substantially faster for the positive case, r e i = 0.5 ps versus 3.3 ps, despite the fact that both systems share the same initial state and solvent characteristics. The reason for this is likely related to the affinity of the posi-tively charged solute for the oxygen site of water as opposed to the hydrogen sites. In the present model, the magnitude of the charge on the oxygen site is twice that of a hydro-gen site, leading to a considerably stronger ion-water interaction for positively charged solutes. The faster electrostriction is in all likelihood due to this stronger interaction. The two remaining water-methanol systems, consisting of 10% water and 90% metha-nol, were selected to investigate dilution effects. The initial and final solute-solvent radial distribution functions are shown in Fig. 4.5. As for the previous systems, the neutral solute shows a moderate preference for methanol, while the negatively charged solute is strongly solvated by water. As before, the effective diameter of the first solvation shell is markedly different for the neutral and charged cases, suggesting that substantial electrostriction occurs. The P(t) and S(t) functions for WM10— are shown in Fig. 4.6. The results are very similar to the 50/50 case, with rapid electrostriction accounting for most of the change in both the total coordination number and the total solute-solvent energy. The subsequent redistribution is substantially slower than at the higher water concentrations (rTe = 104.5 ps versus 27.8 ps) due to the scarcity of water in the immediate environment of the solute. This results in the rate of solvation being more strongly dependent on slower solvent diffusion than on simple exchanges between the first and second solvation shells. This is similar to what was found for dilute Stockmayer systems. It is apparent from these results that in water-methanol systems hydrogen bonding does not result in behaviour which differs greatly from that of the simpler Stockmayer solvents. Any hydrogen bond breaking and reforming most likely occurs on timescales Chapter 4. Solvation Dynamics in Realistic Solvents 78 Figure 4.4: The solvation response functions, Sj(t), and particle solvation response func-tions, Pj(t), for System WM50+. The labels are as in Fig. 4.3. Chapter 4. Solvation Dynamics in Realistic Solvents 79 15.0 10.0 -5.0 g(r) 0.0 10.0 h 5.0 0.0 -WM10N — W a t e r 1 1 1 1 1 1 • • • • M e t h a n o l 1 WM10-1 1 1 1 1 1 ; i T I I I M ^ I i i i I T ^ - I 0.0 2.0 4.0 6.0 8.0 10.0 r (A) Figure 4.5: Solute-solvent radial distribution functions for the systems consisting of 10% water and 90% methanol. Chapter 4. Solvation Dynamics in Realistic Solvents 80 J I I I I I I I I I I I I I I I I I L 0 20 40 60 80 100 t (ps) Figure 4.6: The solvation response functions, Sj(i), and particle solvation response func-tions, Pj(t), for System W M 1 0 - . The labels are as in Fig. 4.3. Chapter 4. Solvation Dynamics in Realistic Solvents 81 which are considerably faster than the overall solvent relaxation, at least for the models considered here. 4.2.2 Water-Dimethylsulfoxide From the discussion above it is clear that, despite the presence of hydrogen bonding, ion solvation dynamics in water-methanol systems strongly resembles that observed in bi-nary Stockmayer solvents. To explore further, water-DMSO mixtures will be considered. Water and DMSO are geometrically less alike than water and methanol and, hence, one might expect more significant deviation from the simple dipolar solvent picture. Unlike HI methanol, the dipole moment of the DMSO model is considerably greater than that of S P C / E water, 4.37 D versus 2.35 D. While this might serve to make DMSO the solvent preferred by charged solutes, its size and bulky geometry might significantly inhibit its ability to solvate relatively small solute particles. Thus, from simple consideration of the models it is not immediately obvious which of the two solvent species will be preferred by either the neutral or charged solute. The solute-solvent radial distribution functions for these systems are shown in Fig. 4.7. It is evident that water is the solvent preferred by the ion, most likely due to its small size and hence its ability to solvate in greater numbers than DMSO. More interesting, however, are some features apparent in these systems that distinguish them from the water-methanol case. While the neutral solute shows the typical moderate preference for one solvent (DMSO), and the charged solutes both show a considerable preference for the other solvent (water), the nature of the DMSO structure in System WD50— is somewhat unexpected. In particular, the first peak in the solute-DMSO radial distribution function remains virtually unchanged from the neutral case. Charging the solute results in only a slight sharpening of the first peak. In contrast, the solute-DMSO distribution function for the positively charged solute in System WD50+ exhibits a peak which has not only Chapter 4. Solvation Dynamics in Realistic Solvents 82 become more well defined, but has also shifted position considerably. This is more akin to previous results, and points once again to an important electrostriction contribution for both species. If the DMSO structure about the solute for System WD50— is examined more closely, it becomes apparent that little or no reorientation occurs when the solute is charged. The solute-DMSO distribution functions for the individual interaction sites of DMSO are shown in Fig. 4.8. It is immediately obvious that the DMSO orientation is virtually the same for the neutral and negatively charged solutes, with the large methyl groups directed inward in both instances. The preference for the methyl groups by the negative ion is not surprising as these sites are the only easily accessible positions on the D M S O model having positive partial charges. The preference of the neutral solute for the methyl groups is likely due to the relatively small magnitude of the site charges (\q\ = 0.16e). The fact that DMSO has no site with a strong positive charge implies that it is not a particularly good solvent for negative ions. Its ability to solvate negative ions is further hindered by the large size of the methyl groups (o = 3.6 A), and the combined effects result in the effective exclusion of DMSO from the first solvation shell as is evident in Figs. 4.7 and 4.8. In contrast, the solute-DMSO site-site distribution functions for System WD50+ differ considerably from the neutral solute case. DMSO solvates positive ions with the oxygen directed inward at a distance of ~ 3.2 A, roughly the same distance as the solute-water peak in Fig. 4.7. The small size (a = 2.94 A) and strong negative charge (q = —0.459e) of the DMSO oxygen serve to make DMSO an effective solvent for positive ions. In addition, unlike the negative ion case, the large methyl groups are now directed outwards and thus do not hinder the ability of DMSO to penetrate into the first solvation shell of the solute. The fact that DMSO solvates positive ions much more strongly than negative ions was deduced long ago from the observation that negative ions are more highly reactive in solution [47]. The higher reactivity is attributed to the fact Chapter 4. Solvation Dynamics in Realistic Solvents 83 12.0 9.0 6.0 3.0 0.0 9.0 g(r) 6.0 3.0 0.0 9.0 6.0 3.0 0.0 WD50N Water DMSO J I 1 I I L ^ 2 ^ „ _ J I I I I I I I I I I J I I I I Zl WD50-J I I I I I L WD50 + J I I I L J I L J I L 0.0 2.0 4.0 6.0 8.0 10.0 r (A) Figure 4.7: Solute-solvent radial distribution functions for the systems consisting of 50% water and 50% DMSO. Chapter 4. Solvation Dynamics in Realistic Solvents 84 that the negative ions are relatively unencumbered by solvent molecules. Equilibrium snapshots taken directly from the simulations showing the immediate environment of the solute particle are also instructive. Fig. 4.9 shows an instantaneous configuration of the solvent particles about the negative ion, while Fig. 4.10 shows one such configuration for the positive ion. It can be seen from Fig. 4.9 that, as was evident in the pair distribution functions [Fig. 4.7], the effective distance at which the DMSO molecules solvate the negative ion is quite large. In fact, it can be seen that the shell of water molecules does indeed appear to be contained within an outer layer of DMSO. This is in stark contrast to the situation shown in Fig. 4.10, where the DMSO solvate the positive ion much more tightly. In this case, water no longer fills the region between the DMSO and the ion, but instead resides in the regions where DMSO is absent. Thus, this bears more in common to the previous Stockmayer and water-methanol systems where the two species directly compete for solvation of the ion. The different behaviour of System W D 5 0 - is also reflected in the solvation dynamics following ionization. The raw first shell coordination numbers, (Nj(t)), and the solute-solvent energy components, (Uj(t)}, are instructive for this system and are shown in Fig. 4.11. From the (ND(t)) curve, one can immediately infer that there is little change in the DMSO structure during the solvation process. Meanwhile, the number of water molecules, (Nw(t)), in the first solvation shell slowly increases as more water is drawn from the surrounding bulk. Interestingly, the degree of fast electrostriction evident for both species appears to be negligible. The result is relaxation on a single timescale determined by the rate at which water is able to penetrate through the essentially static DMSO layer. Consequently, the total coordination number of the solute does not reach equilibrium on the rapid electrostriction timescale as in all systems previously considered, but rather is strongly dependent on the much slower redistribution phase. Chapter 4. Solvation Dynamics in Realistic Solvents 85 Figure 4.8: Solute-solvent site-site pair distribution functions for DMSO in Systems WD50N, W D 5 0 - and WD50+. Figure 4.9: A snapshot showing all (a) DMSO and (b) water molecules present within a distance of 7.0 A of the negatively charged ion. The colours correspond to methyl (black), sulfur (orange), oxygen (blue), hydrogen (red), as well as the ion (green). Figure 4.10: A snapshot showing all (a) DMSO and (b) water molecules present within a distance of 7.0 A of the positively charged ion. The colours are as in Fig. 4.9. Chapter 4. Solvation Dynamics in Realistic Solvents 88 0 100 -200 w I- \ A <Uj(t)> (kJ/mol) -300 -400 -500 10 h <Nj(t)> D l l l l I I I I I I I I l l I' I l l l I l I I I 0 I I I I I ' I I ' I I I I 0 10 20 30 40 50 t (ps) Figure 4.11: Solute-solvent energy components, (Uj(t)), and first shell coordination num-bers, (Nj(i)), for System W D 5 0 - . The labels correspond to water (W), DMSO (D) and the total (T). Chapter 4. Solvation Dynamics in Realistic Solvents 89 Alternatively, one might consider that the relaxation process is solely the slow elec-trostriction of water in a static medium of DMSO. However, this description is not com-pletely satisfactory when one examines the solute-solvent energy components, (Uj(t)) (Fig. 4.11). It is evident from these results that the DMSO is not entirely unaffected by the relaxation process. The small upward drift in (f/rj(i)) indicates that some changes in the DMSO structure are taking place. However, unlike the coordination numbers dis-cussed above, the energy is not dependent on the first shell alone, but rather includes contributions from the entire system. Therefore, the drift in (Ur)(t)) likely results from slight orientational shifts of the DMSO molecules in response to the increasing density of water in the region between them and the solute particle. The peak in the solute-DMSO-methyl distribution function (Fig. 4.8) occurs at a distance of approximately 3.8 A versus the water peak at 3.2 A . Thus, as more water penetrates to the solute, the orientation of a DMSO molecule is probably determined more by its ability to "solvate" the water oxygens than by its interaction with the actual solute particle. This is supported by the snapshot shown in Fig. 4.9, where it is apparent that the DMSO does indeed appear to be separated from the ion by a layer of water. Solvation and particle response functions for System WD50— are shown in Fig. 4.12. Note that, due to the fact that the number of DMSO molecules in the first solvation shell remains effectively unchanged, -Po(i) is not well-defined and has been omitted from the plot. Also, because (NT(t)) depends very strongly on the redistribution process, the simplified phenomenological Eqs. (4.1.4) are no longer adequate. These equations assume that only one timescale, specifically that of the rapid electrostriction phase, determines Pr(t)- This is obviously not the case for the current system and the more general Eqs. (4.1.3) were used for fitting purposes. It is apparent that neither the total coordination number nor the total solute-solvent energy has reached equilibrium in the 50 ps simulated. Again, this is in stark contrast to the situation for all other systems considered, both the Chapter 4. Solvation Dynamics in Realistic Solvents 90 water-methanol case discussed above as well as the Stockmayer models of Chapter 3. It should be pointed out that the cut-off employed for this system, rc = 6.2 A , was chosen so as to include the first peak in the solute-DMSO radial distribution function. This has the effect of including both the first and second peaks of the solute-water radial distribution function. Results obtained using a cut-off of 3.9 A , which includes only the first water peak, were qualitatively similar to those shown here, again exhibiting a constant DMSO density throughout the entire solvation process. The response functions for System WD50+ are shown in Fig. 4.13. The results for this system have much more in common with the water-methanol case than with System WD50—. The total coordination number of the solute reaches equilibrium on the rapid electrostriction timescale, with the subsequent redistribution occurring on an essentially one-to-one basis despite the much larger size of the DMSO molecules. The reason for this can be seen in the site-site distribution functions for the system shown in Fig. 4.8, as well as the equilibrium snapshot shown in Fig. 4.10. As noted above, DMSO solvates the positive ion with its relatively small oxygen site directed inward; the bulk of the DMSO molecule (the methyl groups in particular) are located well outside what one might consider to be the first solvation shell of the ion. Hence, the effective volume which a DMSO molecule occupies in the first shell is very similar to that occupied by a water molecule, leading to the observed one-to-one exchange of water and DMSO. This results in dynamical behaviour comparable with that found in solvents of equal sized Stockmayer particles and in water-methanol systems. Chapter 4. Solvation Dynamics in Realistic Solvents 91 Figure 4.12: The solvation response functions, Sj(t), and particle solvation response functions, Pj(t), for System WD50—. The labels are as in Fig. 4.11. Chapter 4. Solvation Dynamics in Realistic Solvents 92 _ 4 Q r i i i i I i i' i i 1 i i i i I i i i i I i i i i 0 10 20 30 40 50 t (ps) Figure 4.13: The solvation response functions, Sj(t), and particle solvation response functions, Pj(t), for System WD50+. The labels are as in Fig. 4.11. Chapter 4. Solvation Dynamics in Realistic Solvents 93 4.3 Summary In this chapter the results of molecular dynamics simulations of ion solvation in water-methanol and water-DMSO systems have been presented. It was found that the water-methanol binary solvent, despite the presence of strong hydrogen bonding, exhibits re-laxation behaviour qualitatively similar to the Stockmayer results presented in Chapter 3. Apart from the increased importance of electrostriction due to the somewhat looser packing inherent in realistic solvents, the relaxation again consists of rapid electrostric-tion and slow redistribution phases. As in the Stockmayer systems, the total coordination number of the solute remains essentially constant during the entire redistribution phase, implying that the solvent exchange occurs on a one-to-one basis. The water-DMSO systems were considered mainly to examine the influence of sol-vent size and shape on the relaxation dynamics. For these systems, it was immediately apparent that the relaxation behaviour depended very strongly on the sign of the solute charge. Due to its highly asymmetric nature, DMSO is much better able to solvate pos-itive than negative ions. Consequently, the negatively charged case exhibits dynamical solvation behaviour quite unlike that of the water-methanol and simple Stockmayer sys-tems. The DMSO responds only marginally when the solute acquires a negative charge, and the relaxation consists almost exclusively of water penetrating to the solute through a largely static DMSO medium. As a result, for negative solutes the total coordination number does not reach equilibrium rapidly, but rather depends on the slow redistribu-tion timescale. In contrast, for a positively charged solute the relaxation behaviour is similar to that observed in water-methanol and Stockmayer systems. DMSO solvates positive ions with the oxygen directed inward and the larger methyl groups positioned well outside the first solvation shell. Therefore, despite its size, DMSO is able to solvate a positive ion occupying roughly the same effective volume in the first solvation shell as Chapter 4. Solvation Dynamics in Realistic Solvents 94 the considerably smaller water molecules. The result is a one-to-one exchange process much like that observed for solvent species of similar size. Finally, of particular interest are the total solvation response functions, Sr{t), as these are often obtained in experimental studies of solvation dynamics. Somewhat unfor-tunately, the present results once again suggest that very little of the solvent relaxation is actually evident in this function. Much like in the case of the Stockmayer systems investigated in Chapter 3 , the energetic contribution from the preferential solvation is quite small compared to the contributions from reorientation and electrostriction, and hence ST{t) can be very insensitive to the composition of the first shell. This renders 5T (0 less useful for mixtures than for one component solvents. Chapter 5 A Comparison with Time Dependent Density Functional Theory As mentioned previously, computer simulations are but one tool available for the study of liquid phase phenomena. Many other techniques are also available. One such technique which has received considerable attention in recent years is that of time dependent density functional theory (TDDFT) . These approaches have proven very useful in the study of many diverse topics, including solvation in pure solvents. Such studies, and studies based on related theories, have yielded results which agree well with computer simulations and have provided further insights into solvation dynamics. In contrast to the work done on pure solvents, theoretical treatments of mixed solvents have been scarce. Chandra [11] has carried out a preliminary theoretical and computer simulation treatment of mixed Stockmayer solvents. However, this work focussed on the standard solvation response function, S(t), which has now been shown to be somewhat inadequate for characterizing the behaviour of these systems. Additionally, the solvent species considered by Chandra were only moderately different in their characteristics, leading to a relatively weak selective solvation effects. The work in this chapter is the result of collaboration with Dr. Akira Yoshimori of Kyushu University. The aim is to apply T D D F T calculations, as performed by Dr. Yoshimori, to the study of solvation dynamics in mixed solvents and compare the results with those of computer simulations on equivalent systems. This not only gives more insight into the nature of the solvent relaxation, but also provides a means of investigating the validity and limitations of the theory. 95 Chapter 5. A Comparison with Time Dependent Density Functional Theory 96 To allow for the analytical solution of the T D D F T equations with a minimum of ap-proximations, spherically symmetric particles will be of primary interest. These particles will interact via the standard Lennard-Jones potential described previously. For simplic-ity, in the initial equilibrium state the L J interactions between all species are taken to be the same. Then, at time t = 0, the interaction of the solute with one solvent species is in-creased by increasing the energy parameter, e, in the appropriate solute-solvent potential. A l l other interactions remain unchanged. The system then evolves in time reaching its final equilibrium state as t —> oo. This is a very simple model but it allows us to cleanly study the evolution of the solvent densities about the solute as the system evolves. A n extension to the Stockmayer systems presented in Chapter 3 will also be briefly discussed. 5.1 Theoretical Overview Consider a solute particle, X , immersed in a binary solvent. For simplicity, the solute is assumed to be stationary and to be located at the origin of the coordinate system. In the Markovian limit, the time evolution of the solvent densities, pQ(r,t) (a—1,2), is given by [48-55] where F is the interaction portion of the free energy functional, kB is the Boltzmann constant, and T is the temperature. A derivation of Eq. (5.1.1) is given in Appendix D. The above equation is formulated in terms of the parameter, Db. The original work due to Munakata [51,53] referred to this quantity as the "bare diffusion constant". While it is undoubtedly related to the more familiar self-diffusion constant, D, the nature of the relation between the two is not clear. Unfortunately, as it is this parameter which determines the absolute timescale of the theoretical equations, it will be seen that this presents a considerable problem when attempting to compare with simulation or (5.1.1) Chapter 5. A Comparison with Time Dependent Density Functional Theory 97 experimental results. By retaining only the term in the free energy functional which is dependent on pair correlations, thus neglecting higher order terms, one obtains [48-54] 8J-lh*L = l n p a + ^(r) -J2 [dr' cQ / 3(r - r')^(r ' , t ) , (5.1.2) where 5pa = pa — pa is the difference in density for species a at position r at time t from its bulk density, pa. Additionally, </>*(r) = (t>a(v)/kBT where </>Q(r) is the interaction potential of the solute with species a and cap(r) is the direct correlation function for the two solvent species. Combining Eq. (5.1.2) with Eq. (5.1.1) yields ^ = A(v 2 p« + V [ p aV^(r)] - V [ p a V f dr' caf3(r - r') 8Pp(T', t)]). (5.1.3) 0=1 J Due to the spherical symmetry of the present system, this can be simplified to yield 8A.(r,«) _ 1 0^, ( 5 _ 1 - 4 B ) dt r2 dr Ja(r,t) = Db dpa(r, t) + p ^ ^ /d<f>*a(r) + dipa(r,t) (5.1.4b) dr \ dr dr Mr,t) = - £ fdx' ca0(\r - v'\)8Pp{r'\t) . (5.1.4c) These nonlinear equations can then be solved numerically for a given system. Alterna-tively, Eqs. (5.1.4) may be linearized with respect to both density fluctuations as well as the solute-solvent potential, (f>(r), to yield [48-52,56] d 5 p a ^ t ] = Db[V28pa{r,t) + pa^a{r,t)+paV2^a{r)}. (5.1.5) This expression reproduces Debye-Huckel type behaviour [35] in the limit of t —> oo. This is known to be a somewhat poor approximation, however, and an alternative is desirable. This can be accomplished by replacing the bare interaction, (j>*(r), with the Chapter 5. A Comparison with Time Dependent Density Functional Theory 98 solute-solvent direct correlation function, ca(r) [57-61]. This function can be regarded as an "effective interaction". Whereas the use of the bare interaction might be sufficient for predicting the solute-solvent distribution in low density phases, it is a poor approx-imation in condensed phase system such as liquids. In such systems, the solute-solvent distribution is determined by not only the direct interaction, but also by the structure of the liquid itself. The solute-solvent direct correlation function provides a means of accounting for these effects. The use of the direct correlation function yields d S P a £ ' t ] = Db[V25pa(r,t) + paV2Mr,t) - p Q V 2 c Q ( r ) ] . (5.1.6) This linearized simplification of Eqs. (5.1.4) has the advantage of possessing an analytical solution. If one takes the spatial Fourier transform of Eq. (5.1.6) and notes that, for the present choice of models, all solvent-solvent direct correlation functions are equal, one obtains the solution Spa(k,t) = Axe-Xlt + A2e~X2t + 8pa(k, oo) , (5.1.7a) where 2 Ax = Za£[%(M)-<JMfc,oo)], (5.1.7b) 0=1 A2 = xp[5pa(k, 0) - Spa(k, oo)] - xa[8pp(k, 0) - 8p0(k, oo)] , for a ^ (5 ,(5.1.7c) Ai = k2Db[l - pc{k)\ , (5.1.7d) A 2 = k2Db . (5.1.7e) and 8pa(k, 0) and 8pa(k, oo) are the initial and final values of 8pa(k, t), xa = pa/(pi+p2). Additionally, 8pa(k, oo) = paha(k), where ka(k) is the Fourier transform of the solute-solvent pair correlation function at final equilibrium. Interestingly, the solution given by Eqs. (5.1.7) exhibits two timescales inherent in the relaxation, as determined by Ai and A 2 . From the form of the expression for Ai it is apparent that the first timescale Chapter 5. A Comparison with Time Dependent Density Functional Theory 99 corresponds to the change in the total density of the solvent about the solute. Similarly, the second timescale must then correspond to exchange of the two solvent species at constant density. This is precisely what was found in the simulation results of the previous chapters. However, it should be noted that the relaxation times given by Eqs. (5.1.7d) and (5.1.7e) are k dependent. Consequently, when translated back to r space, this clear separation of timescales is no longer evident. Nonetheless, it does indicate that, within the approximations made to obtain this result, the solvent relaxation is indeed dominated by these two effects. In addition to the solution to the linearized equations, numerical solutions to the non-linear form, given by Eq. (5.1.4), were also obtained. In both instances, the equilibrium structure and direct correlation functions were obtained using hypernetted chain (HNC) methods. 5.1.1 The Self-Diffusion Constant As mentioned above, the timescale of the solutions to both the linear and nonlinear theories is determined by the bare diffusion constant, Df,- The nature of the relationship between this quantity and the experimentally measurable self-diffusion constant is not obvious. Unfortunately, much work [48,49,61-69] which followed the original development of the theory has made the implicit assumption that D and Db were, in fact, equivalent within the framework of the linearized theory, without explicitly addressing the issue. The following section is meant to investigate the true nature of the relationship between these two quantities. For simplicity, a one-component system is considered. Of use in the following derivation is the van Hove function, G(r,t), which, in an N particle system, is defined as [35] G(T,t) = (p(r,t)p(O,0)> (5.1.8) P Chapter 5. A Comparison with Time Dependent Density Functional Theory 100 1 / N N \ = ^ ( E E ^ + ^ O l - ^ i ) . (5-1-9) where p(r, t) is the density at position r at time t, (p(v, t)p(0, 0)) is the density correlation function, and 8(x) is the standard delta function, which has the properties 5(x) =0 if x^O (5.1.10) and /oo 6{x)dx = l. (5.1.11) -oo The van Hove function is effectively a much more generalized version of the radial dis-tribution function, g(r), presented earlier. It is a measure of the probability of finding particle i at position r after time t given that particle j was at the origin at time t = 0. This formulation does not exclude the possibility that i = j. Consequently, Eq. (5.1.9) can be separated into two contributions [35] °s(r,t) = jj(jt8[T + Ti(0)-Ti(t)]} (5.1.12a) G*(r,t) = ^(J2ES[r + rJ(0)-^)})- (5-1-126) where Gs(r,t) and Gd(r,t) are the "self" and "distinct" parts of the van Hove function. The self part accounts for those cases where i — j, and is thus the probability of finding particle i at position r after time t assuming it was at the origin at time t = 0. The distinct part accounts for the remaining possibility where i ^ j (and thus, the particles are "distinct"). In addition to the standard van Hove function, also of use are its spatial Fourier transform G ( M ) •= JG(r,t)e-ikrdr = i(PkP-k), (5.1.13) Chapter 5. A Comparison with Time Dependent Density Functional Theory 101 and the subsequent Laplace transform 1 r°° G(k,co) = — G(k,t)elU}tdt. (5.1.14) To obtain an expression for the self-diffusion constant, D, we consider a tagged particle with position r(0) at time t — 0. If, at time t, the particle has diffused to position r(t), then the self-diffusion constant is given by the familiar Einstein relation [35] D= l i m - J - < | r ( i ) - r ( 0 ) | 2 > , (5.1.15) where < \r(t) — r(0) | 2 > is the mean-square displacement of a reference particle at time t. The position at time t, r(t), is related to the particle velocity, v(i) through the simple relation r(t) =r(0) + j\{t')dt\ (5.1.16) Jo which, when squared and averaged over initial points yields (|r(t) - r(0)| 2) = / ' f\v(t") • v(t'))dt'dt". (5.1.17) ./o Jo Substitution into Eq. (5.1.15) yields D = lim - /* f\v(t") • v(t'))dt'dt". (5.1.18) i-+oo Qt Jo Jo A simple change of variables, r = t' — t", followed by integration by parts yields the familiar Green-Kubo formula for the self-diffusion constant D = lim — J->oo Qt 2 / ( t - r ) ( v ( r ) • v(0))d? . Jo 1 r°° = 3Jo < v ( r ) ' V ( 0 ) ) d T - ( 5 - L 1 9 ) Now, consider a set of n tagged particles which are a subset of the entire N particle system. These tagged particles are assumed to be of sufficiently low density such that Chapter 5. A Comparison with Time Dependent Density Functional Theory 102 they don't interact with each other. The movement of these particles over time can be described macroscopically by Fick's Law [35] j(r,t) = -Z)Vp (r , i ) , (5.1.20) where j(r, t) is the velocity density, or current, given by j(r,i) = f>(i)<5(r-r , ( i ) ) . (5.1.21) i Additionally, the system must obey the continuity condition * M = - V J M ) (5.1.22) which simply ensures conservation of particle density. Eqs. (5.1.20) and (5.1.22), when combined, yield the well-known diffusion equation ^ ^ = D V 2 p M ) , (5.1.23) or its reciprocal-space equivalent ^ = -Dk2Pk(t) (5.1.24) which has the solution Pk(t) = P*e-Dk2\ (5.1.25) where pk = Pk(0). Multplying Eq. (5.1.25) by p_fc and averaging over n particles yields [35] n n - e-DkH, (5.1.26) where the fact that the positions of the tagged particles are uncorrelated has been used. Chapter 5. A Comparison with Time Dependent Density Functional Theory 1 0 3 While it may be tempting, based on Eq. ( 5 . 1 . 1 3 ) , to equate Eq. ( 5 . 1 . 2 6 ) directly with the van Hove function, certain considerations must be made. The derivation of Eq. ( 5 . 1 . 2 6 ) is based on the macroscopic diffusion equation, while the van Hove function is a microscopic quantity. Consequently, a direct relation between the two is only possible in the hydrodynamic (long wavelength and low frequency) limit [35]. Thus, lim lim G,(k, t) = e~Dk2\ ( 5 . 1 . 2 7 ) w-H) fc-»0 where the self part of the van Hove function, instead of the van Hove function itself, has been used. Thus is due to the assumption that the tagged particles are at sufficiently low density so as to be completely uncorrelated. Eq. ( 5 . 1 . 2 7 ) can be rewritten it terms of its Laplace transform, Gs(k,co), as [35] 1 Dk2 lim l i m G c f k , c o ) = / T ^ , 0 , 0 . ( 5 . 1 . 2 8 ) w - > o fc^o n ' ' ircu2 + {Dk2)2 y 1 Thus, it can be immediately seen that the self diffusion constant, D, can be related to the self part of the van Hove function by [35] to2 D = lim lim — 7 r G s ( k , to). ( 5 . 1 . 2 9 ) w->ofc->o k2 A similar relationship between the bare diffusion constant, Db, and the van Hove function must also be obtained. For the current system, it has been shown that the Fourier-Laplace transform of the full van Hove function is given by [66] G{k,u) = — Vfl^/l V I , 2 > ( 5 - L 3 ° ) v ' -iw + Db[l - pc(k)]k2 where h(k) is the Fourier transform of the pair correlation function. Also of use is the Kerr approximation, which relates the self part of the van Hove function to the full, given by [35] v ' 1 - pc(k)[iujGsk,uo) +1] v ' Chapter 5. A Comparison with Time Dependent Density Functional Theory 104 Thus, combining Eqs. (5.1.30) and (5.1.31) yields G * < K - > = J T W ( 5 1 - 3 2 ) which, upon substituting into Eq. (5.1.29) and taking the limits yields the simple result D = Db. (5.1.33) Thus, it turns out that, for the linear theory, this equivalence is only valid subject to the Kerr approximation, and is not a result of the linearization itself. An alternative approximation for the relationship between D and Db has also be obtained by Araki and Munakata [54]. The details of the derivation will not be presented here. Only the result D_ _ 3(2TT) 3 Db pjh(k)2dk' { ' is of interest. 5.2 Results and Discussion As discussed above, the system of interest consists of a single solute particle, X, immersed in a binary solvent. The mass, m, is taken to be the same for all species. In the initial state all interactions are identical and hence the solute-solvent distribution is the same for both species. At time t = 0, this symmetry is broken by increasing the interaction between one solvent species and the solute while leaving all other interactions unchanged. This solvent species will be labelled S because it has the stronger interaction with the solute. More explicitly, at t — 0, cXs 1 S changed to e'xs > eXs a n d all other parameters remain the same. The remaining solvent species which interacts more weakly with the solute is labeled W. The evolution to equilibrium after breaking the solvent symmetry may be monitored in various ways, but here we will again focus on the size and composition of the first solvation Chapter 5. A Comparison with Time Dependent Density Functional Theory 105 shell of the solute particle as a function of time. The first solvation shell is defined to include all particles which lie within a specific distance r c of the solute. In the present calculations a value of r c = 1.3a will be used, (note that for the current model oXs = o x w = aSs = oww = cr) which corresponds approximately to the first minimum in the solute-solvent pair distribution function. Once again, the particle solvation response function for species j, Pj(t), can be defined as ^ ( t ) - ( ^ . ( 0 ) - i V , ( o o ) ) ' ( 5 - 2 ' 1 } where (Nj(t)) is the number of Solvent j particles in the first solvation shell at time t. We will set j — W,S, or T to refer to Solvent W, Solvent S, or the total coordination number, respectively. A l l systems considered were characterized by the reduced density p* = po3 = 0.8 and the reduced temperature of T* = kBT/e = 1.35, where e = eXs — ^xw = £ss = cww-The M D simulations were carried out with one solute and 255 solvent particles appropri-ately labeled to reflect the desired composition. The reduced timestep t* = (t/mo2)ll2t = 0.0025 was employed. The initial state for each system was equilibrated for a total of 200000 timesteps at constant temperature. A further 800000 timesteps were calculated during which 160 configurations were extracted at intervals of 5000 timesteps. These configurations were used as initial conditions for the actual nonequilibrium simulations involving the modified solute-solvent interactions. The nonequilibrium calculations were carried out at constant energy. A l l results presented here represent an average over 160 trajectories. The results reported were obtained allowing the solute particle to move con-sistent with physical solvation. Calculations holding the solute particle fixed as assumed in the T D D F T were also carried out and the Pj (t) curves obtained decayed a little more slowly than for the moving solute case. However, the effect of fixing the solute is very small (only marginally distinguishable from statistical error) and is not an important Chapter 5. A Comparison with Time Dependent Density Functional Theory 106 Table 5.1: Mole fractions, xs,w, and interaction parameters for the systems considered. System xs xw e'xs/exs 1 0.25 0.75 4 2 0.05 0.95 4 3 0.25 0.75 6 4 0.05 0.95 6 contribution to the discrepancies between theory and simulation discussed below. Relevant parameters for the four systems studied are summarized in Table 5.1. In all cases the composition of the first shell changed significantly as the system evolved from the initial to the final state. As an example, the initial and final solute-solvent radial distribution functions are shown in Fig. 5.1. We note a rather dramatic change in the first shell composition which is almost entirely Solvent S in the final state. The particle response functions obtained from both the linear and nonlinear versions of T D D F T are compared with simulations in Figs. 5.2-5.5 and 5.6-5.9, respectively. The reduced time scale tD/o2, where D = 0.017(ea 2/m) 1 / ' 2 is the M D self-diffusion constant, is used in the plots. Since the theory is formulated in terms of the "bare diffusion constant", Db, the ratio D/Db must be known in order to determine the absolute timescale for the theoretical curves. For each system considered, three sets of theoretical curves are given representing three different values of the D/Db ratio. The first set of curves assume D/Db = 1, while the second set correspond to the value of D/Db obtained from T D D F T using the approximation of Araki and Munakata. The final set of curves was obtained using a value of D/Db chosen to yield the optimal fit to the simulation results. This value was determined by minimizing the area between corresponding pairs of theoretical and simulation curves. For the purposes of this fitting, all three curves were considered simultaneously except where Pr(t) was too noisy to be useful. Several observations can be made from Figs. 5.2-5.5. Overall, the linear theory does Figure 5.1: The equilibrium pair distribution functions for (a) the initial state and (b) the final state, as obtained from simulation with xs = 0.25 and c'xs/exs = 6. The solid curves are for Solvent W and the dashed curves for Solvent S. Chapter 5. A Comparison with Time Dependent Density Functional Theory 108 0.00 0.05 0.10 0.15 0.20 0.25 tD/cr 2 Figure 5.2: A comparison of particle solvation response functions obtained from simula-tion and linearized theory for System 1. In the upper, middle and lower panels the D/Db ratio is 1, the value given by the theory of Araki and Munakata, and the value obtained by fitting, respectively. In each panel, from top-to-bottom, the curves are Pw{t), Ps(t) and PT(t). Chapter 5. A Comparison with Time Dependent Density Functional Theory 109 0.00 0.05 0.10 0.15 0.20 0.25 tD/o-2 Figure 5.3: A comparison of particle solvation response functions obtained from simula-tion and linearized theory for System 2. Curves and panels are as in Fig. 5.2, although the Pr(t) curve has been omitted, as the statistical noise prevented any useful consideration. Chapter 5. A Comparison with Time Dependent Density Functional Theory 110 Figure 5.4: A comparison of particle solvation response functions obtained from simula-tion and linearized theory for System 3. Curves and panels are as in Fig. 5.2. Chapter 5. A Comparison with Time Dependent Density Functional Theory 111 1.0 0.0 1.0 PL, 0.0 1.0 0.5 0.0 1 1 1 1 1 1 1 D/D b = 1 i i 1 i i i i 1 i i i i 1 i i i i 1 1 1 1 1 1 1 D/D b = 0.421 i i l i i i i l i i i i l i i i i 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0.00 0.05 0.10 0.15 0.20 0.25 ' 2 o tD/< Figure 5.5: A comparison of particle solvation response functions obtained from simula-tion and linearized theory for System 4. Curves and panels are as in Fig. 5.2, although the Pr(t) curve has been omitted, as the statistical noise prevented any useful consideration. Chapter 5. A Comparison with Time Dependent Density Functional Theory 112 0.00 0.05 0.10 0.15 0.20 0.25 tD/cr 2 Figure 5.6: A comparison of particle solvation response functions obtained from simula-tion and nonlinear theory for System 1. Curves and panels are as in Fig. 5.2. Chapter 5. A Comparison with Time Dependent Density Functional Theory 113 0.00 0.05 0.10 0.15 0.20 0.25 tD/a 2 Figure 5.7: A comparison of particle solvation response functions obtained from simula-tion and nonlinear theory for System 2. Curves and panels are as in Fig. 5.2, although the Pr(t) curve has been omitted, as the statistical noise prevented any useful consideration. Chapter 5. A Comparison with Time Dependent Density Functional Theory 114 Q o L- 1 1 1 1 ' 1 1 1 1 ' 1 1 1 1 I 1 1 1 ' 1 1 1 1 I 0.00 0.05 0.10 0.15 0.20 0.25 tD/CT2 Figure 5.8: A comparison of particle solvation response functions obtained from simula-tion and nonlinear theory for System 3. Curves and panels are as in Fig. 5.2. Chapter 5. A Comparison with Time Dependent Density Functional Theory 115 Figure 5.9: A comparison of particle solvation response functions obtained from simula-tion and nonlinear theory for System 4. Curves and panels are as in Fig. 5.2, although the Pr(t) curve has been omitted, as the statistical noise prevented any useful consideration. Chapter 5. A Comparison with Time Dependent Density Functional Theory 116 a relatively poor job of predicting the shapes of the simulation curves. Even the optimal choice oi D/Db obtained by fitting the curves for each system is not entirely satisfactory. This is most evident for System 3 [Fig. 5.4] where the theoretical and simulation curves differ in both shape and relative magnitude. Furthermore, we note that the values of D/Db obtained by fitting the curves from linearized T D D F T vary quite widely from 0.760 for System 1 to 3.295 for System 4. This indicates that the theory is rather poor because D/Db should depend only on the solvent and hence, in principle, should have the same value for all systems considered. Moreover, for Systems 2, 3 and 4 we obtain D/Db > 1 whereas physically we would expect D/Db < 1 since Db is associated with faster momentum decay processes. We note that using D/Db = 1 gives reasonable agreement with simulations for most systems considered (System 4 is an exception) and is always better than using the Araki-Munakata value. The particle response functions given by the nonlinear T D D F T are compared with simulations in Figs. 5.6-5.9. Here we see that if one chooses D/Db to yield the optimal fit, the curves given by the nonlinear theory fit the simulation results extremely well. Also, the values of D/Db obtained for the different systems vary only between 0.086 (System 4) and 0.152 (System 3) and in fact if one uses the average value, 0.11, for all systems the agreement between theory and simulation for the Pj(t) curves is altered only marginally. Furthermore, in this case the values of D/Db are much less than 1, which is as we would physically expect since Db is associated with fast momentum decay processes. Thus, the nonlinear T D D F T appears to give the correct shapes and relative magnitudes of the relaxation curves. We note that assuming D/Db = 1 gives poor results and clearly should not be used with the nonlinear theory. The Araki-Munakata value, D/Db — 0.648, is somewhat better but is still too large to give satisfactory agreement with simulation. An extension of the linearized T D D F T approach to Stockmayer solvents, such as those presented in Chapter 3 was also performed [70]. However, it is well-known that Chapter 5. A Comparison with Time Dependent Density Functional Theory 117 the form of the theory presented above is somewhat inadequate for these particular models. Most importantly, non-Markovian effects which play a significant role in the fast reorientational mode are completely ignored. This can be clearly seen in Fig. 5.10, which shows the particle solvation response and solvation response functions for System 2 of Chapter 3, along with the theoretical curves for the same system. The agreement between the two is rather poor, as has been noted previously for pure systems. A version of T D D F T which includes frequency-dependent diffusion constants has been shown to provide better results for these pure systems, and may do likewise for the mixed solvents presented here. An extension of the non-linear theory to these Stockmayer systems was not attempted, due to the complexity inherent in nonspherically symmetric potentials. 5.3 Summary In this chapter the dynamics of selective solvation in a simple binary solvent have been explored. The main purpose was to investigate the accuracy of both linear and nonlinear versions of T D D F T by comparing with M D simulations. The primary focus has been upon the relaxation of the size and composition of the first solvation shell of a solute particle, as monitored through the particle solvation response functions, Pj(i). A simple model was chosen to allow the T D D F T equations to be solved without further math-ematical approximation. Nevertheless, despite the simplicity of the model, one would fully expect the conclusions arrived at here to hold for more realistic systems. In the course of this work, a problem with the T D D F T approach has become apparent which hasn't received much previous attention. In the general nonlinear formulation of T D D F T the timescale is determined by a "bare diffusion constant", Db. Unfortunately, Db, is not known a priori and there does not appear to be any definition which allows its determination from experiments or simulations. Thus, expressing T D D F T results on an Chapter 5. A Comparison with Time Dependent Density Functional Theory 118 -0.5 1 1 1 1 1 ' 1 1 1 0 1 2 3 4 tDS Figure 5.10: A comparison of particle and solvation response functions obtained from theory and simulation for System 2 of Chapter 3. The solid and dashed curves are simulation and theoretical results, respectively. From top to bottom the curves are Pvv(t) [Ss(t)], Ps(t) [S(t)] and P(t) [SV(i)]. The timescale is plotted in units of the rotation diffusion constant for the weakly polar species, D^ . Chapter 5. A Comparison with Time Dependent Density Functional Theory 119 absolute timescale which can be compared with simulation or experiment is problematic. One possible approach is to find a relationship between the usual self-diffusion constant, D, and Db which, since D can be measured or calculated, would relate T D D F T results to an absolute timescale. However, at present this cannot be done exactly and one must use approximate relationships or resort to "fitting" the theory to the simulation curves. Both approaches have been investigated. In the linear case one can apply the Kerr approximation to obtain D = Db. This is the assumption always made at the outset in calculations using linearized T D D F T and previous calculations have shown that it works reasonably well for solvation in one-component solvents. However, it cannot be shown that D = Db is an exact relationship even in the linear context. For the mixed solvent systems considered here, it has been shown that linearized T D D F T using D/Db = 1 is qualitatively reasonable in predicting correct trends etc., but can be quantitatively very poor for some systems. Another approximate method of relating D and Db recently put forward by Araki and Munakata does not work as well as the Kerr result in the linear case. One can obtain better "fits" between the theoretical and simulation curves by treating D/Db as an adjustable parameter. However, in the linear case this is not very satisfactory because the values of D/Db obtained are obviously unphysical. In particular, the values from the fits show a strong solute dependence, whereas the ratio should depend only on the solvent. Also, the values obtained are often much greater than 1, whereas values less than 1 would be expected on physical grounds. Thus, with linearized T D D F T using D = Db appears to be the best option currently available. However, it should be emphasized that the linear theory does not give an accurate quantitative description of the shapes and relative magnitudes of the Pj(t) curves. The situation is very different for the nonlinear T D D F T . In this case, the assumption that D/Db = 1 gives extremely poor agreement between theory and simulation and it is Chapter 5. A Comparison with Time Dependent Density Functional Theory 120 obvious that the timescales have not been correctly matched. The theory of Araki and Munakata gives a somewhat smaller D/Db ratio and improves things to some extent but the agreement between theory and simulation remains very poor. On the other hand, excellent fits between the theoretical and simulation results can be obtained, showing that the nonlinear theory correctly describes the shapes and relative magnitudes of the Pj(t) curves. Furthermore, the values of D/Db obtained show much less solute dependence than in the linear case and indeed a single average value fits all systems very well. Also, the ratios obtained are much less than 1 as we would expect from an accurate theory. In summary, linearized T D D F T theory gives only a qualitatively accurate description of selective solvation dynamics. The nonlinear theory is much better and gives a very good description of the relaxation of the first solvation shell of the solute. However, the problem of matching theoretical and simulation or experimental timescales (without resorting to fitting procedures) remains unsolved and deserves further thought and investigation. Also, extension of linearized T D D F T in its present form to Stockmayer systems appears to be of limited usefulness due to weaknesses inherent in the formulation. Modifications to the theory will be necessary in order to investigate these systems. Chapter 6 Summary and Conclusions This thesis has been concerned with the study of solvation dynamics in binary sol-vents. The primary technique employed has been that of computer simulation. The initial results concerning the study of relatively simple Stockmayer species yielded im-portant insights into the solvation dynamics in these systems. Of particular interest to experimentalists is the fact that, even for systems with a large difference in solvent po-larities, the standard solvation response function is remarkably similar to that obtained for pure systems. This is true even in situations where a considerable amount of solvent redistribution is taking place. The reason for this is that the solvation response function measures the overall energetic relaxation of the system. The energy gained by the ex-change process, and that which must drive the preferential solvation itself, is small when compared to the total relaxation energy. As the solute-solvent energy is the quantity which is most directly related to the experimentally observed spectrum relaxation, this poses a considerable problem when studying mixed solvents. An alternative, the particle solvation response function P(t), was introduced as a means of studying the compositional relaxation of the first solvation shell of the solute. Application of this function to the Stockmayer systems revealed that the solvent relax-ation occurs in two phases following reorientation. An initial electrostriction phase taking place on short timescales results in an increase in density of both solvent species about the solute. This is followed by a considerably slower redistribution phase by which the solvation shells become richer in the preferred species. This redistribution, at least for 121 Chapter 6. Summary and Conclusions 122 the Stockmayer solvents studied, occurs at constant total coordination number, implying that the solvent species are being exchanged on a one-to-one basis. A simple phenomeno-logical model based on elementary chemical kinetics yielded expressions for P{t) which fit the simulation results extremely well. The application of the information gained from the Stockmayer systems to more realistic solvent models appears to depend greatly on the particular species being con-sidered. For the case of water-methanol, the picture is much the same. A n initial, rapid electrostriction is followed by a considerably slower redistribution phase in which a one-to-one exchange of water and methanol molecules occurs. The strong hydrogen-bonding nature of these solvents does not appear to have a significant effect on the timescale of importance here. Although the results for the water-methanol systems were strikingly similar to those of the Stockmayer systems, there were some differences. The most no-table of these was the greater degree of electrostriction evident in these systems. This is not entirely unexpected, however, as it is well-known that species such as water are somewhat loosely packed when compared to simpler models. This is due to the highly directional nature of the hydrogen bonds. The result is a dramatic increase in density about the solute following ionization as the water and methanol structure is compressed under the ionic field. In contrast to the water-methanol case, the water-DMSO systems exhibited behaviour which differed considerably from the Stockmayer results. Effects of both solvent size as well as geometry were seen to yield results which not only differed from the simpler mod-els, but also differed considerably for positively and negatively charged solutes within the same water-DMSO system. The nature of both the equilibrium structure as well as relax-ation behaviour in the negatively charged system was considerably different than what was observed for the Stockmayer and water-methanol systems. The fact that D M S O is a poor solvent for negatively charged solute resulted in a near exclusion of this species from Chapter 6. Summary and Conclusions 123 the first solvation shell of the solute. Thus, the relaxation consisted mainly of diffusion of water molecules towards the solute, through a medium of DMSO. In contrast, the results for the positively charged system bore much more resemblance to the Stockmayer results. In this case, DMSO acts as a strong solvent for the positively charged solute particle, once again leading to competitive solvation effects. Despite these differences from the simple Stockmayer results, the realistic systems further reinforced the observation that the standard solvation response function S(t) is only of moderate usefulness when studying solvation in mixed solvents. Again, the energetic contribution from the solvent redistribution was shown to be quite small. Thus, experimental studies in such systems may be indistinguishable from related pure solvents. Indeed, this has recently been observed experimentally by Nishiyama and Okada [39], who suggest that the relaxation of the spectral width may be a quantity which proves substantially more sensitive to the redistribution process. Finally, a comparison between computer simulation and time dependent density func-tional theory was made. Interestingly, under certain approximations, the solution to the linearized theoretical equations exhibited two timescales which could be easily identified as electrostriction and redistribution. However, this clear separation of timescales occurs in k space, and does not strictly hold upon translation back to r space. Nonetheless, the observation that the T D D F T results reproduce the simulation results quite well indicates that the r-space behaviour is still dominated by two timescales. It was also found that the common assumption that the bare diffusion constant is equal to the self-diffusion constant is only true subject to the Kerr approximation, and only yields good results for the linearized form of the theory. While the nonlinear theory was much better at pre-dicting the relative shapes of the simulation curves, no means is available for obtaining a good estimate of the bare diffusion constant. Chapter 6. Summary and Conclusions 124 A n application of the linearized theory to Stockmayer systems revealed that the the-ory, in its present form, is incapable of reproducing simulation results for these types of models. It is not obvious whether the nonlinear theory will properly capture the be-haviour observed in the simulation results, or whether an alternative approach, such as the inclusion of frequency dependent diffusion constants, is required. In addition to further insight into the nature of the phenomena itself, the comparison with simulation also yielded important information about the limitations of the theoret-ical approach. The equations themselves are formulated in terms of the "bare diffusion constant" which determines the timescale of the dynamics. Unfortunately, no di-rect physical means of determining this quantity is available, and many workers have made the implicit assumption that this quantity is in fact the self diffusion constant D. From the results presented in Chapter 5, it is clear that this is only an approximation. Additionally, the validity of this approximation appears to be satisfactory only in the linearized form of the theory. Bibliography [1] B. M . Ladanyi and M . Maroncelli. J. Chem. Phys., 109:3204, 1998. [2] J . Zhang, L . L . Lee, and J . F . Brennecke. J. Phys. Chem., 99:9268, 1995. [3] F. Cichos, A . Willert, U . Rempel, and C. von Borczyskowski. J. Phys. Chem., 101:8179, 1997. [4] B. J. Schwartz and P. J. Rossky. J. Phys. Chem., 99:2953, 1994. [5] G. Makov and A. Nitzan. J. Phys. Chem., 98:3459, 1994. [6] M . Maroncelli. J. Chem. Phys., 94:2084, 1991. [7] A . Chandra, D. Wei, and G. N.'Patey. J. Chem. Phys., 99:4926, 1993. [8] D. K . Phelps, M . J. Weaver, and B. M . Ladanyi. Chem. Phys., 176:575, 1993. [9] L. Perera and M . L. Berkowitz. J. Chem. Phys., 96:3092, 1992. [10] E. Neria and A . Nitzan. J. Chem. Phys., 96:5433, 1992. [11] A . Chandra. Chem. Phys. Lett, 235:133, 1995. [12] M . Maroncelli, V . P. Kumar, and A. Papazyan. J. Phys. Chem., 97:13, 1993. [13] M . S. Skaf and B. M . Ladanyi. J. Chem. Phys., 102:6542, 1995. [14] B. M . Ladanyi and R. M . Stratt. J. Phys. Chem., 100:1266, 1996. [15] B. Parbhoo and O. B. Nagy. J. Chem. Soc. Faraday Trans., 82:1789, 1986. [16] R. Jellema, J. Bulthuis, G. van der Zwan, and G. Somsen. J. Chem. Soc. Faraday Trans., 92:2569, 1996. [17] T. Morita, B . M . Ladanyi, and J . T. Hynes. J. Phys. Chem., 93:1286, 1989. [18] T. Fonseca and B. M . Ladanyi. J. Phys. Chem., 95:2116, 1991. [19] R. Biswas and B. Bagchi. J. Phys. Chem., 100:1238, 1995. [20] M . Cho, L. D. Ziegler S. J. Rosenthal, N . F. Scherer, and G. R. Fleming. J. Chem. Phys., 96:5033, 1992. 125 Bibliography 126 [21] A . Chandra and G. N . Patey. J. Chem. Phys., 100:1552, 1994. [22] R. F. Loring and S. Mukamel. J. Chem. Phys., 87:1272, 1987. [23] W. Jarzeba, G. C. Walker, A . E. Johnson, and P. F. Barbara. Chem. Phys., 152:57, 1991. [24] N . K h . Petrov, A . Wiessner, T. Fiebig, and H. Staerk. Chem. Phys. Lett, 241:127, 1995. [25] M . L. Horng, J. A . Gardecki, A . Papazyan, and M . Maroncelli. J. Phys. Chem., 99:17311, 1995. [26] G. R. Fleming and P. G. Wolynes. Physics Today, page 36, 1990. [27] M . Maroncelli, J. Maclnnis, and G. R. Fleming. Science, 243:1674, 1989. [28] G. W. Robinson, R. J. Robbins, G. R. Fleming, J. M . Morris, A . E. W. Knight, and R. J. S. Morrison. J. Am. Chem. Soc, 100:7145, 1978. [29] L . A . Hallidy and M . R. Topp. J. Phys. Chem., 82:2415, 1978. [30] M . P. Allen and D. J. Tildesly. Computer Simulation of Liquids. Clarendon Press, 1989. [31] H. Goldstein. Classical Mechanics. Addison-Wesley, 1980. [32] A . K . Soper and A. Luzar. J. Chem. Phys., 97:1320, 1992. [33] A . Sacco. Chem. Soc. Rev., 23:129, 1994. [34] A . K . Covington and M . Dunn. J. Chem. Soc. Faraday Trans., 85:2827, 1989. [35] J. P. Hansen and I. R. McDonald. Theory of Simple Liquids. Academic Press Limited, 1986. [36] C. J. F. Bottcher. Theory of Electric Polarization. American Elsevier Publishing Company, Inc., 1973. [37] W. Smith. CCP5 Quarterly, 4:13, 1982. [38] D. J. Evans. Mol. Phys., 34:317, 1977. [39] K . Nishiyama and T. Okada. J. Phys. Chem. A, 102:9729, 1998. [40] T. Fonseca and B. M . Ladanyi. J. Mol. Liquids, 60:1, 1994. Bibliography 127 [41] M . S. Skaf and B. M . Ladanyi. J. Phys. Chem., 100:18258, 1996. [42] P. G. Kusalik and G. N . Patey. J. Chem. Phys:, 89:5843, 1988. [43] H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma. J. Phys. Chem., 91:6269, 1987. [44] M . Haughney, M . Ferrario, and I. R. McDonald. J. Phys. Chem., 91:4934, 1987. [45] I. I. Vaisman and M . L. Berkowitz. J. Am. Chem. Soc., 114:7889, 1992. [46] J. Chandrasekhar, D. C. Spellmeyer, and W. L. Jorgensen. J. Am. Chem. Soc, 106:903, 1984. [47] R. T. Morrison and R. N . Boyd. Organic Chemistry. Allyn and Bacon, 1987. [48] B. Bagchi and A. Chandra. Adv. Chem. Phys., 80:1, 1991. [49] A . Chandra and B. Bagchi. J. Chem. Phys., 90:1832, 1989. [50] T. Munakata. J. Phys. Soc. of Japan, 43:1725, 1977. [51] T. Munakata. J. Phys. Soc. of Japan, 58:2434, 1989. [52] T. Munakata. J. Phys. Soc. of Japan, 59:1299, 1990. [53] T. Munakata. Phys. Rev. E, 50:2347, 1994. [54] J. Araki and T. Munakata. Phys. Rev. E, 52:2577, 1995. [55] A . Yoshimori. J. Chem. Phys., 105:5971, 1996. [56] D. F. Calef and P. G. Wolynes. J. Chem. Phys., 78:4145, 1983. [57] F. O. Raineri, H. Resat, B. -C. Perng, F. Hirata, and H. L. Friedman. J. Chem. Phys., 100:1477, 1994. [58] F. O. Raineri, B. -C. Perng, and H. L. Friedman. Chem. Phys., 183:187, 1994. [59] H. L . Friedman, F. O. Raineri, F. Hirata, and B.-C. Perng. J. Stat. Phys., 78:239, 1995. [60] B. Bagchi. J. Chem. Phys., 100:6658, 1994. [61] R. Biswas and B. Bagchi. J. Phys. Chem., 100:1238, 1996. [62] D. Wei and G. N . Patey. J. Chem. Phys., 91:7113, 1989. Bibliography 128 [63] D. Wei and G. N . Patey. J. Chem. Phys., 94:6785, 1991. [64] A . Chandra, D. Wei, and G. N . Patey. J. Chem. Phys., 99:2068, 1993. [65] A . Chandra, D. Wei, and G. N . Patey. J. Chem. Phys., 100:1552, 1994. [66] D. Wei and G. N . Patey. J. Chem. Phys., 93:1399, 1990. [67] S. Roy and B. Bagchi. J. Chem. Phys., 99:9938, 1993. [68] N . Nandi, S. Roy, and B. Bagchi. J. Chem. Phys., 102:1390, 1995. [69] H. J. K im, H. L. Friedman, and F. 0 . Raineri. J. Chem. Phys., 94:1442, 1991. [70] A . Yoshimori, T. J. F. Day, and G. N . Patey. J. Chem. Phys., 109:3222, 1998. [71] S. W. de Leeuw, J. W. Perram, and E. R. Smith. Proc. R. Soc. Lond., 373:27, 1980. [72] C. Kittel . Introduction to Solid State Physics. John Wiley and Sons, Inc., 1986. Appendix A The Potential Due to a Point Dipole For simplicity, consider two point charges, +q and —q, located on the z axis, equidis-tant from the origin at a distance of 1/2, as illustrated in Fig. A . l . Let P be a point in space with position vectors r + and r_ with respect to +o and —q, and r with respect to the origin. The potential at this point is simply the sum of the contributions from both point charges <j>(P) = <j>+{P) + <t>_{P), (A.l) r+ r_ where r± = | r ± | . Minor rearrangement yields r_r+ The potential due to a point dipole can then be obtained by investigating the case where r » I. Conversely, one might imagine replacing I with Is and q with q/s and investigating the limit s —> 0. Either process will yield the same two limiting results r + r _ —»• r 2 r _ - r + -»• I cos 6. (A.3) Consequently, the expression for the potential at point P reduces to where yu = ql is the dipole moment. This expression can be rewritten in the somewhat more useful form m = ^ (A5) 129 Appendix A. The Potential Due to a Point Dipole 130 where p = ql. The interaction of a second point dipole, p', located at P is thus given by UDD = (M-V)0(P) (A.6) p'p 3(/i'-r)(ji-r) y» 3 y»5 which is the same result as that obtained previously using the multipole expansion. Appendix A. The Potential Due to a Point Dipole 131 Figure A . l : Diagram for the derivation of the potential at a point P due to a point dipole moment. Appendix B The Ewald Sum A n infinite sum over an irregular lattice, such as is of interest in computer simulation is absolutely convergent only when the interaction potentials are relatively short-ranged. In the case of charge-charge and dipole-dipole interactions, the sum is only conditionally convergent. That is, the result is dependent on the particular order in which the lattice vectors are included in the summation [71]. The derivation presented here, and that which is most common, assumes that the cells are added on in spherical "shells" of increasing radii, much like the layers of an onion. It is also important to specify the boundary conditions of the system. At any given instant in time, the spherical collection of simulation cells will have a nonzero net dipole moment. Thus, it is necessary to specify the dielectric nature of the medium in which this infinite sphere is imbedded. The result is an extra term in the Ewald formulation which represents the contribution from the overall dipole moment. In the limit as e' —>• oo (i. e. a perfect conductor), this contribution vanishes. As it has been shown [71] that any particular choice of the dielectric constant has very little effect on systems such as those of interest here, it will be assumed that e' —>• oo in the following derivation. Consider an irregular lattice extending in three dimensions whose unit cell consists of N point charges qi with positions such that the cell itself is neutral. Now, suppose that each of these point charges is surrounded by a distributed charge of equal magnitude but opposite sign. The charge density, p, a given distance s from the original point charge may 132 Appendix B. The Ewald Sum 133 be of many forms. For the present work we will choose the familiar Gaussian distribution where a is an arbitrary parameter which determines the width of the distribution, and the leading coefficient ensures that the distributed charge is of the same magnitude as the original point charge. The addition of this diffuse charge screens the interactions of the original point charges, thus effectively reducing the range of the interaction. The exact range over which any significant contribution to the lattice energy is present depends on the particular choice of a. A large value for a will yield a much narrower distributed charge, thus leading to a very short-ranged interaction. A small value will have the opposite effect. The desired lattice energy corresponding to the original point charges can be regained by subsequently adding in a cancelling distribution of similar shape to the screening charge but of the same sign as the original point charge. Consequently, the potential, <f>, at any point r within the cell can be written as the sum of two contributions 4>{T) = (j)screen + 4>canc, .2) which correspond to the screened interactions of the point charges and the contribution from the cancelling distributions, respectively. If one expands both the potential, (j>, and charge density, p, in a Fourier series [72] 0(r) = £ c k e l k r , (B.3a) k p(r) = £ p k e ' k r , (B.3b) k where k are reciprocal lattice vectors, and relates the two using Poisson's equation, V 2^ = 4-7rp, one obtains * = % (BA) Appendix B. The Ewald Sum 134 To obtain an expression for the density coefficient, pk/, corresponding to a particular lattice vector k', one can multiply both sides of Eq. (B.3b) by exp — ik' • r and integrate over the unit cell / e- k 'Xr)dr = / £ p k e i k r e - i k ' r dr . (B.5) J cell J cell k which, due to the orthogonality of the Fourier terms, simplifies to / e-k'rp(r)dr = pk, f dr J cell J cell = PwV, (B.6) where V is the volume of the unit cell. Unfortunately, the integration on the right-hand side of Eq. (B.6) over the charge density p(r) requires that all contributions to the density in the primary cell be accounted for. This includes not only those charges located within the cell itself, but also the tails of the infinite number of charge distributions originating from outside the cell. However, due to the periodicity inherent in the system, the integration of the total charge density over the single cell can be replaced by an integration over all space considering only charges which originate in the single cell itself. Thus, Eq. (B.6) can be rewritten as [72] > .Qi — I allspace ( £ ^ e " i k ' ' r i ) e - * 2 a 2 / 4 . (B.7) PwV = I ^ q ^ e - ^ - ^ e - ^ ^ d r J , 7T ' Substituting into Eq. (B.4), one can then obtain an expression for Ck' C « V = 7^(E^e- k'- r0e- f c 2 Q 2 / 4 , (B.8) which yields the following expression for the contribution from the cancelling charges Knc(r) = f E E | e - ^ k - r - ^ / 4 ; ' , ( R 9 ) Appendix B. The Ewald Sum 135 or, if one considers the origin, so that r = 0 ^c(r) = y £ £ | e - ^ e - * 2 " 2 / 4 . (£.10) When considering the potential at one of the point charges itself, the above formulation includes the contribution from the cancelling distribution centered at the location being considered. This contribution, (f)seif, must be corrected for and can be calculated as simply an integral over the Gaussian charge distribution <f>self(v) = r^rp{v)dr = (B.ll) Jo V 7 1" The remaining contribution to the overall potential is that due to the original point charges along with their screening charge distributions. This is given by ^ee,(r) = E ^ r / c M , (B.12) i n where er /c is the complementary error function given by erfc(x) = -= \ e~s ds. (B.13) y/TT Jx Thus, the potential at a lattice point is given by Eqs. (B.10), ( B . l l ) , and (B.12) *M = v E E * * ' / f t - ' v " + - 2*2! (B.i4) k % i r v ^ The value of Eq. (B.14) must obviously be independent of any particular choice of a. However, a carefully chosen value can optimize the rate of convergence of the two sums. A typical choice is made such that only the charges in the primary simulation cell need be considered in the second term of Eq. (B.14). A n appropriately large number of reciprocal lattice vectors, k, are then included to guarantee the correct value for the potential is obtained. Appendix B. The Ewald Sum 136 The above treatment of point charges can be extended to point dipoles simply by replacing with (p{ • V ) in Eq. (B.14) to yield the potential at r due to a lattice of point dipoles [30] = y E E ( ^ - V ) e - ^ A 2 e - f e 2 a 2 / 4 k i + 2>, • V f l ^ l - 2ff. (B.15) where Hi = \pi\. A particular problem arises in relation to the work in this thesis. The above formulation of the Ewald method assumes that the unit cell from which the infinite lattice is generated carries no net charge. However, if one wishes to study the case where the solute particle has been ionized, this is clearly not true. Unfortunately, if each periodic image of the primary cell contains such an ion, the lattice as a whole will carry an infinite charge, and thus the lattice sum will be divergent. This problem can be circumvented by specifying that only the primary cell carries a net charge. The solute particles in the periodic images, while still present, are defined to be uncharged. Thus, the lattice as a whole will only carry a net charge of +le. As a result, the desired energy is simply the energy of insertion of a single point charge into the existing lattice. For the system of point charges, this is given by k i + £ ^ ^ ^ , (B.16) i r while for the dipolar system it is given by k i + E 9 ( M , - V ) ^ M . ( B . 1 7 ) Appendix B. The Ewald Sum 137 Since the point at which the field is being sampled, namely the location of the solute point charge, is no longer a regular lattice site, the correction term which was required in Eqs. (B.14) and (B.15) is no longer required in Eqs. (B.16) and (B.17). Appendix C The Reaction Field Method Consider a point dipole, p,, located at the center of a cavity surrounded by a polariz-able continuum. The field due to the dipole will polarize the surroundings, which will in turn produce a "reaction" field, E#, which acts on the original dipole. If it is assumed that the magnitude of the dipole moment is not unusually large and that the cavity is spherical, this field should be proportional to the dipole moment itself, ER OC p. (Cl) The proportionality constant can be determined from the solution of Laplace's equation, which, for this configuration, yields [36] _ 1 2 ( 6 - 1 ) E f l - ^ " 2 7 T i ~ ' 4 ' ( } where a is the radius of the cavity and e is the dielectric constant of the surrounding continuum. Thus, the energy of interaction of the original dipole with this reaction field is given simply by URF = » - ER. (C.3) However, it must be noted that this expression assumes that the surrounding dielectric is already polarized. If one desires to calculate the total interaction energy of the dipole with the surrounding medium, one must take into account the work required to polarize the dielectric continuum. The work, W\, required to bring a point dipole Hi to position r i in a cavity within a dielectric can be obtained by slowly increasing the magnitude of the dipole moment from 138 Appendix C. The Reaction Field Method 139 zero to its actual value in infinitesimally small steps. Thus, at any point the moment can be written as A ux where A ranges from 0 to 1. Similarly, the potential at r i due to the polarization of the dielectric can be written as Xcj)x (ri). Consequently, the work required to increase the moment by SX is given by SWX = Xtf (n) • pxSX, (CA) and thus, the work for the overall charging process is given by [36] Wx = jf 1 A0f (n) • pxSX = px • \ct>f(vx). (C.5) Consequently, the actual energy of the dipole within the cavity is only half of the inter-action between the dipole and the corresponding reaction field [Eq. (C.3)], and is given by V = 1 (C . 6 ) 2a3 2e+l K ' If a second dipole is brought into the cavity, the energy required will then consist of three contributions; the direct interaction with the first dipole, the interaction with the reaction field of the first dipole, and the interaction with its own reaction field. In a system consisting of N such point dipoles located within the cavity, the reaction field acting on any given molecule is the sum of contributions from all dipoles, including the one under consideration Due to the additive nature of this interaction, it can be reconstructed as a more familiar pairwise contribution to the overall interaction. Thus, the total energy for a system of point dipoles calculated using the reaction field method may be written as [36] u-ypu™ i 2 ( e ~ 1 } 1 Tu u i fo-1) 1 W (CS) \ i j + 2e + l atf" ^ + 2e + l a * \ ^ { C * ] where ufjD is the standard dipole-dipole interaction given by Eq. (2.2.12). Appendix D The Generalized Langevin Equation Presented here is a phenomenological derivation of the T D D F T equations employed in Chapter 5. A more rigourous derivation is also available [51,53], but is not essential to the work presented in this thesis. We begin with the familiar form for the generalized Langevin equation [35] where v(t) is the velocity of a tagged particle, m is its mass, and R(i) is a stochastic force acting on it. The integral, involving the "memory function" M(t), couples the behaviour of the system at time t to the history of the system prior to that time. If one makes the Markovian approximation, which assumes that the behaviour of the system is independent of its history and depends only on the current state, then the memory function can be replaced with where £ is often referred to as the "friction" coefficient. Substituting into Eq. (D.l) yields the Langevin equation M(t)=^5(t), (D.2) m dvjt) dt m£y{t) + R(r,*) (D.3) along with which the usual definition of the diffusion constant is made D = (DA) where kb is Boltzmann's constant and T is the temperature. 140 Appendix D. The Generalized Langevin Equation 141 If one now considers a single solute particle immersed in a solvent, then under equi-librium conditions, a given solvent particle at position r must be governed by a force balance equation m^jS. =B(r,*)-rmy(t) + R(r,t), (D.5) Cut / where the function B(r, t) is the force due to the direct interaction between the solute and the solvent particle, as well as effective interactions due to other factors such as the solvent structure. Due to the fact that B(r,t) contains effects which might normally be accounted for to some extent in the standard friction coefficient £, a different friction coefficient, T, has been introduced. If we now introduce the particle current j(r,*) = Ev i(J(r-r i(t)) (D.6) as well as the particle density p(r,t) = Yi5(r-vi(t)), (D.7) then Eq. (D.5) can be generalized to mJ-M = p ( r ; t ) B (r , t) - rmj(r, t) + R(r, t). (D.8) This form of the Langevin equation extends the original formulation to describe the diffusion of a number of tagged particles. If one makes the assumption that the velocities relax much more rapidly than does the density, one can set dt in Eq. (D .8 ) , which yields = 0 (D.9) j(r, t) = ^ p ( r , r)B(r, t) + p^R(r, t). (D.1Q) Appendix D. The Generalized Langevin Equation 142 If the usual continuity requirement P(r,t) = - V - j ( r , t ) (DAI) dt is then applied, the resulting expression becomes ^ = i v . W M ) B ( M ) ) . (D.12) As mentioned above, B ( r , t ) includes the direct solute-solvent interaction as well as effective interactions arising from structural considerations. Consequently, B ( r , t) may be related to the interaction part of the free energy functional, F[p(r,t)], via so that ^ M = ^V.fp(r , f )- V ^<f) 'V (D.U) dt Tm V ip(r,t) / Here it has been assumed that the time derivative of the free energy functional, which is an equilibrium quantity, may be employed. The "bare diffusion constant" Db, which is analogous to the self diffusion constant in the Langevin equation, is then defined as D i = k£ (ZX15) so that
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Ion solvation dynamics in binary solvents
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Ion solvation dynamics in binary solvents Day, Tyler James Frederick 1999
pdf
Page Metadata
Item Metadata
Title | Ion solvation dynamics in binary solvents |
Creator |
Day, Tyler James Frederick |
Date Issued | 1999 |
Description | The dynamics of selective ion solvation in binary solvents is investigated using molecular dynamics simulations. Initial studies of Stockmayer solvents focus on the dependence of the usual solvation response function, S(t), on solvent composition and on the relative polarity of the solvent species. Particle solvation response functions, P(t), are also introduced which describe the compositional relaxation of the first solvation shell. It is shown that the selective solvation process can be well described by a simple phenomenological model based on the ideas of elementary chemical kinetics. This model is useful and helps in the identification of two distinct timescales present in the selective solvation process. These are associated with a rapid electrostriction step during which the total number of particles in the first shell increases to its equilibrium value, and a slower spatial redistribution process during which the composition of the first shell achieves equilibrium. The redistribution phase depends on the rate of mutual ion-solvent diffusion and also on the rate of particle exchange between the first and second shells. A detailed analysis of the exchange process indicates that exchanges occur on virtually a one-to-one basis with the insertion of a stronger dipole into the first shell being mirrored by an almost immediate ejection of a weaker one. An extension to more realistic systems consisting of both water-methanol and waterdimethylsulfoxide (DMSO) mixtures was also carried out. It is found that, despite the presence of hydrogen-bonding, the dynamical behaviour in water-methanol systems is very similar to that observed for simple Stockmayer solvents. For the water-DMSO mixtures, however, it is found that the relative sizes and geometries of the solvent species can have a substantial influence on the preferential solvation process. For positively charged solutes in water-DMSO the physical picture does not differ greatly from the water-methanol case. However, for negative solutes the DMSO component shows only a weak response to the charge, and the solvation process consists largely of water molecules moving slowly to the solute through an essentially static DMSO medium. The results also illustrate that the usual solvation response function, which depends on the total solute-solvent energy, is not a very sensitive probe of selective solvation dynamics. This insensitivity has since been noticed in recent experimental studies, and an alternative function that appears to be more sensitive to the solvent composition near the solute has been suggested. A comparison of simulation results with time dependent density functional theory is also carried out. It is shown that the nonlinear version of the theory is necessary in order to obtain a good quantitative description of selective solvation dynamics. The linear theory is only of qualitative value. Also, attention is drawn to a previously little appreciated problem which arises when one attempts to compare time dependent density functional theory with computer simulation or experimental results. The difficulty involves matching the theoretical and absolute timescales. Finally, an extension of the theory to Stockmayer systems reveals that the theory in its present form is insufficient to capture the behaviour of these systems. |
Extent | 6932033 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-07-02 |
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. |
DOI | 10.14288/1.0061536 |
URI | http://hdl.handle.net/2429/9960 |
Degree |
Doctor of Philosophy - PhD |
Program |
Chemistry |
Affiliation |
Science, Faculty of Chemistry, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1999-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1999-463389.pdf [ 6.61MB ]
- Metadata
- JSON: 831-1.0061536.json
- JSON-LD: 831-1.0061536-ld.json
- RDF/XML (Pretty): 831-1.0061536-rdf.xml
- RDF/JSON: 831-1.0061536-rdf.json
- Turtle: 831-1.0061536-turtle.txt
- N-Triples: 831-1.0061536-rdf-ntriples.txt
- Original Record: 831-1.0061536-source.json
- Full Text
- 831-1.0061536-fulltext.txt
- Citation
- 831-1.0061536.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-0061536/manifest