ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/chrisDissertation/Water.tex
(Generate patch)

Comparing trunk/chrisDissertation/Water.tex (file contents):
Revision 2987 by chrisfen, Wed Aug 30 22:36:06 2006 UTC vs.
Revision 3028 by chrisfen, Tue Sep 26 23:15:24 2006 UTC

# Line 5 | Line 5 | performed involve interactions with or between solvent
5   molecules of interest.  In some cases (such as in the simulation of
6   phospholipid bilayers), the majority of the calculations that are
7   performed involve interactions with or between solvent molecules.
8 < Thus, the motion and behavior of molecules in biochemical simulations
9 < are highly dependent on the properties of the water model that is
10 < chosen as the solvent.
8 > Thus, the motion and behavior of solute molecules in biochemical
9 > simulations are highly dependent on the properties of the water model
10 > that is chosen as the solvent.
11   \begin{figure}
12   \includegraphics[width=\linewidth]{./figures/waterModels.pdf}
13   \caption{Partial-charge geometries for the TIP3P, TIP4P, TIP5P, and
# Line 19 | Line 19 | interactions that need to be calculated.}
19   \label{fig:waterModels}
20   \end{figure}
21  
22 < As discussed in the previous chapter, water it typically modeled with
22 > As discussed in the previous chapter, water is typically modeled with
23   fixed geometries of point charges shielded by the repulsive part of a
24   Lennard-Jones interaction. Some of the common water models are shown
25   in figure \ref{fig:waterModels}. The various models all have their
# Line 31 | Line 31 | come from the extra point charges in the model descrip
31   computational cost.\cite{Mahoney00,Mahoney01} This cost is entirely
32   due to the additional distance and electrostatic calculations that
33   come from the extra point charges in the model description. Thus, the
34 < criteria for choosing a water model are capturing the liquid state
35 < properties and having the fewest number of points to insure efficient
36 < performance. As researchers have begun to study larger systems, such
37 < as entire viruses, the choice readily shifts towards efficiency over
38 < accuracy in order to make the calculations
39 < feasible.\cite{Freddolino06}
34 > main criteria for choosing a water model are:
35  
36 + \begin{enumerate}
37 + \item capturing the liquid state properties and
38 + \item having the fewest number of points to insure efficient performance.
39 + \end{enumerate}
40 + As researchers have begun to study larger systems, such as entire
41 + viruses, the choice has shifted towards efficiency over accuracy in
42 + order to make the calculations feasible.\cite{Freddolino06}
43 +
44   \section{Soft Sticky Dipole Model for Water}
45  
46 < One recently developed model that largely succeeds in retaining the
46 > One recently-developed model that largely succeeds in retaining the
47   accuracy of bulk properties while greatly reducing the computational
48   cost is the Soft Sticky Dipole (SSD) water
49   model.\cite{Liu96,Liu96b,Chandra99,Tan03} The SSD model was developed
# Line 50 | Line 53 | local structure in liquid water, a short ranged ``stic
53   Lennard-Jones core.  However, since the normal aligned and
54   anti-aligned geometries favored by point dipoles are poor mimics of
55   local structure in liquid water, a short ranged ``sticky'' potential
56 < is also added.  The sticky potential directs the molecules to assume
56 > is also added.  This sticky potential directs the molecules to assume
57   the proper hydrogen bond orientation in the first solvation shell.
58  
59   The interaction between two SSD water molecules \emph{i} and \emph{j}
# Line 133 | Line 136 | properties and behavior under the more computationally
136   electrostatics.\cite{Ewald21} In applying this water model in a
137   variety of molecular systems, it would be useful to know its
138   properties and behavior under the more computationally efficient
139 < reaction field (RF) technique, the correction techniques discussed in
140 < the previous chapter, or even a simple
141 < cutoff.\cite{Onsager36,Fennell06} This study addresses these issues by
142 < looking at the structural and transport behavior of SSD over a variety
143 < of temperatures with the purpose of utilizing the RF correction
144 < technique.  We then suggest modifications to the parameters that
145 < result in more realistic bulk phase behavior.  It should be noted that
146 < in a recent publication, some of the original investigators of the SSD
147 < water model have suggested adjustments to the SSD water model to
148 < address abnormal density behavior (also observed here), calling the
149 < corrected model SSD1.\cite{Tan03} In the later sections of this
150 < chapter, we compare our modified variants of SSD with both the
151 < original SSD and SSD1 models and discuss how our changes improve the
149 < depiction of water.
139 > reaction field (RF) technique (and later, the {\sc sf} method), the
140 > correction techniques discussed in the previous chapter, or even a
141 > simple cutoff.\cite{Onsager36,Fennell06} This chapter addresses these
142 > issues by looking at the structural and transport behavior of SSD over
143 > a variety of temperatures with the purpose of utilizing the RF
144 > correction technique.  We then suggest modifications to the parameters
145 > that result in more realistic bulk phase behavior.  It should be noted
146 > that in a recent publication, some of the original investigators of
147 > the SSD water model have suggested adjustments to address abnormal
148 > density behavior (also observed here), calling the corrected model
149 > SSD1.\cite{Tan03} In the later sections of this chapter, we compare
150 > our modified variants of SSD with both the original SSD and SSD1
151 > models and discuss how our changes improve the depiction of water.
152  
153 < \section{Simulation Methods}
153 > \section{Simulation Methods}{\label{sec:waterSimMethods}
154  
155   Most of the calculations in this particular study were performed using
156 < a internally developed simulation code that was one of the precursors
156 > an internally developed simulation code that was one of the precursors
157   of the {\sc oopse} molecular dynamics (MD) package.\cite{Meineke05}
158 < All of the capabilities of this code have been efficiently
159 < incorporated into {\sc oopse}, and calculation results are consistent
160 < between the two simulation packages. The later calculations involving
161 < the damped shifted force ({\sc sf}) techniques were performed using
160 < {\sc oopse}.
158 > All of the capabilities of this code have been incorporated into {\sc
159 > oopse}, and calculation results are consistent between the two
160 > simulation packages. The later calculations involving the damped
161 > shifted force ({\sc sf}) techniques were performed using {\sc oopse}.
162  
163   In the primary simulations of this study, long-range dipole-dipole
164 < interaction corrections were accounted for by using either the
165 < reaction field technique or a simple cubic switching function at the
166 < cutoff radius. Interestingly, one of the early applications of a
167 < reaction field was in Monte Carlo simulations of liquid
168 < water.\cite{Barker73} In this method, the magnitude of the reaction
169 < field acting on dipole $i$ is
164 > corrections were accounted for by using either the reaction field
165 > technique or a simple cubic switching function at the cutoff
166 > radius. Interestingly, one of the early applications of a reaction
167 > field was in Monte Carlo simulations of liquid water.\cite{Barker73}
168 > In this method, the magnitude of the reaction field acting on dipole
169 > $i$ is
170   \begin{equation}
171   \mathcal{E}_{i} = \frac{2(\varepsilon_{s} - 1)}{2\varepsilon_{s} + 1}
172   \frac{1}{r_{c}^{3}} \sum_{j\in{\mathcal{R}}} {\bf \mu}_{j} s(r_{ij}),
# Line 175 | Line 176 | and $s(r_{ij})$ is a cubic switching function.\cite{Al
176   ($r_{c}$), $\varepsilon_{s}$ is the dielectric constant imposed on the
177   system, ${\bf\mu}_{j}$ is the dipole moment vector of particle $j$,
178   and $s(r_{ij})$ is a cubic switching function.\cite{Allen87} The
179 < reaction field contribution to the total energy by particle $i$ is
179 > reaction field contribution to the total energy from particle $i$ is
180   given by $-\frac{1}{2}{\bf\mu}_{i}\cdot\mathcal{E}_{i}$ and the torque
181   on dipole $i$ by ${\bf\mu}_{i}\times\mathcal{E}_{i}$.\cite{Allen87} An
182   applied reaction field will alter the bulk orientational properties of
183   simulated water, and there is particular sensitivity of these
184   properties on changes in the length of the cutoff
185 < radius.\cite{vanderSpoel98} This variable behavior makes reaction
186 < field a less attractive method than the Ewald sum; however, for very
187 < large systems, the computational benefit of reaction field is is
187 < significant.
185 > radius.\cite{vanderSpoel98} This behavior makes the reaction field a
186 > less attractive method than the Ewald sum; however, for very large
187 > systems, the computational benefit of a reaction field is significant.
188  
189   In contrast to the simulations with a reaction field, we have also
190   performed a companion set of simulations {\it without} a surrounding
# Line 193 | Line 193 | To determine the preferred densities of the models, we
193   which can be used either with or without an active reaction field.
194  
195   To determine the preferred densities of the models, we performed
196 < simulations in the isobaric-isothermal ({\it NPT}) ensemble. All
196 > simulations in the isobaric-isothermal ($NPT$) ensemble. All
197   dynamical properties for these models were then obtained from
198 < microcanonical ({\it NVE}) simulations done at densities matching the
199 < {\it NPT} density for a particular target temperature.  The constant
198 > microcanonical ($NVE$) simulations done at densities matching the
199 > $NPT$ density for a particular target temperature.  The constant
200   pressure simulations were implemented using an integral thermostat and
201   barostat as outlined by Hoover.\cite{Hoover85,Hoover86} All molecules
202   were treated as non-linear rigid bodies. Vibrational constraints are
# Line 209 | Line 209 | method.\cite{Dullweber97,Evans77,Evans77b,Allen87} The
209   to carry out the integration of the equations of motion in place of
210   the more prevalent quaternion
211   method.\cite{Dullweber97,Evans77,Evans77b,Allen87} The reason behind
212 < this decision was that, in {\it NVE} simulations, the energy drift
212 > this decision was that, in $NVE$ simulations, the energy drift
213   when using quaternions was substantially greater than when using the
214 < {\sc dlm} method (Fig. \ref{fig:timeStepIntegration}). This steady
214 > {\sc dlm} method (figure \ref{fig:timeStepIntegration}). This steady
215   drift in the total energy has also been observed in other
216   studies.\cite{Kol97}
217  
# Line 233 | Line 233 | use time steps that are more that twice as long and st
233   simulations using the quaternion method and an identical time
234   step. This additional expense is justified because of the ability to
235   use time steps that are more that twice as long and still achieve the
236 < same energy conservation.
236 > same degree of energy conservation.
237  
238   Figure \ref{fig:timeStepIntegration} shows the resulting energy drift
239   at various time steps for both {\sc dlm} and quaternion
# Line 252 | Line 252 | of 4~fs. As expected, this drift increases dramaticall
252   unnoticeable for time steps up to 3~fs. We observed a slight energy
253   drift on the order of 0.012~kcal/mol per nanosecond with a time step
254   of 4~fs. As expected, this drift increases dramatically with increasing
255 < time step. To insure accuracy in our {\it NVE} simulations, time steps
256 < were set at 2~fs and were also kept at this value for {\it NPT}
255 > time step. To insure accuracy in our $NVE$ simulations, time steps
256 > were set at 2~fs and were also kept at this value for $NPT$
257   simulations.
258  
259   Proton-disordered ice crystals in both the I$_\textrm{h}$ and
# Line 261 | Line 261 | arranging the centers of mass of the SSD particles int
261   simulations. The I$_\textrm{h}$ crystals were formed by first
262   arranging the centers of mass of the SSD particles into a
263   ``hexagonal'' ice lattice of 1024 particles. Because of the crystal
264 < structure of I$_\textrm{h}$ ice, the simulation boxes were
264 > structure of ice I$_\textrm{h}$, the simulation boxes were
265   orthorhombic in shape with an edge length ratio of approximately
266   1.00$\times$1.06$\times$1.23. We then allowed the particles to orient
267   freely about their fixed lattice positions with angular momenta values
268   randomly sampled at 400~K. The rotational temperature was then scaled
269   down in stages to slowly cool the crystals to 25~K. The particles were
270   then allowed to translate with fixed orientations at a constant
271 < pressure of 1atm for 50~ps at 25~K. Finally, all constraints were
271 > pressure of 1~atm for 50~ps at 25~K. Finally, all constraints were
272   removed and the ice crystals were allowed to equilibrate for 50~ps at
273 < 25~K and a constant pressure of 1atm.  This procedure resulted in
274 < structurally stable I$_\textrm{h}$ ice crystals that obey the
273 > 25~K and a constant pressure of 1~atm.  This procedure resulted in
274 > structurally stable ice I$_\textrm{h}$ crystals that obey the
275   Bernal-Fowler rules.\cite{Bernal33,Rahman72} This method was also
276 < utilized in the making of diamond lattice I$_\textrm{c}$ ice crystals,
276 > utilized in the making of diamond lattice ice I$_\textrm{c}$ crystals,
277   with each cubic simulation box consisting of either 512 or 1000
278   particles. Only isotropic volume fluctuations were performed under
279   constant pressure, so the ratio of edge lengths remained constant
# Line 282 | Line 282 | Melting studies were performed on the randomized ice c
282   \section{SSD Density Behavior}
283  
284   Melting studies were performed on the randomized ice crystals using
285 < the {\it NPT} ensemble. During melting simulations, the melting
285 > the $NPT$ ensemble. During melting simulations, the melting
286   transition and the density maximum can both be observed, provided that
287   the density maximum occurs in the liquid and not the supercooled
288   regime. It should be noted that the calculated melting temperature
# Line 329 | Line 329 | densities of several other models and experiment obtai
329   The density maximum for SSD compares quite favorably to other simple
330   water models. Figure \ref{fig:ssdDense} also shows calculated
331   densities of several other models and experiment obtained from other
332 < sources.\cite{Jorgensen98b,Baez94,CRC80} Of the listed simple water
333 < models, SSD has a temperature closest to the experimentally observed
334 < density maximum. Of the {\it charge-based} models in figure
332 > sources.\cite{Jorgensen98b,Baez94,CRC80} Of the water models shown,
333 > SSD has a temperature closest to the experimentally observed density
334 > maximum. Of the {\it charge-based} models in figure
335   \ref{fig:ssdDense}, TIP4P has a density maximum behavior most like
336   that seen in SSD. Though not included in this plot, it is useful to
337   note that TIP5P has a density maximum nearly identical to the
338   experimentally measured temperature (see section
339 < \ref{sec:t5peDensity}.
339 > \ref{sec:t5peDensity}).
340  
341   Liquid state densities in water have been observed to be dependent on
342   the cutoff radius ($R_\textrm{c}$), both with and without the use of a
# Line 345 | Line 345 | within error and showed no significant trend toward lo
345   radius of 12~\AA\, complementing the 9~\AA\ $R_\textrm{c}$ used in the
346   previous SSD simulations. All of the resulting densities overlapped
347   within error and showed no significant trend toward lower or higher
348 < densities in simulations both with and without reaction field.
348 > densities in simulations both with and without a reaction field.
349  
350   The key feature to recognize in figure \ref{fig:ssdDense} is the
351   density scaling of SSD relative to other common models at any given
352   temperature. SSD assumes a lower density than any of the other listed
353 < models at the same pressure, behavior which is especially apparent at
353 > models at the same pressure. This behavior is especially apparent at
354   temperatures greater than 300~K. Lower than expected densities have
355   been observed for other systems using a reaction field for long-range
356   electrostatic interactions, so the most likely reason for the reduced
# Line 364 | Line 364 | decreasing densities at higher temperatures; however,
364   water. The shape of the curve is similar to the curve produced from
365   SSD simulations using reaction field, specifically the rapidly
366   decreasing densities at higher temperatures; however, a shift in the
367 < density maximum location, down to 245~K, is observed. This is a more
368 < accurate comparison to the other listed water models, in that no long
369 < range corrections were applied in those
367 > density maximum location, down to 245~K, is also observed. This is a
368 > more accurate comparison to the other listed water models, in that no
369 > long-range corrections were applied in those
370   simulations.\cite{Baez94,Jorgensen98b} However, even without the
371   reaction field, the density around 300~K is still significantly lower
372   than experiment and comparable water models. This anomalous behavior
373   was what lead Tan {\it et al.} to recently reparametrize
374 < SSD.\cite{Tan03} Throughout the remainder of the paper our
374 > SSD.\cite{Tan03} Throughout the remainder of this chapter our
375   reparametrizations of SSD will be compared with their newer SSD1
376   model.
377  
# Line 379 | Line 379 | important when using the model to study permeation or
379  
380   Accurate dynamical properties of a water model are particularly
381   important when using the model to study permeation or transport across
382 < biological membranes.  In order to probe transport in bulk water, {\it
383 < NVE} simulations were performed at the average densities obtained from
384 < the {\it NPT} simulations at an identical target
385 < temperature. Simulations started with randomized velocities and
386 < underwent 50~ps of temperature scaling and 50~ps of constant energy
387 < equilibration before a 200~ps data collection run. Diffusion constants
388 < were calculated via linear fits to the long-time behavior of the
389 < mean-square displacement as a function of time.\cite{Allen87} The
390 < averaged results from five sets of {\it NVE} simulations are displayed
391 < in figure \ref{fig:ssdDiffuse}, alongside experimental, SPC/E, and TIP5P
392 < results.\cite{Gillen72,Holz00,Baez94,Mahoney01}
382 > biological membranes.  In order to probe transport in bulk water,
383 > $NVE$ simulations were performed at the average densities obtained
384 > from the $NPT$ simulations at the target temperature. Simulations
385 > started with randomized velocities and underwent 50~ps of temperature
386 > scaling and 50~ps of constant energy equilibration before a 200~ps
387 > data collection run. Diffusion constants were calculated via linear
388 > fits to the long-time behavior of the mean-square displacement as a
389 > function of time.\cite{Allen87} The average results from the five sets
390 > of $NVE$ simulations are displayed in figure \ref{fig:ssdDiffuse},
391 > alongside experimental, SPC/E, and TIP5P
392 > results.\cite{Gillen72,Holz00,Baez94,Mahoney01}
393  
394   \begin{figure}
395   \centering
# Line 414 | Line 414 | since the densities of both TIP5P and SSD are lower th
414   temperatures and tend to diffuse too rapidly at higher temperatures.
415   This behavior at higher temperatures is not particularly surprising
416   since the densities of both TIP5P and SSD are lower than experimental
417 < water densities at higher temperatures.  When calculating the
417 > water densities at these higher temperatures.  When calculating the
418   diffusion coefficients for SSD at experimental densities (instead of
419 < the densities from the {\it NPT} simulations), the resulting values
420 < fall more in line with experiment at these temperatures.
419 > the densities from the $NPT$ simulations), the resulting values
420 > fall more in line with experiment at the respective temperatures.
421  
422   \section{Structural Changes and Characterization}
423  
424   By starting the simulations from the crystalline state, we can get an
425   estimation of the $T_\textrm{m}$ of the ice structure, and beyond the
426 < melting point, we study the phase behavior of the liquid. The constant
427 < pressure heat capacity ($C_\textrm{p}$) was monitored to locate
428 < $T_\textrm{m}$ in each of the simulations. In the melting simulations
429 < of the 1024 particle ice I$_\textrm{h}$ simulations, a large spike in
430 < $C_\textrm{p}$ occurs at 245~K, indicating a first order phase
431 < transition for the melting of these ice crystals (see figure
432 < \ref{fig:ssdCp}. When the reaction field is turned off, the melting
426 > melting point, we can study the phase behavior of the liquid. The
427 > constant pressure heat capacity ($C_\textrm{p}$) was monitored to
428 > locate $T_\textrm{m}$ in each of the simulations. In the melting
429 > simulations of the 1024 particle ice I$_\textrm{h}$ crystals, a large
430 > spike in $C_\textrm{p}$ occurs at 245~K, indicating a first order
431 > phase transition for the melting of these ice crystals (see figure
432 > \ref{fig:ssdCp}). When the reaction field is turned off, the melting
433   transition occurs at 235~K.  These melting transitions are considerably
434   lower than the experimental value of 273~K, indicating that the solid
435   ice I$_\textrm{h}$ is not thermodynamically preferred relative to the
# Line 459 | Line 459 | correlation values below 0.5 and black areas have valu
459   \centering
460   \includegraphics[width=2.5in]{./figures/corrDiag.pdf}
461   \caption{ An illustration of angles involved in the correlations
462 < observed in figure \ref{fig:contour}.}
462 > displayed in figure \ref{fig:contour}.}
463   \label{fig:corrAngle}
464   \end{figure}
465  
466   Additional analysis of the melting process was performed using
467   two-dimensional structure and dipole angle correlations. Expressions
468   for these correlations are as follows:
469 –
469   \begin{equation}
470 < g_{\text{AB}}(r,\cos\theta) = \frac{V}{N_\text{A}N_\text{B}}\langle\sum_{i\in\text{A}}\sum_{j\in\text{B}}\delta(\cos\theta-\cos\theta_{ij})\delta(r-\left|{\bf r}_{ij}\right|)\rangle\ ,
470 > g_{\text{AB}}(r,\cos\theta) = \frac{V}{N_\text{A}N_\text{B}}\left\langle\sum_{i\in\text{A}}\sum_{j\in\text{B}}\delta(\cos\theta-\cos\theta_{ij})\delta(r-\left|{\bf r}_{ij}\right|)\right\rangle\ ,
471   \end{equation}
472   \begin{equation}
473   g_{\text{AB}}(r,\cos\omega) =
474 < \frac{V}{N_\text{A}N_\text{B}}\langle\sum_{i\in\text{A}}\sum_{j\in\text{B}}\delta(\cos\omega-\cos\omega_{ij})\delta(r-\left|{\bf r}_{ij}\right|)\rangle\ ,
474 > \frac{V}{N_\text{A}N_\text{B}}\left\langle\sum_{i\in\text{A}}\sum_{j\in\text{B}}\delta(\cos\omega-\cos\omega_{ij})\delta(r-\left|{\bf r}_{ij}\right|)\right\rangle\ ,
475   \end{equation}
476   where $\theta$ and $\omega$ refer to the angles shown in figure
477   \ref{fig:corrAngle}. By binning over both distance and the cosine of the
# Line 492 | Line 491 | include a fraction of aligned dipole peaks within the
491   rapidly. The first solvation shell still shows the strong effect of
492   the sticky-potential, although it covers a larger area, extending to
493   include a fraction of aligned dipole peaks within the first solvation
494 < shell. The latter peaks lose due to thermal motion and as the
494 > shell. The latter peaks are lost due to thermal motion and as the
495   competing dipole force overcomes the sticky potential's tight
496   tetrahedral structuring of the crystal.
497  
# Line 505 | Line 504 | tightly bunched group of axially arranged dipoles that
504   this split character alters to show the leading 4~\AA\ peak dominated
505   by equatorial anti-parallel dipole orientations. There is also a
506   tightly bunched group of axially arranged dipoles that most likely
507 < consist of the smaller fraction of aligned dipole pairs. The trailing
507 > consist of a smaller fraction of aligned dipole pairs. The trailing
508   component of the split peak at 5~\AA\ is dominated by aligned dipoles
509   that assume hydrogen bond arrangements similar to those seen in the
510   first solvation shell. This evidence indicates that the dipole pair
# Line 579 | Line 578 | $r_u^\prime$ & (\AA) & 4.00 & 4.00 & 3.35 & 3.35\\
578   \begin{figure}
579   \centering
580   \includegraphics[width=4.5in]{./figures/newGofRCompare.pdf}
581 < \caption{ Plots showing the experimental $g(r)$ (Ref. \cite{Hura00})
581 > \caption{ Plots showing the experimental $g(r)$ (reference \cite{Hura00})
582   with SSD/E and SSD1 without reaction field (top), as well as SSD/RF
583   and SSD1 with reaction field turned on (bottom). The changes in
584   parameters have lowered and broadened the first peak of SSD/E and
# Line 636 | Line 635 | improves the densities, these changes alone are insuff
635  
636   While adjusting the location and shape of the first peak of $g(r)$
637   improves the densities, these changes alone are insufficient to bring
638 < the system densities up to the values observed experimentally.  To
639 < further increase the densities, the dipole moments were increased in
640 < both of our adjusted models. Since SSD is a dipole based model, the
641 < structure and transport are very sensitive to changes in the dipole
642 < moment. The original SSD simply used the dipole moment calculated from
643 < the TIP3P water model, which at 2.35~D is significantly greater than
644 < the experimental gas phase value of 1.84~D. The larger dipole moment
645 < is a more realistic value and improves the dielectric properties of
646 < the fluid. Both theoretical and experimental measurements indicate a
647 < liquid phase dipole moment ranging from 2.4~D to values as high as
648 < 3.11~D, providing a substantial range of reasonable values for a
649 < dipole moment.\cite{Sprik91,Gubskaya02,Badyal00,Barriol64} Moderately
638 > the system densities up to the experimental values.  To further
639 > increase the densities, the dipole moments were increased in both of
640 > our adjusted models. Since SSD is a dipole based model, the structure
641 > and transport are very sensitive to changes in the dipole moment. The
642 > original SSD used the dipole moment of the TIP3P water model, which at
643 > 2.35~D is significantly greater than the experimental gas phase value
644 > of 1.84~D. The larger dipole moment is a more realistic value and
645 > improves the dielectric properties of the fluid. Both theoretical and
646 > experimental measurements indicate a liquid phase dipole moment
647 > ranging from 2.4~D to values as high as 3.11~D, providing a
648 > substantial range of reasonable values for a dipole
649 > moment.\cite{Sprik91,Gubskaya02,Badyal00,Barriol64} Moderately
650   increasing the dipole moments to 2.42 and 2.48~D for SSD/E and SSD/RF,
651   respectively, leads to significant changes in the density and
652   transport of the water models.
# Line 655 | Line 654 | In order to demonstrate the benefits of these reparame
654   \subsection{Density Behavior}
655  
656   In order to demonstrate the benefits of these reparametrizations, we
657 < performed a series of {\it NPT} and {\it NVE} simulations to probe the
658 < density and transport properties of the adapted models and compare the
657 > performed a series of $NPT$ and $NVE$ simulations to probe the density
658 > and transport properties of the adapted models and compared the
659   results to the original SSD model. This comparison involved full {\it
660 < NPT} melting sequences for both SSD/E and SSD/RF, as well as {\it NVE}
660 > NPT} melting sequences for both SSD/E and SSD/RF, as well as $NVE$
661   transport calculations at the calculated self-consistent
662   densities. Again, the results were obtained from five separate
663   simulations of 1024 particle systems, and the melting sequences were
664   started from different ice I$_\textrm{h}$ crystals constructed as
665 < described previously. Each {\it NPT} simulation was equilibrated for
665 > described previously. Each $NPT$ simulation was equilibrated for
666   100~ps before a 200~ps data collection run at each temperature step,
667   and the final configuration from the previous temperature simulation
668 < was used as a starting point. All {\it NVE} simulations had the same
668 > was used as a starting point. All $NVE$ simulations had the same
669   thermalization, equilibration, and data collection times as stated
670   previously.
671  
# Line 714 | Line 713 | simulations.}
713   \label{fig:ssdrfDense}
714   \end{figure}
715  
716 < Including the reaction field long-range correction results in a more
716 > Including the reaction field long-range correction results is a more
717   interesting comparison.  A density profile including SSD/RF and SSD1
718   with an active reaction field is shown in figure \ref{fig:ssdrfDense}.
719   As observed in the simulations without a reaction field, the densities
# Line 739 | Line 738 | SSD1 (both without a reaction field) along with experi
738   \includegraphics[width=\linewidth]{./figures/ssdeDiffuse.pdf}
739   \caption{ The diffusion constants calculated from SSD/E and
740   SSD1 (both without a reaction field) along with experimental
741 < results.\cite{Gillen72,Holz00} The {\it NVE} calculations were
742 < performed at the average densities from the {\it NPT} simulations for
741 > results.\cite{Gillen72,Holz00} The $NVE$ calculations were
742 > performed at the average densities from the $NPT$ simulations for
743   the respective models. SSD/E is slightly more mobile than experiment
744   at all of the temperatures, but it is closer to experiment at
745   biologically relevant temperatures than SSD1 without a long-range
# Line 755 | Line 754 | to SSD1 without an active reaction field at the densit
754   of SSD be maintained or improved. Figure \ref{fig:ssdeDiffuse}
755   compares the temperature dependence of the diffusion constant of SSD/E
756   to SSD1 without an active reaction field at the densities calculated
757 < from their respective {\it NPT} simulations at 1 atm. The diffusion
757 > from their respective $NPT$ simulations at 1 atm. The diffusion
758   constant for SSD/E is consistently higher than experiment, while SSD1
759   remains lower than experiment until relatively high temperatures
760   (around 360~K). Both models follow the shape of the experimental curve
# Line 778 | Line 777 | SSD1 (both with an active reaction field) along with e
777   \includegraphics[width=\linewidth]{./figures/ssdrfDiffuse.pdf}
778   \caption{ The diffusion constants calculated from SSD/RF and
779   SSD1 (both with an active reaction field) along with experimental
780 < results.\cite{Gillen72,Holz00} The {\it NVE} calculations were
781 < performed at the average densities from the {\it NPT} simulations for
780 > results.\cite{Gillen72,Holz00} The $NVE$ calculations were
781 > performed at the average densities from the $NPT$ simulations for
782   both of the models. SSD/RF captures the self-diffusion of water
783   throughout most of this temperature range. The increasing diffusion
784   constants at high temperatures for both models can be attributed to
# Line 836 | Line 835 | where $\rho$ is the number density of the specified pa
835   n_H = 4\pi\rho_{\text{OH}}\int_{0}^{b}r^2g_{\textrm{OH}}(r)dr,
836   \end{equation}
837   where $\rho$ is the number density of the specified pair interactions,
838 < $a$ and $b$ are the radial locations of the minima following the first
839 < peak in $g_\textrm{OO}(r)$ or $g_\textrm{OH}(r)$ respectively. The
840 < number of hydrogen bonds stays relatively constant across all of the
841 < models, but the coordination numbers of SSD/E and SSD/RF show an
842 < improvement over SSD1.  This improvement is primarily due to extension
843 < of the first solvation shell in the new parameter sets.  Because $n_H$
844 < and $n_C$ are nearly identical in SSD1, it appears that all molecules
845 < in the first solvation shell are involved in hydrogen bonds.  Since
846 < $n_H$ and $n_C$ differ in the newly parameterized models, the
847 < orientations in the first solvation shell are a bit more ``fluid''.
848 < Therefore SSD1 over-structures the first solvation shell and our
849 < reparametrizations have returned this shell to more realistic
850 < liquid-like behavior.
838 > and $a$ and $b$ are the radial locations of the minima following the
839 > first peak in $g_\textrm{OO}(r)$ or $g_\textrm{OH}(r)$
840 > respectively. The number of hydrogen bonds stays relatively constant
841 > across all of the models, but the coordination numbers of SSD/E and
842 > SSD/RF show an improvement over SSD1.  This improvement is primarily
843 > due to extension of the first solvation shell in the new parameter
844 > sets.  Because $n_H$ and $n_C$ are nearly identical in SSD1, it
845 > appears that all molecules in the first solvation shell are involved
846 > in hydrogen bonds.  Since $n_H$ and $n_C$ differ in the newly
847 > parameterized models, the orientations in the first solvation shell
848 > are a bit more ``fluid''.  Therefore SSD1 over-structures the first
849 > solvation shell and our reparametrizations have returned this shell to
850 > more realistic liquid-like behavior.
851  
852   The time constants for the orientational autocorrelation functions
853   are also displayed in Table \ref{tab:liquidProperties}. The dipolar
# Line 865 | Line 864 | regime ($t > \tau_l$).\cite{Rothschild84} Calculation
864   correlation functions, the orientational relaxation time of the dipole
865   vector can be calculated from an exponential fit in the long-time
866   regime ($t > \tau_l$).\cite{Rothschild84} Calculation of these time
867 < constants were averaged over five detailed {\it NVE} simulations
867 > constants were averaged over five detailed $NVE$ simulations
868   performed at the ambient conditions for each of the respective
869   models. It should be noted that the commonly cited value of 1.9~ps for
870 < $\tau_2$ was determined from the NMR data in Ref. \cite{Krynicki66} at
870 > $\tau_2$ was determined from the NMR data in reference \cite{Krynicki66} at
871   a temperature near 34$^\circ$C.\cite{Rahman71} Because of the strong
872   temperature dependence of $\tau_2$, it is necessary to recalculate it
873   at 298~K to make proper comparisons. The value shown in Table
874   \ref{tab:liquidProperties} was calculated from the same NMR data in the
875 < fashion described in Ref. \cite{Krynicki66}. Similarly, $\tau_1$ was
876 < recomputed for 298~K from the data in Ref. \cite{Eisenberg69}.
875 > fashion described in reference \cite{Krynicki66}. Similarly, $\tau_1$ was
876 > recomputed for 298~K from the data in reference \cite{Eisenberg69}.
877   Again, SSD/E and SSD/RF show improved behavior over SSD1, both with
878   and without an active reaction field. Turning on the reaction field
879   leads to much improved time constants for SSD1; however, these results
# Line 886 | Line 885 | In section \ref{sec:dampingMultipoles}, a method was d
885   \subsection{SSD/RF and Damped Electrostatics}\label{sec:ssdrfDamped}
886  
887   In section \ref{sec:dampingMultipoles}, a method was described for
888 < applying the damped {\sc sf} or {\sc sp} techniques to for systems
888 > applying the damped {\sc sf} or {\sc sp} techniques to systems
889   containing point multipoles. The SSD family of water models is the
890   perfect test case because of the dipole-dipole (and
891   charge-dipole/quadrupole) interactions that are present. The {\sc sf}
# Line 925 | Line 924 | about because we neglect the pressure effect due to th
924   surprisingly well. The average density shows a modest increase when
925   using damped electrostatics in place of the reaction field. This comes
926   about because we neglect the pressure effect due to the surroundings
927 < outside of the cuttoff, instead relying on screening effects to
927 > outside of the cutoff, instead relying on screening effects to
928   neutralize electrostatic interactions at long distances. The $C_p$
929   also shows a slight increase, indicating greater fluctuation in the
930   enthalpy at constant pressure. The only other differences between the
931   damped and reaction field results are the dipole reorientational time
932   constants, $\tau_1$ and $\tau_2$. When using damped electrostatics,
933 < the water molecules relax more quickly and are almost identical to the
934 < experimental values. These results indicate that not only is it
935 < reasonable to use damped electrostatics with SSD/RF, it is recommended
936 < if capturing realistic dynamics is of primary importance. This is an
937 < encouraging result because of the more varied applicability of damping
938 < over the reaction field technique. Rather than be limited to
939 < homogeneous systems, SSD/RF can be used effectively with mixed
940 < systems, such as dissolved ions, dissolved organic molecules, or even
942 < proteins.
933 > the water molecules relax more quickly and exhibit behavior very
934 > similar to the experimental values. These results indicate that not
935 > only is it reasonable to use damped electrostatics with SSD/RF, it is
936 > recommended if capturing realistic dynamics is of primary
937 > importance. This is an encouraging result because the damping methods
938 > are more generally applicable than reaction field. Using damping,
939 > SSD/RF can be used effectively with mixed systems, such as dissolved
940 > ions, dissolved organic molecules, or even proteins.
941  
942   In addition to the properties tabulated in table
943   \ref{tab:dampedSSDRF}, we calculated the static dielectric constant
# Line 1025 | Line 1023 | peak of $g_\textrm{OO}(r)$ (see figure \ref{fig:tredGo
1023   peaks in $g_\textrm{OO}(r)$ and $g_\textrm{OH}(r)$, while the $\sigma$
1024   and $\epsilon$ values were varied to adjust the location of the first
1025   peak of $g_\textrm{OO}(r)$ (see figure \ref{fig:tredGofR}) and the
1026 < density. The $\sigma$ and $epsilon$ optimization was carried out by
1026 > density. The $\sigma$ and $\epsilon$ optimization was carried out by
1027   separating the Lennard-Jones potential into near entirely repulsive
1028   ($A$) and attractive ($C$) parts, such that
1029   \begin{equation}
# Line 1040 | Line 1038 | translational diffusion constant at 298~K and 1~atm.
1038   were made with the goal of capturing the experimental density and
1039   translational diffusion constant at 298~K and 1~atm.
1040  
1041 < Being that TRED is a two-point water model, we expect its
1042 < computational efficiency to fall some place in between the single and
1043 < three-point water models. In comparative simulations, TRED was
1044 < approximately 85\% slower than SSD/RF, while SPC/E was 225\% slower
1045 < than SSD/RF. While TRED loses some of the performance advantage of
1046 < SSD, it is still nearly twice as computationally efficient as SPC/E
1049 < and TIP3P.
1041 > Since TRED is a two-point water model, we expect its computational
1042 > efficiency to fall some place in between the one-point and three-point
1043 > water models. In comparative simulations, TRED was approximately 85\%
1044 > slower than SSD/RF, while SPC/E was 225\% slower than SSD/RF. While
1045 > TRED loses some of the performance advantage of SSD, it is still
1046 > nearly twice as computationally efficient as SPC/E and TIP3P.
1047  
1048   \begin{table}
1049   \caption{PROPERTIES OF TRED COMPARED WITH SSD/RF AND EXPERIMENT}
# Line 1107 | Line 1104 | electrostatic damping act in an equivalent fashion. Bo
1104   \ref{fig:dielectricMap}D) highlights the similarities in how these
1105   models respond to dielectric damping and how the dipolar and monopolar
1106   electrostatic damping act in an equivalent fashion. Both these
1107 < dielectric maps span a larger range than the 3, 4, and 5 point-charge
1108 < water models; however, the SSD/RF range is greater than TRED,
1109 < indicating that multipoles are a little more sensitive to damping than
1110 < monopoles.
1107 > dielectric maps span a larger range than the three-, four-, and
1108 > five-point water models; however, the SSD/RF range is greater than
1109 > TRED, indicating that multipoles are more sensitive to damping than
1110 > monopoles.
1111  
1112   The final dielectric comparison comes through the Debye relaxation
1113   time ($\tau_D$) or the collective dipolar relaxation time when
1114   assuming a Debye treatment for the dielectric
1115   relaxation.\cite{Chandra99,Kindt96} This value is calculated through
1116   equation (\ref{eq:reorientCorr}) applied to the total system dipole
1117 < moment. The values for both of the models are a slower than the
1117 > moment. The values for both of the models are slower than the
1118   experimental relaxation; however, they compare compare very well to
1119   experiment considering the Debye relaxation times calculated for the
1120   original SSD (11.95~ps) and the SPC/E (6.95~ps) and TIP3P (6.1~ps)
# Line 1129 | Line 1126 | of the self-diffusion constant were studied for the SS
1126  
1127   In the above sections, the density maximum and temperature dependence
1128   of the self-diffusion constant were studied for the SSD water model,
1129 < both with and without the use of reaction field, via a series of {\it
1130 < NPT} and {\it NVE} simulations. The constant pressure simulations
1131 < showed a density maximum near 260~K. In most cases, the calculated
1132 < densities were significantly lower than the densities obtained from
1133 < other water models (and experiment). Analysis of self-diffusion showed
1134 < SSD to capture the transport properties of water well in both the
1135 < liquid and supercooled liquid regimes.
1129 > both with and without the use of reaction field, via a series of $NPT$
1130 > and $NVE$ simulations. The constant pressure simulations showed a
1131 > density maximum near 260~K. In most cases, the calculated densities
1132 > were significantly lower than the densities obtained from other water
1133 > models (and experiment). Analysis of self-diffusion showed SSD to
1134 > capture the transport properties of water well in both the liquid and
1135 > supercooled liquid regimes.
1136  
1137   In order to correct the density behavior, we reparametrized the
1138   original SSD model for use both with and without a reaction field
# Line 1148 | Line 1145 | electrostatic conditions, validating the multipolar da
1145  
1146   We also showed that SSD/RF performs well under the alternative damped
1147   electrostatic conditions, validating the multipolar damping work in
1148 < the previous chapter. To improve the modeling of water when {\sc sf},
1149 < the TRED water model was developed. This model maintains improves upon
1150 < the thermodynamic properties of SSD/RF using damped electrostatics
1151 < while maintaining the exceptional depiction of water dynamics.
1148 > the previous chapter. To improve the modeling of water with the {\sc
1149 > sf} technique, the TRED water model was developed. This model improves
1150 > upon the thermodynamic properties of SSD/RF using damped
1151 > electrostatics while maintaining the exceptional depiction of water
1152 > dynamics.
1153  
1154   The simple water models investigated here are excellent choices for
1155 < representing explicit water in large scale simulations of biochemical
1155 > representing water in large scale simulations of biochemical
1156   systems. They are more computationally efficient than the common
1157 < charge based water models, and, in many cases, exhibit more realistic
1158 < bulk phase fluid properties. These models are one of the few cases in
1159 < which maximizing efficiency does not result in a loss in realistic
1160 < liquid water representation.
1157 > charge based water models, and often exhibit more realistic bulk phase
1158 > fluid properties. These models are one of the few cases in which
1159 > maximizing efficiency does not result in a loss in realistic
1160 > representation.

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines