Water is an essential component in simulations of biomolecular systems. Water makes up the majority of the molecules involved, and the hydrogen bond networking of the solvent is an essential aspect dictating the phase behavior of lipid systems. In choosing a water model for our simulations, we wanted to a model that was computationally efficient, yet was capable of capturing the electrostatic and hydrogen bonding interactions that exist in real water. We selected the soft-sticky dipole (SSD) model as the explicit solvent that we would utilize in our research.[1] The SSD model differs from other simple water models in that it is a single point model, and it requires one distance calculation for each pair of molecules. For a comparison, TIP3P, the most commonly used water model in biomolecular simulations, is a three site model that requires nine distance calculations per pair. The interaction potential is given as follows:



Thus, the potential is the sum of three pair-interactions (from left to right): Lennard-Jones for van der Waals interactions, dipole instead of charge-charge interactions to represent the electrostatic interactions, and a "sticky" potential to mimic the tetrahedral hydrogen-bonding. The sticky potential consists of the sum of two spherical harmonics as well as an empirical term to counter dipole alignment in the first solvation shell. Switching functions turn these interactions off outside the first solvation shell. These functions can be visualized as follows:




When we performed our first SSD simulations, we discovered that the water model exhibited densities significantly lower than the experimental density of 0.997 g/cm³ at 298 K. These initial results showed that a complete reparametrization of SSD was required, particularly for use with a reaction field (an interaction correction that compensates for truncation effects).[2] Our newly parametrized potentials (SSD/E and SSD/RF) involve a contraction of the Lennard-Jones surface, strengthening of the dipole moment, and modification of the sticky potential to reduce particle repulsion. The effects of our modifications to the sticky potential can be seen in the isosurfaces below.

SSD & SSD1 Isosurface	SSD/E & SSD/RF Isosurface

The gold regions represent the attractive portions of the potential while the blue represents repulsive regions. These reparametrizations resulted in significant improvements in various static and dynamic liquid properties.

Property	SSD1	SSD/E	SSD1 (RF)	SSD/RF	Expt.
Density (g/cm3)	0.999(0.001)	0.996(0.001)	0.972(0.002)	0.997(0.001)	0.997
Cp (cal/mol K)	28.80(0.11)	25.45(0.09)	28.28(0.06)	23.83(0.16)	17.98
D (10-5 cm2/s)	1.78(0.07)	2.51(0.18)	2.00(0.17)	2.32(0.06)	2.299
Coordination Number (nC)	3.9	4.3	3.8	4.4	4.7
H-bonds per particle (nH)	3.7	3.6	3.7	3.7	3.5
tau1 (ps)	10.9(0.6)	7.3(0.4)	7.5(0.7)	7.2(0.4)	5.7
tau2 (ps)	4.7(0.4)	3.1(0.2)	3.5(0.3)	3.2(0.2)	2.3

(SSD1 is a recently developed model that was created solely to correct for the density anomaly in the original SSD.[3])

Crystallization and Ice-i

While parametrizing our new models, several systems of SSD/E crystallized from a super-cooled liquid state in isothermal-isobaric (NPT) simulations. Crystallization events are marked by sudden shifts in the simulation properties, such as a decrease in density or drop in the system potential energy. It is visually apparent that ordered regions developed from the initially disordered system - note the diagonal rows in the crystal on the right in the image below.



Note that the diagonal rows do not propagate through the entire crystal and seem to stop about halfway across the cell. This is due to the fact that the system crystallized well below the melting point, and it became locked in a fractured state - looking along another plane, you would be able to pick out some order in the regions that looked disordered at this angle. From these fractured solid systems, we were able to piece together an ideal crystal. It was determined to be an ice polymorph that has never been observed experimentally (or in other simulations). The following rendering shows what this new ice polymorph, which we are calling Ice-i, looks like when viewing it down on the (001) crystal face:



The crystal is composed of rows of water tetramers that interlock in such a fashion to result in open octagonal cavities. The diameter of these open cavities is on average approximately 6.3 Å. These large empty volumes make the crystal considerably less dense than any other known or observed ice polymorph. In fact, it is ~7.5% less dense than ice I_h, the lowest density variant of ice on record. The unit cell for this structure is quite simple, consisting of 8 molecules in a stacked square arrangement:

Free Energies and Phase Diagrams

Spontaneous crystallization of a solid in constant pressure simulations was a good indicator that the structure is the minimum free energy structure for the SSD/E water model. In order to verify this, as well as to evaluate the significance of this crystal structure for all water models, we needed to determine the free energy and compare it to the free energies of other common ice crystals. By comparing the free energies of the various polymorphs under common pressure and temperature conditions, we can understand which crystalline phase is thermodynamically preferred for that particular model and thus, which form we would expect it to adopt under unrestricted dynamics.

