--- trunk/chrisDissertation/Introduction.tex 2006/09/26 01:30:32 3022 +++ trunk/chrisDissertation/Introduction.tex 2006/09/26 03:07:59 3023 @@ -1,5 +1,61 @@ -\chapter{\label{chap:intro}INTRODUCTION AND BACKGROUND} +\chapter{\label{chap:intro}INTRODUCTION AND BACKGROUND} +The following dissertation presents the primary aspects of the +research I have performed and been involved with over the last several +years. Rather than presenting the topics in a chronological fashion, +they were arranged to form a series where the later topics apply and +extend the findings of the former topics. This layout does lead to +occasional situations where knowledge gleaned from earlier chapters +(particularly chapter \ref{chap:electrostatics}) is not applied +outright in the later work; however, I feel that this organization is +more instructive and provides a more cohesive progression of research +efforts. + +This chapter gives a general overview of molecular simulations, with +particular emphasis on considerations that need to be made in order to +apply the technique properly. This leads quite naturally into chapter +\ref{chap:electrostatics}, where we investigate correction techniques +for proper handling of long-ranged electrostatic interactions. In +particular we develop and evaluate some new simple pairwise +methods. These techniques make an appearance in the later chapters, as +they are applied to specific systems of interest, showing how it they +can improve the quality of various molecular simulations. + +In chapter \ref{chap:water}, we focus on simple water models, +specifically the single-point soft sticky dipole (SSD) family of water +models. These single-point models are more efficient than the common +multi-point partial charge models and, in many cases, better capture +the dynamic properties of water. We discuss improvements to these +models in regards to long-range electrostatic corrections and show +that these models can work well with the techniques discussed in +chapter \ref{chap:electrostatics}. By investigating and improving +simple water models, we are extending the limits of computational +efficiency for systems that depend heavily on water calculations. + +In chapter \ref{chap:ice}, we study a unique polymorph of ice that we +discovered while performing water simulations with the fast simple +water models discussed in the previous chapter. This form of ice, +which we called ``imaginary ice'' (Ice-$i$), has a low-density +structure which is different from any known polymorph from either +experiment or other simulations. In this study, we perform a free +energy analysis and see that this structure is in fact the +thermodynamically preferred form of ice for both the single-point and +commonly used multi-point water models under the chosen simulation +conditions. We also consider electrostatic corrections, again +including the techniques discussed in chapter +\ref{chap:electrostatics}, to see how they influence the free energy +results. This work, to some degree, addresses the appropriateness of +using these simplistic water models outside of the conditions for +which they were developed. + +Finally, in chapter \ref{chap:conclusion}, we summarize the work +presented in the previous chapters and connect ideas together with +some final comments. The supporting information follows in the +appendix, and it gives a more detailed look at systems discussed in +chapter \ref{chap:electrostatics}. + +\section{On Molecular Simulation} + In order to advance our understanding of natural chemical and physical processes, researchers develop explanations for events observed experimentally. These hypotheses, supported by a body of corroborating @@ -35,10 +91,9 @@ molecular dynamics near exclusively, so we will presen molecular dynamics can be used to investigate dynamical quantities. The research presented in the following chapters utilized molecular dynamics near exclusively, so we will present a general -introduction to molecular dynamics and not Monte Carlo. There are -several resources available for those desiring a more in-depth -presentation of either of these -techniques.\cite{Allen87,Frenkel02,Leach01} +introduction to molecular dynamics. There are several resources +available for those desiring a more in-depth presentation of either of +these techniques.\cite{Allen87,Frenkel02,Leach01} \section{\label{sec:MolecularDynamics}Molecular Dynamics} @@ -48,8 +103,8 @@ given configuration of particles at time $t_1$, basica configuration to evolve in a manner that mimics real molecular systems. To do this, we start with clarifying what we know about a given configuration of particles at time $t_1$, basically that each -particle in the configuration has a position ($q$) and velocity -($\dot{q}$). We now want to know what the configuration will be at +particle in the configuration has a position ($\mathbf{q}$) and velocity +($\dot{\mathbf{q}}$). We now want to know what the configuration will be at time $t_2$. To find out, we need the classical equations of motion, and one useful formulation of them is the Lagrangian form. @@ -65,25 +120,31 @@ over time is zero,\cite{Tolman38} motion, to say that the variation of the integral of the Lagrangian over time is zero,\cite{Tolman38} \begin{equation} -\delta\int_{t_1}^{t_2}L(q,\dot{q})dt = 0. +\delta\int_{t_1}^{t_2}L(\mathbf{q},\dot{\mathbf{q}})dt = 0. \end{equation} -This can be written as a summation of integrals to give +The variation can be transferred to the variables that make up the +Lagrangian, \begin{equation} -\int_{t_1}^{t_2}\sum_{i=1}^{3N}\left(\frac{\partial L}{\partial q_i}\delta q_i - + \frac{\partial L}{\partial \dot{q}_i}\delta \dot{q}_i\right)dt = 0. +\int_{t_1}^{t_2}\sum_{i=1}^{3N}\left( + \frac{\partial L}{\partial \mathbf{q}_i}\delta \mathbf{q}_i + + \frac{\partial L}{\partial \dot{\mathbf{q}}_i}\delta + \dot{\mathbf{q}}_i\right)dt = 0. \end{equation} -Using fact that $\dot q$ is the derivative with respect to time of $q$ -and integrating the second partial derivative in the parenthesis by -parts, this equation simplifies to +Using the fact that $\dot{\mathbf{q}}$ is the derivative of +$\mathbf{q}$ with respect to time and integrating the second partial +derivative in the parenthesis by parts, this equation simplifies to \begin{equation} \int_{t_1}^{t_2}\sum_{i=1}^{3N}\left( - \frac{d}{dt}\frac{\partial L}{\partial \dot{q}_i} - - \frac{\partial L}{\partial q_i}\right)\delta {q}_i dt = 0, + \frac{d}{dt}\frac{\partial L}{\partial \dot{\mathbf{q}}_i} + - \frac{\partial L}{\partial \mathbf{q}_i}\right) + \delta {\mathbf{q}}_i dt = 0, \end{equation} -and the above equation is only valid for +and since each variable is independent, we can separate the +contribution from each of the variables: \begin{equation} -\frac{d}{dt}\frac{\partial L}{\partial \dot{q}_i} - - \frac{\partial L}{\partial q_i} = 0\quad\quad(i=1,2,\dots,3N). +\frac{d}{dt}\frac{\partial L}{\partial \dot{\mathbf{q}}_i} + - \frac{\partial L}{\partial \mathbf{q}_i} = 0 + \quad\quad(i=1,2,\dots,3N). \label{eq:formulation} \end{equation} To obtain the classical equations of motion for the particles, we can @@ -126,8 +187,8 @@ differential equation. The velocities are necessary fo interesting to note that equation (\ref{eq:verlet}) does not include velocities, and this makes sense since they are not present in the differential equation. The velocities are necessary for the -calculation of the kinetic energy, and they can be calculated at each -time step with the following equation: +calculation of the kinetic energy and can be calculated at each time +step with the equation: \begin{equation} \mathbf{v}_i(t) = \frac{\mathbf{r}_i(t+\delta t)-\mathbf{r}_i(t-\delta t)} {2\delta t}. @@ -141,14 +202,14 @@ symplectic, meaning that it preserves area in phase-sp higher order information that is discarded after integrating steps. Another interesting property of this algorithm is that it is symplectic, meaning that it preserves area in phase-space. Symplectic -algorithms system stays in the region of phase-space dictated by the -ensemble and enjoy greater long-time energy -conservations.\cite{Frenkel02} +algorithms keep the system evolving in the region of phase-space +dictated by the ensemble and enjoy a greater degree of energy +conservation.\cite{Frenkel02} While the error in the positions calculated using the Verlet algorithm is small ($\mathcal{O}(\delta t^4)$), the error in the velocities is -quite a bit larger ($\mathcal{O}(\delta t^2)$).\cite{Allen87} Swope -{\it et al.} developed a corrected for of this algorithm, the +substantially larger ($\mathcal{O}(\delta t^2)$).\cite{Allen87} Swope +{\it et al.} developed a corrected version of this algorithm, the `velocity Verlet' algorithm, which improves the error of the velocity calculation and thus the energy conservation.\cite{Swope82} This algorithm involves a full step of the positions, @@ -161,8 +222,8 @@ and a half step of the velocities, \mathbf{v}\left(t+\frac{1}{2}\delta t\right) = \mathbf{v}(t) + \frac{1}{2}\delta t\mathbf{a}(t). \end{equation} -After calculating new forces, the velocities can be updated to a full -step, +After forces are calculated at the new positions, the velocities can +be updated to a full step, \begin{equation} \mathbf{v}(t+\delta t) = \mathbf{v}\left(t+\frac{1}{2}\delta t\right) + \frac{1}{2}\delta t\mathbf{a}(t+\delta t). @@ -230,7 +291,7 @@ $\mathsf{A}$ can be expressed as arithmetic operations representation of the orientation of a rigid body.\cite{Evans77,Evans77b,Allen87} Thus, the elements of $\mathsf{A}$ can be expressed as arithmetic operations involving the -four quaternions ($q_w, q_x, q_y,$ and $q_z$), +four quaternions ($q_0, q_1, q_2,$ and $q_3$), \begin{equation} \mathsf{A} = \left( \begin{array}{l@{\quad}l@{\quad}l} q_0^2+q_1^2-q_2^2-q_3^2 & 2(q_1q_2+q_0q_3) & 2(q_1q_3-q_0q_2) \\ @@ -241,11 +302,11 @@ performance enhancements over Euler angles, particular Integration of the equations of motion involves a series of arithmetic operations involving the quaternions and angular momenta and leads to performance enhancements over Euler angles, particularly for very -small systems.\cite{Evans77} This integration methods works well for +small systems.\cite{Evans77} This integration method works well for propagating orientational motion in the canonical ensemble ($NVT$); however, energy conservation concerns arise when using the simple quaternion technique under the microcanonical ($NVE$) ensemble. An -earlier implementation of {\sc oopse} utilized quaternions for +earlier implementation of our simulation code utilized quaternions for propagation of rotational motion; however, a detailed investigation showed that they resulted in a steady drift in the total energy, something that has been observed by Kol {\it et al.} (also see @@ -334,7 +395,7 @@ propagator.\cite{Trotter59} This has three effects: {\it symplectic}), \item the integrator is time-{\it reversible}, and \item the error for a single time step is of order -$\mathcal{O}\left(\delta t^4\right)$ for time steps of length $\delta t$. +$\mathcal{O}\left(\delta t^3\right)$ for time steps of length $\delta t$. \end{enumerate} After the initial half-step ({\tt moveA}), the forces and torques are @@ -421,8 +482,8 @@ integrator, while the quaternion scheme will require ~ 0.001~kcal~mol$^{-1}$ per particle is desired, a nanosecond of simulation time will require ~19 hours of CPU time with the {\sc dlm} integrator, while the quaternion scheme will require ~154 hours of CPU -time. This demonstrates the computational advantage of the integration -scheme utilized in {\sc oopse}. +time. This demonstrates the computational advantage of the {\sc dlm} +integration scheme. \section{Accumulating Interactions} @@ -432,12 +493,12 @@ interactions increases by $\mathcal{O}(N(N-1)/2)$ if y multipole) on each particle from their surroundings. This can quickly become a cumbersome task for large systems since the number of pair interactions increases by $\mathcal{O}(N(N-1)/2)$ if you accumulate -interactions between all particles in the system, note the utilization -of Newton's third law to reduce the interaction number from -$\mathcal{O}(N^2)$. The case of periodic boundary conditions further -complicates matters by turning the finite system into an infinitely -repeating one. Fortunately, we can reduce the scale of this problem by -using spherical cutoff methods. +interactions between all particles in the system. (Note the +utilization of Newton's third law to reduce the interaction number +from $\mathcal{O}(N^2)$.) The case of periodic boundary conditions +further complicates matters by turning the finite system into an +infinitely repeating one. Fortunately, we can reduce the scale of this +problem by using spherical cutoff methods. \begin{figure} \centering @@ -450,8 +511,8 @@ interactions between all particles in the simulation, \end{figure} With spherical cutoffs, rather than accumulating the full set of interactions between all particles in the simulation, we only -explicitly consider interactions between local particles out to a -specified cutoff radius distance, $R_\textrm{c}$, (see figure +explicitly consider interactions between particles separated by less +than a specified cutoff radius distance, $R_\textrm{c}$, (see figure \ref{fig:sphereCut}). This reduces the scaling of the interaction to $\mathcal{O}(N\cdot\textrm{c})$, where `c' is a value that depends on the size of $R_\textrm{c}$ (c $\approx R_\textrm{c}^3$). Determination @@ -460,17 +521,17 @@ list radius $R_\textrm{l}$, which is larger than $R_\t the this expense, we can use neighbor lists.\cite{Verlet67,Thompson83} With neighbor lists, we have a second list of particles built from a list radius $R_\textrm{l}$, which is larger than $R_\textrm{c}$. Once -any of the particles within $R_\textrm{l}$ move a distance of +any particle within $R_\textrm{l}$ moves half the distance of $R_\textrm{l}-R_\textrm{c}$ (the ``skin'' thickness), we rebuild the list with the full $N(N-1)/2$ loop.\cite{Verlet67} With an appropriate -skin thickness, these updates are only performed every $\sim$20 -time steps, significantly reducing the time spent on pair-list -bookkeeping operations.\cite{Allen87} If these neighbor lists are -utilized, it is important that these list updates occur -regularly. Incorrect application of this technique leads to -non-physical dynamics, such as the ``flying block of ice'' behavior -for which improper neighbor list handling was identified a one of the -possible causes.\cite{Harvey98,Sagui99} +skin thickness, these updates are only performed every $\sim$20 time +steps, significantly reducing the time spent on pair-list bookkeeping +operations.\cite{Allen87} If these neighbor lists are utilized, it is +important that these list updates occur regularly. Incorrect +application of this technique leads to non-physical dynamics, such as +the ``flying block of ice'' behavior for which improper neighbor list +handling was identified a one of the possible +causes.\cite{Harvey98,Sagui99} \subsection{Correcting Cutoff Discontinuities} \begin{figure} @@ -485,15 +546,16 @@ region in order to smooth out the discontinuity.} region in order to smooth out the discontinuity.} \label{fig:ljCutoff} \end{figure} -As particle move in and out of $R_\textrm{c}$, there will be sudden -discontinuous jumps in the potential (and forces) due to their -appearance and disappearance. In order to prevent heating and poor -energy conservation in the simulations, this discontinuity needs to be -smoothed out. There are several ways to modify the function so that it -crosses $R_\textrm{c}$ in a continuous fashion, and the easiest -methods is shifting the potential. To shift the potential, we simply -subtract out the value we calculate at $R_\textrm{c}$ from the whole -potential. For the shifted form of the Lennard-Jones potential is +As the distance between a pair of particles fluctuates around +$R_\textrm{c}$, there will be sudden discontinuous jumps in the +potential (and forces) due to their inclusion and exclusion from the +interaction loop. In order to prevent heating and poor energy +conservation in the simulations, this discontinuity needs to be +smoothed out. There are several ways to modify the potential function +to eliminate these discontinuties, and the easiest methods is shifting +the potential. To shift the potential, we simply subtract out the +value we calculate at $R_\textrm{c}$ from the whole potential. The +shifted form of the Lennard-Jones potential is \begin{equation} V_\textrm{SLJ} = \left\{\begin{array}{l@{\quad\quad}l} V_\textrm{LJ}(r_{ij}) - V_\textrm{LJ}(R_\textrm{c}) & r_{ij} < R_\textrm{c}, \\ @@ -502,15 +564,16 @@ where \end{equation} where \begin{equation} -V_\textrm{LJ} = 4\epsilon\left[\left(\frac{\sigma}{r_{ij}}\right)^{12} - - \left(\frac{\sigma}{r_{ij}}\right)^6\right]. +V_\textrm{LJ}(r_{ij}) = + 4\epsilon\left[\left(\frac{\sigma}{r_{ij}}\right)^{12} - + \left(\frac{\sigma}{r_{ij}}\right)^6\right]. \end{equation} In figure \ref{fig:ljCutoff}, the shifted form of the potential reaches zero at the cutoff radius at the expense of the correct -magnitude below $R_\textrm{c}$. This correction method also does +magnitude inside $R_\textrm{c}$. This correction method also does nothing to correct the discontinuity in the forces. We can account for this force discontinuity by shifting in the by applying the shifting -in the forces rather than just the potential via +in the forces as well as in the potential via \begin{equation} V_\textrm{SFLJ} = \left\{\begin{array}{l@{\quad\quad}l} V_\textrm{LJ}({r_{ij}}) - V_\textrm{LJ}(R_\textrm{c}) - @@ -544,25 +607,25 @@ switching transition. \subsection{\label{sec:LJCorrections}Long-Range Interaction Corrections} -While a good approximation, accumulating interaction only from the -nearby particles can lead to some issues, because the further away -surrounding particles do still have a small effect. For instance, -while the strength of the Lennard-Jones interaction is quite weak at -$R_\textrm{c}$ in figure \ref{fig:ljCutoff}, we are discarding all of -the attractive interaction that extends out to extremely long +While a good approximation, accumulating interactions only from nearby +particles can lead to some issues, because particles at distances +greater than $R_\textrm{c}$ do still have a small effect. For +instance, while the strength of the Lennard-Jones interaction is quite +weak at $R_\textrm{c}$ in figure \ref{fig:ljCutoff}, we are discarding +all of the attractive interactions that extend out to extremely long distances. Thus, the potential is a little too high and the pressure on the central particle from the surroundings is a little too low. For -homogeneous Lennard-Jones systems, we can correct for this neglect by +homogeneous Lennard-Jones systems, we can correct for this effect by assuming a uniform density and integrating the missing part, \begin{equation} V_\textrm{full}(r_{ij}) \approx V_\textrm{c} + 2\pi N\rho\int^\infty_{R_\textrm{c}}r^2V_\textrm{LJ}(r)dr, \end{equation} where $V_\textrm{c}$ is the truncated Lennard-Jones -potential.\cite{Allen87} Like the potential, other Lennard-Jones -properties can be corrected by integration over the relevant -function. Note that with heterogeneous systems, this correction begins -to break down because the density is no longer uniform. +potential.\cite{Allen87} Like the potential, other properties can be +corrected by integration over the relevant function. Note that with +heterogeneous systems, this correction breaks down because the density +is no longer uniform. Correcting long-range electrostatic interactions is a topic of great importance in the field of molecular simulations. There have been @@ -607,14 +670,14 @@ unrestricted fashion while calculating interactions us region in space. To avoid this ``soft'' confinement, it is common practice to allow the particle coordinates to expand in an unrestricted fashion while calculating interactions using a wrapped -set of ``minimum image'' coordinates. These minimum image coordinates -need not be stored because they are easily calculated on the fly when -determining particle distances. +set of ``minimum image'' coordinates. These coordinates need not be +stored because they are easily calculated while determining particle +distances. \section{Calculating Properties} -In order to use simulations to model experimental process and evaluate -theories, we need to be able to extract properties from the +In order to use simulations to model experimental processes and +evaluate theories, we need to be able to extract properties from the results. In experiments, we can measure thermodynamic properties such as the pressure and free energy. In computer simulations, we can calculate properties from the motion and configuration of particles in @@ -623,17 +686,17 @@ isobaric-isothermal ($NPT$), and microcanonical ($NVE$ The work presented in the later chapters use the canonical ($NVT$), isobaric-isothermal ($NPT$), and microcanonical ($NVE$) statistical -mechanics ensembles. The different ensembles lend themselves to more +mechanical ensembles. The different ensembles lend themselves to more effectively calculating specific properties. For instance, if we concerned ourselves with the calculation of dynamic properties, which are dependent upon the motion of the particles, it is better to choose -an ensemble that does not add system motions to keep properties like -the temperature or pressure constant. In this case, the $NVE$ ensemble -would be the most appropriate choice. In chapter \ref{chap:ice}, we -discuss calculating free energies, which are non-mechanical -thermodynamic properties, and these calculations also depend on the -chosen ensemble.\cite{Allen87} The Helmholtz free energy ($A$) depends -on the $NVT$ partition function ($Q_{NVT}$), +an ensemble that does not add artificial motions to keep properties +like the temperature or pressure constant. In this case, the $NVE$ +ensemble would be the most appropriate choice. In chapter +\ref{chap:ice}, we discuss calculating free energies, which are +non-mechanical thermodynamic properties, and these calculations also +depend on the chosen ensemble.\cite{Allen87} The Helmholtz free energy +($A$) depends on the $NVT$ partition function ($Q_{NVT}$), \begin{equation} A = -k_\textrm{B}T\ln Q_{NVT}, \end{equation} @@ -652,25 +715,27 @@ calculated from an average over the densities for the chosen property like the density in the $NPT$ ensemble, where the volume is allowed to fluctuate. The density for the simulation is calculated from an average over the densities for the individual -configurations. Being that the configurations from the integration -process are related to one another by the time evolution of the -interactions, this average is technically a time average. In -calculating thermodynamic properties, we would actually prefer an -ensemble average that is representative of all accessible states of -the system. We can calculate thermodynamic properties from the time -average by taking advantage of the ergodic hypothesis, which states -that over a long period of time, the time average and the ensemble -average are the same. +configurations. Since the configurations from the integration process +are related to one another by the time evolution of the interactions, +this average is technically a time average. In calculating +thermodynamic properties, we would actually prefer an ensemble average +that is representative of all accessible states of the system. We can +calculate thermodynamic properties from the time average by taking +advantage of the ergodic hypothesis, which states that for a +sufficiently chaotic system, and over a long enough period of time, +the time and ensemble averages are the same. -In addition to the average, the fluctuations of that particular -property can be determined via the standard deviation. Fluctuations -are useful for measuring various thermodynamic properties in computer +In addition to the average, the fluctuations of a particular property +can be determined via the standard deviation. Not only are +fluctuations useful for determining the spread of values around the +average and the error in the calculation of the value, they are also +useful for measuring various thermodynamic properties in computer simulations. In section \ref{sec:t5peThermo}, we use fluctuations in -properties like the enthalpy and volume to calculate various +properties like the enthalpy and volume to calculate other thermodynamic properties, such as the constant pressure heat capacity and the isothermal compressibility. -\section{Application of the Techniques} +\section{OOPSE} In the following chapters, the above techniques are all utilized in the study of molecular systems. There are a number of excellent @@ -678,50 +743,17 @@ Because of our interest in rigid body dynamics, point incorporate many of these methods.\cite{Brooks83,MacKerell98,Pearlman95,Berendsen95,Lindahl01,Smith96,Ponder87} Because of our interest in rigid body dynamics, point multipoles, and -systems where the orientational degrees cannot be handled using the -common constraint procedures (like {\sc shake}), the majority of the -following work was performed using {\sc oopse}, the object-oriented -parallel simulation engine.\cite{Meineke05} The {\sc oopse} package -started out as a collection of separate programs written within our -group, and has developed into one of the few parallel molecular -dynamics packages capable of accurately integrating rigid bodies and -point multipoles. +systems where the orientational degrees of freedom cannot be handled +using the common constraint procedures (like {\sc shake}), the +majority of the following work was performed using {\sc oopse}, the +object-oriented parallel simulation engine.\cite{Meineke05} The {\sc +oopse} package started out as a collection of separate programs +written within our group, and has developed into one of the few +parallel molecular dynamics packages capable of accurately integrating +rigid bodies and point multipoles. This simulation package is +open-source software, available to anyone interested in performing +molecular dynamics simulations. More information about {\sc oopse} can +be found in reference \cite{Meineke05} or at the {\tt +http://oopse.org} website. -In chapter \ref{chap:electrostatics}, we investigate correction -techniques for proper handling of long-ranged electrostatic -interactions. In particular we develop and evaluate some new pairwise -methods which we have incorporated into {\sc oopse}. These techniques -make an appearance in the later chapters, as they are applied to -specific systems of interest, showing how it they can improve the -quality of various molecular simulations. -In chapter \ref{chap:water}, we focus on simple water models, -specifically the single-point soft sticky dipole (SSD) family of water -models. These single-point models are more efficient than the common -multi-point partial charge models and, in many cases, better capture -the dynamic properties of water. We discuss improvements to these -models in regards to long-range electrostatic corrections and show -that these models can work well with the techniques discussed in -chapter \ref{chap:electrostatics}. By investigating and improving -simple water models, we are extending the limits of computational -efficiency for systems that heavily depend on water calculations. - -In chapter \ref{chap:ice}, we study a unique polymorph of ice that we -discovered while performing water simulations with the fast simple -water models discussed in the previous chapter. This form of ice, -which we called ``imaginary ice'' (Ice-$i$), has a low-density -structure which is different from any known polymorph from either -experiment or other simulations. In this study, we perform a free -energy analysis and see that this structure is in fact the -thermodynamically preferred form of ice for both the single-point and -commonly used multi-point water models under the chosen simulation -conditions. We also consider electrostatic corrections, again -including the techniques discussed in chapter -\ref{chap:electrostatics}, to see how they influence the free energy -results. This work, to some degree, addresses the appropriateness of -using these simplistic water models outside of the conditions for -which they were developed. - -In the final chapter we summarize the work presented previously. We -also close with some final comments before the supporting information -presented in the appendices.