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

Comparing trunk/ssdePaper/nptSSD.tex (file contents):
Revision 856 by chrisfen, Thu Nov 6 23:00:00 2003 UTC vs.
Revision 861 by chrisfen, Wed Nov 12 13:37:15 2003 UTC

# Line 1 | Line 1
1 < \documentclass[prb,aps,times,twocolumn,tabularx]{revtex4}
2 < %\documentclass[prb,aps,times,tabularx,preprint]{revtex4}
1 > %\documentclass[prb,aps,times,twocolumn,tabularx]{revtex4}
2 > \documentclass[prb,aps,times,tabularx,preprint]{revtex4}
3   \usepackage{amsmath}
4   \usepackage{graphicx}
5  
# Line 42 | Line 42 | experimental very well in both the normal and super-co
42   cases, the calculated densities were significantly lower than the
43   densities calculated in simulations of other water models. Analysis of
44   particle diffusion showed SSD to capture the transport properties of
45 < experimental very well in both the normal and super-cooled liquid
46 < regimes. In order to correct the density behavior, SSD was
45 > experimental water very well in both the normal and super-cooled
46 > liquid regimes. In order to correct the density behavior, SSD was
47   reparameterized for use both with and without a long-range interaction
48 < correction, SSD/RF and SSD/E respectively. In addition to correcting
49 < the abnormally low densities, these new versions were show to maintain
50 < or improve upon the transport and structural features of the original
51 < water model.
48 > correction, SSD/RF and SSD/E respectively. Compared to the density
49 > corrected version of SSD (SSD1), these modified models were shown to
50 > maintain or improve upon the structural and transport properties.
51   \end{abstract}
52  
53   \maketitle
# Line 145 | Line 144 | comparable models while maintaining the structural beh
144   calculations are simplified significantly. In the original Monte Carlo
145   simulations using this model, Ichiye \emph{et al.} reported a
146   calculation speed up of up to an order of magnitude over other
147 < comparable models while maintaining the structural behavior of
147 > comparable models, while maintaining the structural behavior of
148   water.\cite{Ichiye96} In the original molecular dynamics studies, it
149   was shown that SSD improves on the prediction of many of water's
150   dynamical properties over TIP3P and SPC/E.\cite{Ichiye99} This
151   attractive combination of speed and accurate depiction of solvent
152   properties makes SSD a model of interest for the simulation of large
153 < scale biological systems, such as membrane phase behavior, a specific
155 < interest within our group.
153 > scale biological systems, such as membrane phase behavior.
154  
155   One of the key limitations of this water model, however, is that it
156   has been parameterized for use with the Ewald Sum technique for the
# Line 166 | Line 164 | utilizing the RF correction technique. Towards the end
164   a cutoff that lacks any form of long range correction. This study
165   addresses these issues by looking at the structural and transport
166   behavior of SSD over a variety of temperatures, with the purpose of
167 < utilizing the RF correction technique. Towards the end, we suggest
167 > utilizing the RF correction technique. Toward the end, we suggest
168   alterations to the parameters that result in more water-like
169   behavior. It should be noted that in a recent publication, some the
170   original investigators of the SSD water model have put forth
171 < adjustments to the original SSD water model to address abnormal
172 < density behavior (also observed here), calling the corrected model
173 < SSD1.\cite{Ichiye03} This study will consider this new model's
174 < behavior as well, and hopefully improve upon its depiction of water
175 < under conditions without the Ewald Sum.
171 > adjustments to the SSD water model to address abnormal density
172 > behavior (also observed here), calling the corrected model
173 > SSD1.\cite{Ichiye03} This study will make comparisons with this new
174 > model's behavior with the goal of improving upon the depiction of
175 > water under conditions without the Ewald Sum.
176  
177   \section{Methods}
178  
# Line 204 | Line 202 | SSD more compatible with a reaction field.
202   is dramatic. To address some of the dynamical property alterations due
203   to the use of reaction field, simulations were also performed without
204   a surrounding dielectric and suggestions are proposed on how to make
205 < SSD more compatible with a reaction field.
205 > SSD more accurate both with and without a reaction field.
206  
207   Simulations were performed in both the isobaric-isothermal and
208   microcanonical ensembles. The constant pressure simulations were
# Line 310 | Line 308 | the heat capacity, in addition to determining the dens
308   Melting studies were performed on the randomized ice crystals using
309   constant pressure and temperature dynamics. By performing melting
310   simulations, the melting transition can be determined by monitoring
311 < the heat capacity, in addition to determining the density maximum,
311 > the heat capacity, in addition to determining the density maximum -
312   provided that the density maximum occurs in the liquid and not the
313 < supercooled regimes. An ensemble average from five separate melting
313 > supercooled regime. An ensemble average from five separate melting
314   simulations was acquired, each starting from different ice crystals
315   generated as described previously. All simulations were equilibrated
316   for 100 ps prior to a 200 ps data collection run at each temperature
317 < setting, ranging from 25 to 400 K, with a maximum degree increment of
318 < 25 K. For regions of interest along this stepwise progression, the
319 < temperature increment was decreased from 25 K to 10 and then 5 K. The
320 < above equilibration and production times were sufficient in that the
321 < system volume fluctuations dampened out in all but the very cold
322 < simulations (below 225 K).
317 > setting. The temperature range of study spanned from 25 to 400 K, with
318 > a maximum degree increment of 25 K. For regions of interest along this
319 > stepwise progression, the temperature increment was decreased from 25
320 > K to 10 and 5 K. The above equilibration and production times were
321 > sufficient in that the system volume fluctuations dampened out in all
322 > but the very cold simulations (below 225 K).
323  
324   \subsection{Density Behavior}
325 < In the initial average density versus temperature plot, the density
326 < maximum appears between 255 and 265 K. The calculated densities within
327 < this range were nearly indistinguishable, as can be seen in the zoom
328 < of this region of interest, shown in figure
329 < \ref{dense1}. The greater certainty of the average value at 260 K makes
330 < a good argument for the actual density maximum residing at this
331 < midpoint value. Figure \ref{dense1} was constructed using ice $I_h$
332 < crystals for the initial configuration; and though not pictured, the
333 < simulations starting from ice $I_c$ crystal configurations showed
334 < similar results, with a liquid-phase density maximum in this same
335 < region (between 255 and 260 K). In addition, the $I_c$ crystals are
336 < more fragile than the $I_h$ crystals, leading them to deform into a
337 < dense glassy state at lower temperatures. This resulted in an overall
338 < low temperature density maximum at 200 K, but they still retained a
339 < common liquid state density maximum with the $I_h$ simulations.
325 > Initial simulations focused on the original SSD water model, and an
326 > average density versus temperature plot is shown in figure
327 > \ref{dense1}. Note that the density maximum when using a reaction
328 > field appears between 255 and 265 K, where the calculated densities
329 > within this range were nearly indistinguishable. The greater certainty
330 > of the average value at 260 K makes a good argument for the actual
331 > density maximum residing at this midpoint value. Figure \ref{dense1}
332 > was constructed using ice $I_h$ crystals for the initial
333 > configuration; and though not pictured, the simulations starting from
334 > ice $I_c$ crystal configurations showed similar results, with a
335 > liquid-phase density maximum in this same region (between 255 and 260
336 > K). In addition, the $I_c$ crystals are more fragile than the $I_h$
337 > crystals, leading them to deform into a dense glassy state at lower
338 > temperatures. This resulted in an overall low temperature density
339 > maximum at 200 K, but they still retained a common liquid state
340 > density maximum with the $I_h$ simulations.
341  
342   \begin{figure}
343   \includegraphics[width=65mm,angle=-90]{dense2.eps}
# Line 348 | Line 347 | common liquid state density maximum with the $I_h$ sim
347   change in densities observed when turning off the reaction field. The
348   the lower than expected densities for the SSD model were what
349   prompted the original reparameterization.\cite{Ichiye03}}
350 < \label{dense2}
350 > \label{dense1}
351   \end{figure}
352  
353   The density maximum for SSD actually compares quite favorably to other
354 < simple water models. Figure \ref{dense2} shows a plot of these
355 < findings with the density progression of several other models and
357 < experiment obtained from other
354 > simple water models. Figure \ref{dense1} also shows calculated
355 > densities of several other models and experiment obtained from other
356   sources.\cite{Jorgensen98b,Clancy94,CRC80} Of the listed simple water
357   models, SSD has results closest to the experimentally observed water
358   density maximum. Of the listed water models, TIP4P has a density
359 < maximum behavior most like that seen in SSD. Though not shown, it is
360 < useful to note that TIP5P has a water density maximum nearly identical
361 < to experiment.
359 > maximum behavior most like that seen in SSD. Though not included in
360 > this particular plot, it is useful to note that TIP5P has a water
361 > density maximum nearly identical to experiment.
362  
365 – Possibly of more importance is the density scaling of SSD relative to
366 – other common models at any given temperature (Fig. \ref{dense2}). Note
367 – that the SSD model assumes a lower density than any of the other
368 – listed models at the same pressure, behavior which is especially
369 – apparent at temperatures greater than 300 K. Lower than expected
370 – densities have been observed for other systems with the use of a
371 – reaction field for long-range electrostatic interactions, so the most
372 – likely reason for these significantly lower densities in these
373 – simulations is the presence of the reaction field.\cite{Berendsen98}
374 – In order to test the effect of the reaction field on the density of
375 – the systems, the simulations were repeated for the temperature region
376 – of interest without a reaction field present. The results of these
377 – simulations are also displayed in figure \ref{dense2}. Without
378 – reaction field, these densities increase considerably to more
379 – experimentally reasonable values, especially around the freezing point
380 – of liquid water. The shape of the curve is similar to the curve
381 – produced from SSD simulations using reaction field, specifically the
382 – rapidly decreasing densities at higher temperatures; however, a slight
383 – shift in the density maximum location, down to 245 K, is
384 – observed. This is probably a more accurate comparison to the other
385 – listed water models in that no long range corrections were applied in
386 – those simulations.\cite{Clancy94,Jorgensen98b}
387 –
363   It has been observed that densities are dependent on the cutoff radius
364   used for a variety of water models in simulations both with and
365   without the use of reaction field.\cite{Berendsen98} In order to
# Line 396 | Line 371 | use of a longer cutoff radius results in a near doubli
371   radius, both for simulations with and without reaction field. These
372   results indicate that there is no major benefit in choosing a longer
373   cutoff radius in simulations using SSD. This is comforting in that the
374 < use of a longer cutoff radius results in a near doubling of the time
375 < required to compute a single trajectory.
374 > use of a longer cutoff radius results in significant increases in the
375 > time required to obtain a single trajectory.
376  
377 + The most important thing to recognize in figure \ref{dense1} is the
378 + density scaling of SSD relative to other common models at any given
379 + temperature. Note that the SSD model assumes a lower density than any
380 + of the other listed models at the same pressure, behavior which is
381 + especially apparent at temperatures greater than 300 K. Lower than
382 + expected densities have been observed for other systems with the use
383 + of a reaction field for long-range electrostatic interactions, so the
384 + most likely reason for these significantly lower densities in these
385 + simulations is the presence of the reaction
386 + field.\cite{Berendsen98,Nezbeda02} In order to test the effect of the
387 + reaction field on the density of the systems, the simulations were
388 + repeated without a reaction field present. The results of these
389 + simulations are also displayed in figure \ref{dense1}. Without
390 + reaction field, these densities increase considerably to more
391 + experimentally reasonable values, especially around the freezing point
392 + of liquid water. The shape of the curve is similar to the curve
393 + produced from SSD simulations using reaction field, specifically the
394 + rapidly decreasing densities at higher temperatures; however, a shift
395 + in the density maximum location, down to 245 K, is observed. This is
396 + probably a more accurate comparison to the other listed water models,
397 + in that no long range corrections were applied in those
398 + simulations.\cite{Clancy94,Jorgensen98b} However, even without a
399 + reaction field, the density around 300 K is still significantly lower
400 + than experiment and comparable water models. This anomalous behavior
401 + was what lead Ichiye \emph{et al.} to recently reparameterize SSD and
402 + make SSD1.\cite{Ichiye03} In discussing potential adjustments later in
403 + this paper, all comparisons were performed with this new model.
404 +
405   \subsection{Transport Behavior}
406   Of importance in these types of studies are the transport properties
407   of the particles and how they change when altering the environmental
# Line 434 | Line 437 | experiment at these temperatures, albeit not at standa
437   SSD are lower than experimental water at temperatures higher than room
438   temperature. When calculating the diffusion coefficients for SSD at
439   experimental densities, the resulting values fall more in line with
440 < experiment at these temperatures, albeit not at standard
438 < pressure. Results under these conditions can be found later in this
439 < paper.
440 > experiment at these temperatures, albeit not at standard pressure.
441  
442   \subsection{Structural Changes and Characterization}
443   By starting the simulations from the crystalline state, the melting
# Line 449 | Line 450 | surprising in that SSD is a simple rigid body model wi
450   melting of these ice crystals. When the reaction field is turned off,
451   the melting transition occurs at 235 K.  These melting transitions are
452   considerably lower than the experimental value, but this is not
453 < surprising in that SSD is a simple rigid body model with a fixed
453 < dipole.
453 > surprising when considering the simplicity of the SSD model.
454  
455   \begin{figure}
456   \includegraphics[width=85mm]{fullContours.eps}
# Line 462 | Line 462 | Additional analyses for understanding the melting phas
462   \label{contour}
463   \end{figure}
464  
465 – Additional analyses for understanding the melting phase-transition
466 – process were performed via two-dimensional structure and dipole angle
467 – correlations. Expressions for the correlations are as follows:
468 –
465   \begin{figure}
466   \includegraphics[width=45mm]{corrDiag.eps}
467   \caption{Two dimensional illustration of the angles involved in the
# Line 473 | Line 469 | correlations observed in figure \ref{contour}.}
469   \label{corrAngle}
470   \end{figure}
471  
472 + Additional analysis of the melting phase-transition process was
473 + performed by using two-dimensional structure and dipole angle
474 + correlations. Expressions for these correlations are as follows:
475 +
476   \begin{multline}
477   g_{\text{AB}}(r,\cos\theta) = \\
478   \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|\mathbf{r}_{ij}\right|)\rangle\ ,
# Line 481 | Line 481 | where $\theta$ and $\omega$ refer to the angles shown
481   g_{\text{AB}}(r,\cos\omega) = \\
482   \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|\mathbf{r}_{ij}\right|)\rangle\ ,
483   \end{multline}
484 < where $\theta$ and $\omega$ refer to the angles shown in the above
485 < illustration. By binning over both distance and the cosine of the
484 > where $\theta$ and $\omega$ refer to the angles shown in figure
485 > \ref{corrAngle}. By binning over both distance and the cosine of the
486   desired angle between the two dipoles, the g(\emph{r}) can be
487   dissected to determine the common dipole arrangements that constitute
488   the peaks and troughs. Frames A and B of figure \ref{contour} show a
489   relatively crystalline state of an ice $I_c$ simulation. The first
490 < peak of the g(\emph{r}) primarily consists of the preferred hydrogen
491 < bonding arrangements as dictated by the tetrahedral sticky potential,
490 > peak of the g(\emph{r}) consists primarily of the preferred hydrogen
491 > bonding arrangements as dictated by the tetrahedral sticky potential -
492   one peak for the donating and the other for the accepting hydrogen
493   bonds. Due to the high degree of crystallinity of the sample, the
494   second and third solvation shells show a repeated peak arrangement
495   which decays at distances around the fourth solvation shell, near the
496   imposed cutoff for the Lennard-Jones and dipole-dipole interactions.
497 < In the higher temperature simulation shown in frames C and D, the
498 < repeated peak features are significantly blurred. The first solvation
499 < shell still shows the strong effect of the sticky-potential, although
500 < it covers a larger area, extending to include a fraction of aligned
501 < dipole peaks within the first solvation shell. The latter peaks lose
502 < definition as thermal motion and the competing dipole force overcomes
503 < the sticky potential's tight tetrahedral structuring of the fluid.
497 > In the higher temperature simulation shown in frames C and D, these
498 > longer-ranged repeated peak features deteriorate rapidly. The first
499 > solvation shell still shows the strong effect of the sticky-potential,
500 > although it covers a larger area, extending to include a fraction of
501 > aligned dipole peaks within the first solvation shell. The latter
502 > peaks lose definition as thermal motion and the competing dipole force
503 > overcomes the sticky potential's tight tetrahedral structuring of the
504 > fluid.
505  
506   This complex interplay between dipole and sticky interactions was
507   remarked upon as a possible reason for the split second peak in the
# Line 517 | Line 518 | immediately outside of it's range (around 4 \AA), and
518   indicates that the dipole pair interaction begins to dominate outside
519   of the range of the dipolar repulsion term, with the primary
520   energetically favorable dipole arrangements populating the region
521 < immediately outside of it's range (around 4 \AA), and arrangements
522 < that seek to ideally satisfy both the sticky and dipole forces locate
523 < themselves just beyond this region (around 5 \AA).
521 > immediately outside this repulsion region (around 4 \AA), and
522 > arrangements that seek to ideally satisfy both the sticky and dipole
523 > forces locate themselves just beyond this initial buildup (around 5
524 > \AA).
525  
526   From these findings, the split second peak is primarily the product of
527 < the dipolar repulsion term of the sticky potential. In fact, the
528 < leading of the two peaks can be pushed out and merged with the outer
529 < split peak just by extending the switching function cutoff
530 < ($s^\prime(r_{ij})$) from its normal 4.0 \AA\ to values of 4.5 or even
531 < 5 \AA. This type of correction is not recommended for improving the
532 < liquid structure, because the second solvation shell will still be
533 < shifted too far out. In addition, this would have an even more
534 < detrimental effect on the system densities, leading to a liquid with a
535 < more open structure and a density considerably lower than the normal
536 < SSD behavior shown previously. A better correction would be to include
537 < the quadrupole-quadrupole interactions for the water particles outside
538 < of the first solvation shell, but this reduces the simplicity and
539 < speed advantage of SSD, so it is not the most desirable path to take.
527 > the dipolar repulsion term of the sticky potential. In fact, the inner
528 > peak can be pushed out and merged with the outer split peak just by
529 > extending the switching function cutoff ($s^\prime(r_{ij})$) from its
530 > normal 4.0 \AA\ to values of 4.5 or even 5 \AA. This type of
531 > correction is not recommended for improving the liquid structure,
532 > because the second solvation shell will still be shifted too far
533 > out. In addition, this would have an even more detrimental effect on
534 > the system densities, leading to a liquid with a more open structure
535 > and a density considerably lower than the normal SSD behavior shown
536 > previously. A better correction would be to include the
537 > quadrupole-quadrupole interactions for the water particles outside of
538 > the first solvation shell, but this reduces the simplicity and speed
539 > advantage of SSD.
540  
541 < \subsection{Adjusted Potentials: SSD/E and SSD/RF}
541 > \subsection{Adjusted Potentials: SSD/RF and SSD/E}
542   The propensity of SSD to adopt lower than expected densities under
543   varying conditions is troubling, especially at higher temperatures. In
544 < order to correct this behavior, it's necessary to adjust the force
545 < field parameters for the primary intermolecular interactions. In
546 < undergoing a reparameterization, it is important not to focus on just
547 < one property and neglect the other important properties. In this case,
548 < it would be ideal to correct the densities while maintaining the
549 < accurate transport properties.
544 > order to correct this model for use with a reaction field, it is
545 > necessary to adjust the force field parameters for the primary
546 > intermolecular interactions. In undergoing a reparameterization, it is
547 > important not to focus on just one property and neglect the other
548 > important properties. In this case, it would be ideal to correct the
549 > densities while maintaining the accurate transport properties.
550  
551   The possible parameters for tuning include the $\sigma$ and $\epsilon$
552   Lennard-Jones parameters, the dipole strength ($\mu$), and the sticky
# Line 570 | Line 572 | parameterized potentials - soft sticky dipole enhanced
572   potentials. The results of the reparameterizations are shown in table
573   \ref{params}. Note that both the tetrahedral attractive and dipolar
574   repulsive don't share the same lower cutoff ($r_l$) in the newly
575 < parameterized potentials - soft sticky dipole enhanced (SSD/E) and
576 < soft sticky dipole reaction field (SSD/RF).
575 > parameterized potentials - soft sticky dipole reaction field (SSD/RF -
576 > for use with a reaction field) and soft sticky dipole enhanced (SSD/E
577 > - an attempt to improve the liquid structure in simulations without a
578 > long-range correction).
579  
580   \begin{table}
581   \caption{Parameters for the original and adjusted models}
# Line 622 | Line 626 | g(\emph{r}) by comparing the original SSD (with and wi
626   adjustments to SSD were made while taking into consideration the new
627   experimental findings.\cite{Head-Gordon00_1} Figure \ref{grcompare}
628   shows the relocation of the first peak of the oxygen-oxygen
629 < g(\emph{r}) by comparing the original SSD (with and without reaction
630 < field), SSD-E, and SSD-RF to the new experimental results. Both the
631 < modified water models have shorter peaks that are brought in more
632 < closely to the experimental peak (as seen in the insets of figure
633 < \ref{grcompare}). This structural alteration was accomplished by a
634 < reduction in the Lennard-Jones $\sigma$ variable as well as adjustment
635 < of the sticky potential strength and cutoffs. The cutoffs for the
636 < tetrahedral attractive and dipolar repulsive terms were nearly swapped
637 < with each other. Isosurfaces of the original and modified sticky
638 < potentials are shown in figure \cite{isosurface}. In these
639 < isosurfaces, it is easy to see how altering the cutoffs changes the
640 < repulsive and attractive character of the particles. With a reduced
641 < repulsive surface (the darker region), the particles can move closer
642 < to one another, increasing the density for the overall system. This
643 < change in interaction cutoff also results in a more gradual
644 < orientational motion by allowing the particles to maintain preferred
645 < dipolar arrangements before they begin to feel the pull of the
646 < tetrahedral restructuring. Upon moving closer together, the dipolar
647 < repulsion term becomes active and excludes the unphysical
648 < arrangements. This compares with the original SSD's excluding dipolar
649 < before the particles feel the pull of the ``hydrogen bonds''. Aside
650 < from improving the shape of the first peak in the g(\emph{r}), this
651 < improves the densities considerably by allowing the persistence of
652 < full dipolar character below the previous 4.0 \AA\ cutoff.
629 > g(\emph{r}) by comparing the revised SSD model (SSD1), SSD-E, and
630 > SSD-RF to the new experimental results. Both the modified water models
631 > have shorter peaks that are brought in more closely to the
632 > experimental peak (as seen in the insets of figure \ref{grcompare}).
633 > This structural alteration was accomplished by the combined reduction
634 > in the Lennard-Jones $\sigma$ variable and adjustment of the sticky
635 > potential strength and cutoffs. As can be seen in table \ref{params},
636 > the cutoffs for the tetrahedral attractive and dipolar repulsive terms
637 > were nearly swapped with each other. Isosurfaces of the original and
638 > modified sticky potentials are shown in figure \cite{isosurface}. In
639 > these isosurfaces, it is easy to see how altering the cutoffs changes
640 > the repulsive and attractive character of the particles. With a
641 > reduced repulsive surface (the darker region), the particles can move
642 > closer to one another, increasing the density for the overall
643 > system. This change in interaction cutoff also results in a more
644 > gradual orientational motion by allowing the particles to maintain
645 > preferred dipolar arrangements before they begin to feel the pull of
646 > the tetrahedral restructuring. Upon moving closer together, the
647 > dipolar repulsion term becomes active and excludes unphysical
648 > nearest-neighbor arrangements. This compares with how SSD and SSD1
649 > exclude preferred dipole alignments before the particles feel the pull
650 > of the ``hydrogen bonds''. Aside from improving the shape of the first
651 > peak in the g(\emph{r}), this improves the densities considerably by
652 > allowing the persistence of full dipolar character below the previous
653 > 4.0 \AA\ cutoff.
654  
655   While adjusting the location and shape of the first peak of
656 < g(\emph{r}) improves the densities to some degree, these changes alone
657 < are insufficient to bring the system densities up to the values
658 < observed experimentally. To finish bringing up the densities, the
659 < dipole moments were increased in both the adjusted models. Being a
660 < dipole based model, the structure and transport are very sensitive to
661 < changes in the dipole moment. The original SSD simply used the dipole
662 < moment calculated from the TIP3P water model, which at 2.35 D is
656 > g(\emph{r}) improves the densities, these changes alone are
657 > insufficient to bring the system densities up to the values observed
658 > experimentally. To finish bringing up the densities, the dipole
659 > moments were increased in both the adjusted models. Being a dipole
660 > based model, the structure and transport are very sensitive to changes
661 > in the dipole moment. The original SSD simply used the dipole moment
662 > calculated from the TIP3P water model, which at 2.35 D is
663   significantly greater than the experimental gas phase value of 1.84
664 < D. The larger dipole moment is a more realistic value and improve the
664 > D. The larger dipole moment is a more realistic value and improves the
665   dielectric properties of the fluid. Both theoretical and experimental
666   measurements indicate a liquid phase dipole moment ranging from 2.4 D
667 < to values as high as 3.11 D, so there is quite a range available for
668 < adjusting the dipole
667 > to values as high as 3.11 D, so there is quite a range of available
668 > values for a reasonable dipole
669   moment.\cite{Sprik91,Kusalik02,Badyal00,Barriol64} The increasing of
670 < the dipole moments to 2.418 and 2.48 D for SSD/E and SSD/RF
671 < respectively is moderate in the range of the experimental values;
672 < however, it leads to significant changes in the density and transport
668 < of the water models.
670 > the dipole moments to 2.42 and 2.48 D for SSD/E and SSD/RF
671 > respectively is moderate in this range; however, it leads to
672 > significant changes in the density and transport of the water models.
673  
674 < In order to demonstrate the benefits of this reparameterization, a
674 > In order to demonstrate the benefits of these reparameterizations, a
675   series of NPT and NVE simulations were performed to probe the density
676   and transport properties of the adapted models and compare the results
677   to the original SSD model. This comparison involved full NPT melting
678   sequences for both SSD/E and SSD/RF, as well as NVE transport
679 < calculations at both self-consistent and experimental
680 < densities. Again, the results come from five separate simulations of
681 < 1024 particle systems, and the melting sequences were started from
682 < different ice $I_h$ crystals constructed as stated previously. Like
683 < before, all of the NPT simulations were equilibrated for 100 ps before
684 < a 200 ps data collection run, and they used the previous temperature's
685 < final configuration as a starting point. All of the NVE simulations
686 < had the same thermalization, equilibration, and data collection times
687 < stated earlier in this paper.
679 > calculations at the calculated self-consistent densities. Again, the
680 > results come from five separate simulations of 1024 particle systems,
681 > and the melting sequences were started from different ice $I_h$
682 > crystals constructed as stated earlier. Like before, each NPT
683 > simulation was equilibrated for 100 ps before a 200 ps data collection
684 > run at each temperature step, and they used the final configuration
685 > from the previous temperature simulation as a starting point. All of
686 > the NVE simulations had the same thermalization, equilibration, and
687 > data collection times stated earlier in this paper.
688  
689   \begin{figure}
690   \includegraphics[width=62mm, angle=-90]{ssdeDense.epsi}
691 < \caption{Comparison of densities calculated with SSD/E to SSD without a
691 > \caption{Comparison of densities calculated with SSD/E to SSD1 without a
692   reaction field, TIP3P\cite{Jorgensen98b}, TIP5P\cite{Jorgensen00},
693   SPC/E\cite{Clancy94}, and Experiment\cite{CRC80}. The window shows a
694   expansion around 300 K with error bars included to clarify this region
# Line 693 | Line 697 | Figure \ref{ssdedense} shows the density profile for t
697   \label{ssdedense}
698   \end{figure}
699  
700 < Figure \ref{ssdedense} shows the density profile for the SSD/E water
701 < model in comparison to the original SSD without a reaction field,
702 < experiment, and the other common water models considered
703 < previously. The calculated densities have increased significantly over
704 < the original SSD model and match the experimental value just below 298
705 < K. At 298 K, the density of SSD/E is 0.995$\pm$0.001 g/cm$^3$, which
706 < compares well with the experimental value of 0.997 g/cm$^3$ and is
707 < considerably better than the SSD value of 0.967$\pm$0.003
708 < g/cm$^3$. The increased dipole moment in SSD/E has helped to flatten
709 < out the curve at higher temperatures, only the improvement is marginal
710 < at best. This steep drop in densities is due to the dipolar rather
711 < than charge based interactions which decay more rapidly at longer
712 < distances.
713 <
714 < By monitoring C$\text{p}$ throughout these simulations, the melting
715 < transition for SSD/E was observed at 230 K, about 5 degrees lower than
716 < SSD. The resulting density maximum is located at 240 K, again about 5
717 < degrees lower than the SSD value of 245 K. Though there is a decrease
718 < in both of these values, the corrected densities near room temperature
715 < justify the modifications taken.
700 > Figure \ref{ssdedense} shows the density profile for the SSD/E model
701 > in comparison to SSD1 without a reaction field, experiment, and other
702 > common water models. The calculated densities for both SSD/E and SSD1
703 > have increased significantly over the original SSD model (see figure
704 > \ref{dense1} and are in significantly better agreement with the
705 > experimental values. At 298 K, the density of SSD/E and SSD1 without a
706 > long-range correction are 0.996$\pm$0.001 g/cm$^3$ and 0.999$\pm$0.001
707 > g/cm$^3$ respectively.  These both compare well with the experimental
708 > value of 0.997 g/cm$^3$, and they are considerably better than the SSD
709 > value of 0.967$\pm$0.003 g/cm$^3$. The changes to the dipole moment
710 > and sticky switching functions have improved the structuring of the
711 > liquid (as seen in figure \ref{grcompare}, but they have shifted the
712 > density maximum to much lower temperatures. This comes about via an
713 > increase of the liquid disorder through the weakening of the sticky
714 > potential and strengthening of the dipolar character. However, this
715 > increasing disorder in the SSD/E model has little affect on the
716 > melting transition. By monitoring C$\text{p}$ throughout these
717 > simulations, the melting transition for SSD/E occurred at 235 K, the
718 > same transition temperature observed with SSD and SSD1.
719  
720   \begin{figure}
721   \includegraphics[width=62mm, angle=-90]{ssdrfDense.epsi}
722 < \caption{Comparison of densities calculated with SSD/RF to SSD with a
722 > \caption{Comparison of densities calculated with SSD/RF to SSD1 with a
723   reaction field, TIP3P\cite{Jorgensen98b}, TIP5P\cite{Jorgensen00},
724   SPC/E\cite{Clancy94}, and Experiment\cite{CRC80}. The inset shows the
725   necessity of reparameterization when utilizing a reaction field
# Line 725 | Line 728 | Figure \ref{ssdrfdense} shows a density comparison bet
728   \label{ssdrfdense}
729   \end{figure}
730  
731 < Figure \ref{ssdrfdense} shows a density comparison between SSD/RF and
732 < SSD with an active reaction field. Like in the simulations of SSD/E,
733 < the densities show a dramatic increase over normal SSD. At 298 K,
734 < SSD/RF has a density of 0.997$\pm$0.001 g/cm$^3$, right in line with
735 < experiment and considerably better than the SSD value of
736 < 0.941$\pm$0.001 g/cm$^3$. The melting point is observed at 240 K,
737 < which is 5 degrees lower than SSD with a reaction field, and the
738 < density maximum at 255 K, again 5 degrees lower than SSD. The density
739 < at higher temperature still drops off more rapidly than the charge
740 < based models but is in better agreement than SSD/E.
731 > Including the reaction field long-range correction results in a more
732 > interesting comparison. A density profile including SSD/RF and SSD1
733 > with an active reaction field is shown in figure \ref{ssdrfdense}.  As
734 > observed in the simulations without a reaction field, the densities of
735 > SSD/RF and SSD1 show a dramatic increase over normal SSD (see figure
736 > \ref{dense1}). At 298 K, SSD/RF has a density of 0.997$\pm$0.001
737 > g/cm$^3$, right in line with experiment and considerably better than
738 > the SSD value of 0.941$\pm$0.001 g/cm$^3$ and the SSD1 value of
739 > 0.972$\pm$0.002 g/cm$^3$. These results further emphasize the
740 > importance of reparameterization in order to model the density
741 > properly under different simulation conditions. Again, these changes
742 > don't have that profound an effect on the melting point which is
743 > observed at 245 K for SSD/RF, identical to SSD and only 5 K lower than
744 > SSD1 with a reaction field. However, the difference in density maxima
745 > is not quite as extreme with SSD/RF showing a density maximum at 255
746 > K, fairly close to 260 and 265 K, the density maxima for SSD and SSD1
747 > respectively.
748  
749 + \begin{figure}
750 + \includegraphics[width=65mm, angle=-90]{ssdeDiffuse.epsi}
751 + \caption{Plots of the diffusion constants calculated from SSD/E and SSD1,
752 + both without a reaction field, along with experimental results are
753 + from Gillen \emph{et al.}\cite{Gillen72} and Mills\cite{Mills73}. The
754 + NVE calculations were performed at the average densities observed in
755 + the 1 atm NPT simulations for the respective models. SSD/E is
756 + slightly more fluid than experiment at all of the temperatures, but
757 + it is closer than SSD1 without a long-range correction.}
758 + \label{ssdediffuse}
759 + \end{figure}
760 +
761   The reparameterization of the SSD water model, both for use with and
762   without an applied long-range correction, brought the densities up to
763   what is expected for simulating liquid water. In addition to improving
764   the densities, it is important that particle transport be maintained
765   or improved. Figure \ref{ssdediffuse} compares the temperature
766 < dependence of the diffusion constant of SSD/E to SSD without an active
767 < reaction field, both at the densities calculated at 1 atm and at the
768 < experimentally calculated densities for super-cooled and liquid
766 > dependence of the diffusion constant of SSD/E to SSD1 without an
767 > active reaction field, both at the densities calculated at 1 atm and
768 > at the experimentally calculated densities for super-cooled and liquid
769   water. In the upper plot, the diffusion constant for SSD/E is
770 < consistently a little faster than experiment, while SSD starts off
771 < slower than experiment and crosses to merge with SSD/E at high
772 < temperatures. Both models follow the experimental trend well, but
773 < diffuse too rapidly at higher temperatures. This abnormally fast
774 < diffusion is caused by the decreased system density. Since the
775 < densities of SSD/E don't deviate as much from experiment as those of
776 < SSD, it follows the experimental trend more closely. This observation
777 < is backed up by looking at the lower plot. The diffusion constants for
778 < SSD/E track with the experimental values while SSD deviates on the low
779 < side of the trend with increasing temperature. This is again a product
780 < of SSD/E having densities closer to experiment, and not deviating to
781 < lower densities with increasing temperature as rapidly.
770 > consistently a little faster than experiment, while SSD1 remains
771 > slower than experiment until relatively high temperatures (greater
772 > than 330 K). Both models follow the shape of the experimental trend
773 > well below 300 K, but the trend leans toward diffusing too rapidly at
774 > higher temperatures, something that is especially apparent with
775 > SSD1. This accelerated increasing of diffusion is caused by the
776 > rapidly decreasing system density with increasing temperature. Though
777 > it is difficult to see in figure \ref{ssdedense}, the densities of SSD1
778 > decay more rapidly with temperature than do those of SSD/E, leading to
779 > more visible deviation from the experimental diffusion trend. Thus,
780 > the changes made to improve the liquid structure may have had an
781 > adverse affect on the density maximum, but they improve the transport
782 > behavior of the water model.
783  
784   \begin{figure}
785   \includegraphics[width=65mm, angle=-90]{ssdrfDiffuse.epsi}
# Line 772 | Line 795 | lower densities with increasing temperature as rapidly
795   \label{ssdrfdiffuse}
796   \end{figure}
797  
775 – \begin{figure}
776 – \includegraphics[width=65mm, angle=-90]{ssdeDiffuse.epsi}
777 – \caption{Plots of the diffusion constants calculated from SSD/E and SSD1,
778 – both without a reaction field, along with experimental results are
779 – from Gillen \emph{et al.}\cite{Gillen72} and Mills\cite{Mills73}. The
780 – NVE calculations were performed at the average densities observed in
781 – the 1 atm NPT simulations for the respective models. SSD/E is
782 – slightly more fluid than experiment at all of the temperatures, but
783 – it is closer than SSD1 without a long-range correction.}
784 – \label{ssdediffuse}
785 – \end{figure}
786 –
798   In figure \ref{ssdrfdiffuse}, the diffusion constants for SSD/RF are
799 < compared with SSD with an active reaction field. In the upper plot,
800 < SSD/RF tracks with the experimental results incredibly well, identical
801 < within error throughout the temperature range and only showing a
802 < slight increasing trend at higher temperatures. SSD also tracks
803 < experiment well, only it tends to diffuse a little more slowly at low
804 < temperatures and deviates to diffuse too rapidly at high
805 < temperatures. As was stated in the SSD/E comparisons, this deviation
806 < away from the ideal trend is due to a rapid decrease in density at
807 < higher temperatures. SSD/RF doesn't suffer from this problem as much
808 < as SSD, because the calculated densities are more true to
809 < experiment. This is again emphasized in the lower plot, where SSD/RF
810 < tracks the experimental diffusion exactly while SSD's diffusion
800 < constants are slightly too low due to its need for a lower density at
801 < the specified temperature.
799 > compared with SSD1 with an active reaction field. Note that SSD/RF
800 > tracks the experimental results incredibly well, identical within
801 > error throughout the temperature range shown and only showing a slight
802 > increasing trend at higher temperatures. SSD1 tends to diffuse more
803 > slowly at low temperatures and deviates to diffuse too rapidly at
804 > temperatures greater than 330 K. As was stated in the SSD/E
805 > comparisons, this deviation away from the ideal trend is due to a
806 > rapid decrease in density at higher temperatures. SSD/RF doesn't
807 > suffer from this problem as much as SSD1, because the calculated
808 > densities are more true to experiment. These results again emphasize
809 > the importance of careful reparameterization when using an altered
810 > long-range correction.
811  
812   \subsection{Additional Observations}
813  
805 – While performing the melting sequences of SSD/E, some interesting
806 – observations were made. After melting at 230 K, two of the systems
807 – underwent crystallization events near 245 K. As the heating process
808 – continued, the two systems remained crystalline until finally melting
809 – between 320 and 330 K. These simulations were excluded from the data
810 – set shown in figure \ref{ssdedense} and replaced with two additional
811 – melting sequences that did not undergo this anomalous phase
812 – transition, while this crystallization event was investigated
813 – separately.
814 –
814   \begin{figure}
815   \includegraphics[width=85mm]{povIce.ps}
816 < \caption{Crystal structure of an ice 0 lattice shown from the (001) face.}
816 > \caption{A water lattice built from the crystal structure that SSD/E
817 > assumed when undergoing an extremely restricted temperature NPT
818 > simulation. This form of ice is referred to as ice 0 to emphasize its
819 > simulation origins. This image was taken of the (001) face of the
820 > crystal.}
821   \label{weirdice}
822   \end{figure}
823  
824 < The final configurations of these two melting sequences shows an
825 < expanded zeolite-like crystal structure that does not correspond to
826 < any known form of ice. For convenience and to help distinguish it from
827 < the experimentally observed forms of ice, this crystal structure will
828 < henceforth be referred to as ice-zero (ice 0). The crystallinity was
829 < extensive enough than a near ideal crystal structure could be
830 < obtained. Figure \ref{weirdice} shows the repeating crystal structure
831 < of a typical crystal at 5 K. The unit cell contains eight molecules,
832 < and figure \ref{unitcell} shows a unit cell built from the water
833 < particle center of masses that can be used to construct a repeating
834 < lattice, similar to figure \ref{weirdice}. Each molecule is hydrogen
835 < bonded to four other water molecules; however, the hydrogen bonds are
836 < flexed rather than perfectly straight. This results in a skewed
837 < tetrahedral geometry about the central molecule. Looking back at
838 < figure \ref{isosurface}, it is easy to see how these flexed hydrogen
839 < bonds are allowed in that the attractive regions are conical in shape,
840 < with the greatest attraction in the central region. Though not ideal,
841 < these flexed hydrogen bonds are favorable enough to stabilize an
842 < entire crystal generated around them. In fact, the imperfect ice 0
843 < crystals were so stable that they melted at greater than room
844 < temperature.
824 > While performing restricted temperature melting sequences of SSD/E not
825 > discussed earlier in this paper, some interesting observations were
826 > made. After melting at 235 K, two of five systems underwent
827 > crystallization events near 245 K. As the heating process continued,
828 > the two systems remained crystalline until finally melting between 320
829 > and 330 K. The final configurations of these two melting sequences
830 > show an expanded zeolite-like crystal structure that does not
831 > correspond to any known form of ice. For convenience and to help
832 > distinguish it from the experimentally observed forms of ice, this
833 > crystal structure will henceforth be referred to as ice-zero (ice
834 > 0). The crystallinity was extensive enough that a near ideal crystal
835 > structure could be obtained. Figure \ref{weirdice} shows the repeating
836 > crystal structure of a typical crystal at 5 K. Each water molecule is
837 > hydrogen bonded to four others; however, the hydrogen bonds are flexed
838 > rather than perfectly straight. This results in a skewed tetrahedral
839 > geometry about the central molecule. Looking back at figure
840 > \ref{isosurface}, it is easy to see how these flexed hydrogen bonds
841 > are allowed in that the attractive regions are conical in shape, with
842 > the greatest attraction in the central region. Though not ideal, these
843 > flexed hydrogen bonds are favorable enough to stabilize an entire
844 > crystal generated around them. In fact, the imperfect ice 0 crystals
845 > were so stable that they melted at temperatures nearly 100 K greater
846 > than both ice I$_c$ and I$_h$.
847  
848 < \begin{figure}
844 < \includegraphics[width=65mm]{ice0cell.eps}
845 < \caption{Simple unit cell for constructing ice 0. In this cell, $c$ is
846 < equal to $0.4714\times a$, and a typical value for $a$ is 8.25 \AA.}
847 < \label{unitcell}
848 < \end{figure}
849 <
850 < The initial simulations indicated that ice 0 is the preferred ice
848 > These initial simulations indicated that ice 0 is the preferred ice
849   structure for at least SSD/E. To verify this, a comparison was made
850   between near ideal crystals of ice $I_h$, ice $I_c$, and ice 0 at
851 < constant pressure with SSD/E, SSD/RF, and SSD. Near ideal versions of
852 < the three types of crystals were cooled to ~1 K, and the potential
851 > constant pressure with SSD/E, SSD/RF, and SSD1. Near ideal versions of
852 > the three types of crystals were cooled to 1 K, and the potential
853   energies of each were compared using all three water models. With
854   every water model, ice 0 turned out to have the lowest potential
855 < energy: 5\% lower than $I_h$ with SSD, 6.5\% lower with SSD/E, and
856 < 7.5\% lower with SSD/RF. In all three of these water models, ice $I_c$
859 < was observed to be ~2\% less stable than ice $I_h$. In addition to
860 < having the lowest potential energy, ice 0 was the most expanded of the
861 < three ice crystals, ~5\% less dense than ice $I_h$ with all of the
862 < water models. In all three water models, ice $I_c$ was observed to be
863 < ~2\% more dense than ice $I_h$.
855 > energy: 5\% lower than $I_h$ with SSD1, 6.5\% lower with SSD/E, and
856 > 7.5\% lower with SSD/RF.
857  
858 < In addition to the low temperature comparisons, melting sequences were
859 < performed with ice 0 as the initial configuration using SSD/E, SSD/RF,
860 < and SSD both with and without a reaction field. The melting
861 < transitions for both SSD/E and SSD without a reaction field occurred
862 < at temperature in excess of 375 K. SSD/RF and SSD with a reaction
858 > In addition to these low temperature comparisons, melting sequences
859 > were performed with ice 0 as the initial configuration using SSD/E,
860 > SSD/RF, and SSD1 both with and without a reaction field. The melting
861 > transitions for both SSD/E and SSD1 without a reaction field occurred
862 > at temperature in excess of 375 K. SSD/RF and SSD1 with a reaction
863   field had more reasonable melting transitions, down near 325 K. These
864   melting point observations emphasize how preferred this crystal
865   structure is over the most common types of ice when using these single
# Line 874 | Line 867 | models, it is interesting to speculate on the favorabi
867  
868   Recognizing that the above tests show ice 0 to be both the most stable
869   and lowest density crystal structure for these single point water
870 < models, it is interesting to speculate on the favorability of this
871 < crystal structure with the different charge based models. As a quick
870 > models, it is interesting to speculate on the relative stability of
871 > this crystal structure with charge based water models. As a quick
872   test, these 3 crystal types were converted from SSD type particles to
873   TIP3P waters and read into CHARMM.\cite{Karplus83} Identical energy
874   minimizations were performed on all of these crystals to compare the
# Line 884 | Line 877 | continuation on work studing ice 0 with multipoint wat
877   $I_h$, which was in turn ~3\% lower than ice $I_c$. From these initial
878   results, we would not be surprised if results from the other common
879   water models show ice 0 to be the lowest energy crystal structure. A
880 < continuation on work studing ice 0 with multipoint water models will
880 > continuation on work studying ice 0 with multi-point water models will
881   be published in a coming article.
882  
883   \section{Conclusions}
# Line 897 | Line 890 | density behavior, SSD was reparameterized for use both
890   of other water models. Analysis of particle diffusion showed SSD to
891   capture the transport properties of experimental very well in both the
892   normal and super-cooled liquid regimes. In order to correct the
893 < density behavior, SSD was reparameterized for use both with and
894 < without a long-range interaction correction, SSD/RF and SSD/E
895 < respectively. In addition to correcting the abnormally low densities,
896 < these new versions were show to maintain or improve upon the transport
897 < and structural features of the original water model, all while
898 < maintaining the fast performance of the SSD water model. This work
899 < shows these simple water models, and in particular SSD/E and SSD/RF,
900 < to be excellent choices to represent explicit water in future
901 < simulations of biochemical systems.
893 > density behavior, the original SSD model was reparameterized for use
894 > both with and without a reaction field (SSD/RF and SSD/E), and
895 > comparison simulations were performed with SSD1, the density corrected
896 > version of SSD. Both models improve the liquid structure, density
897 > values, and diffusive properties under their respective conditions,
898 > indicating the necessity of reparameterization when altering the
899 > long-range correction specifics. When taking the appropriate
900 > considerations, these simple water models are excellent choices for
901 > representing explicit water in large scale simulations of biochemical
902 > systems.
903  
904   \section{Acknowledgments}
905   Support for this project was provided by the National Science

Diff Legend

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