In order to perform this free energy comparison, the computational technique of thermodynamic integration was used on several different ice polymorphs. Thermodynamic integration is an established technique that has been successfully used in the past to evaluate the free energies of both liquid and solid phases of water.[4-6] Since our dynamics code is already based on rigid body orientational propagation, applying the method used by Báez and Clancy was relatively straightforward.[5] Basically, our system is transformed into a system where the partition function can be solved for analytically (an Einstein crystal for solids). We know the free energy at the transformed state and can determine the change in free energy during the transformation; thus, with this information, it is easy to back out the free energy of the initial system. This process was used on five different low pressure ice polymorphs (ice I_h, ice I_c, ice B, Ice-i, and Ice-i´) using six different common water models (TIP3P, TIP4P, TIP5P, SPC/E, SSD/E, and SSD/RF). Results for the initial 9 Å cutoff simulations are shown in the following table:

Water Model	Ih	Ic	B	Ice-i	Ice-i'	Tm (*Ts)	Tb
TIP3P	-11.41(2)	-11.23(3)	-11.82(3)	-12.30(3)	-	269(4)	357(2)
TIP4P	-11.84(3)	-12.04(2)	-12.08(3)	-	-12.33(3)	266(5)	354(2)
TIP5P	-11.85(3)	-11.86(2)	-11.96(2)	-	-12.29(2)	271(4)	337(2)
SPC/E	-12.87(2)	-13.05(2)	-13.26(3)	-	-13.55(2)	296(3)	396(2)
SSD/E	-11.27(2)	-11.19(4)	-12.09(2)	-12.54(2)	-	*355(2)	-
SSD/RF	-11.51(2)	-11.47(2)	-12.08(3)	-12.29(2)	-	278(4)	349(2)

Click on the crystal images for a more detailed view.

In all of these simulations, Ice-i (or Ice-i´) has the lowest free energy for all of the water models studied. Another interesting point is the fact that the melting temperatures (with the obvious exception of SSD/E) show reasonably good agreement with the experimental value of 273 K, i.e. much better than the values observed for TIP4P from ice I_h (from 214 to 238 K).[7-9] An important distinction to note is that these melting transitions are calculated from the most stable phase (Ice-i or Ice-i´ in this case), and calculation of the T_m from ice I_h from this work places TIP4P's melting transition at ~210 K - in line with the estimates from other groups.

With the free energies at one or more state points, we can project out to other temperatures using the Gibbs-Helmholtz equation and other pressures using the pressure dependence of the free energy.[6] This allows us to generate phase diagrams for the various models. With phase diagrams, we can get an idea of what transition is in store for your system under the chosen conditions. Below are phase diagrams for TIP3P and SSD/RF as determined from this data:

TIP3P Phase Diagram	SSD/RF Phase Diagram

Note that the experimentally expected low-density ice polymorph, ice I_h, is not present in either of these phase diagrams. In fact, ice I_h is not present (only metastable) in the phase diagrams for any of the water models studied using a 9.0 Å cutoff.

Additional tests were performed in order to gauge the effect of the long-range interaction truncation on the free energy results. The cutoff radii in simulations of the most computationally efficient water models (3 particles or less) were increased systematically to determine the free energy dependence. In addition, a long range correction (in the form of a reaction field or Ewald summation) was incorporated to estimate the effect of infinite periodicity. The results of these tests are plotted in the figure on the right.

Note that Ice-i loses some of its relative stability when moving out to longer cutoff radii in nearly all cases. In fact, with the SPC/E model, each of the studied polymorphs spends some time as the minimum energy state at differing cutoff radii. In spite of this initial look, Ice-i looks to be the minimum energy polymorph regardless of conditions for TIP3P and the single point water models. It is also a stable polymorph for SPC/E, and the differences of free energy between the various polymorphs that simulation conditions (such as the volume in canonical ensemble simulations) may have an effect in what polymorph is expressed during crystallization.

References

[1] Liu, Y.; Ichiye, T. J. Phys. Chem. 1996, 100, 2723-2730.

[2] Fennell, C. J.; Gezelter, J. D. J. Chem. Phys. 2004, 120, 9175-9184.

[3] Tan, M.-L.; Fischer, J. T.; Chandra, A.; Brooks, B. R.; Ichiye, T. Chem. Phys. Lett. 2003, 376, 646-652.

[4] Hermans, J.; Pathiaseril, A.; Anderson, A.; J. Am. Chem. Soc. 1988, 110, 5982-5986.

[5] Báez, L. A.; Clancy, P. Mol. Phys. 1995, 86, 385-396.

[6] Báez, L. A.; Clancy, P. J. Chem. Phys. 1995, 103, 9744-9755.

[7] Sanz, E.; Vega, C.; Abascal, J. L. F.; MacDowell, L. G. Phys. Rev. Lett. 2004, 92, 255701.

[8] Vlot, M. J.; Huinink, J.; van der Eerden, J. P. J. Chem. Phys. 1999, 110, 55-61.

[9] Gao, G. T.; Zeng, X. C.; Tanaka, H. J. Chem. Phys. 2000, 112, 8534-8538.