Dynamical properties of the soft sticky dipole-quadrupole-octupole water model: a molecular dynamics study

The dynamical properties of the soft sticky dipole-quadrupole-octupole (SSDQO) water model using SPC/E moments are calculated utilizing molecular dynamics simulations. This new potential for liquid water describes the water-water interactions by a Lennard-Jones term and a sticky potential, which is an approximate moment expansion with point dipole, quadrupole, and octupole moments, and reproduces radial distribution functions of pure liquid water using the moments of SPC/E [Ichiye and Tan, J. Chem. Phys. 124, 134504 (2006)]. The forces and torques of SSDQO water for the dipole-quadrupole, quadrupole-quadrupole, and dipole-octupole interactions are derived here. The simulations are carried out at 298 K in the microcanonical ensemble employing the Ewald method for the long-range dipole-dipole interactions. Here, various dynamical properties associated with translational and rotational motions of SSDQO water using the moments of SPC/E (SSDQO:SPC/E) water are compared with the results from SPC/E and also experiment. The self-diffusion coefficient of SSDQO:SPC/E water is found to be in excellent agreement with both SPC/E and experiment whereas the single particle orientational relaxation time for dipole vector is better than SPC/E water but it is somewhat smaller than experiment. The dielectric constant of SSDQO:SPC/E is essentially identical to SPC/E, and both are slightly lower than experiment. Also, molecular dynamics simulations of the SSDQO water model are found to be about twice as fast as three-site models such as SPC/E.     

Study of multipole contributions to the structure of water around ions in solution using the soft sticky dipole-quadrupole-octupole (SSDQO) model of water

The solvation of ions in the soft sticky dipole-quadrupole-octupole (SSDQO) model for liquid water is presented here. This new potential energy function for liquid water describes water-water interactions by a Lennard-Jones term plus a sticky potential consisting of an approximate moment expansion with point dipole, quadrupole, and octupole moments. The SSDQO potential energy function using the moments from extended simple point charge (SPC/E), TIP3P, or TIP5P reproduces the pair potential energy functions and radial distribution functions of the respective multipoint model but it is much faster than even the three-point models. Here, the solvation of ions in SSDQO water is studied using ion-water potential energy functions consisting of moment expansions up to the charge-quadrupole term, up to the charge-octupole term, and up to an approximate charge-hexadecapole term using the moments of SPC/E water. The radial distributions from Monte Carlo simulations show the best agreement with the results for ions in SPC/E water for the expansion up to the charge-hexadecapole term. Thus, the best results are obtained when the water-water and ion-water potentials are exact up to the 1r(4) term and also contain an approximate 1r(5) term. Overall, the simplicity, efficiency, and accuracy of the SSDQO potential make it potentially very useful for computer simulations of aqueous solvation.     

Molecular dynamics simulation of the density and surface tension of water by particle-particle particle-mesh method

In this work, molecular dynamics simulation is performed to study the density and surface tension of water for a range of temperatures from 300 to 600 K. The extended simple point charge interaction potential for water is used. The particle-particle particle-mesh method, which automatically includes untruncated long-range terms, is used for the Lennard-Jones and the Coulombic terms. The results show that the long-range correction for the Lennard-Jones term is very important for the calculation of surface tension. It is found that the calculated density and surface tension of water fit well with experimental data for temperatures less than 500 K. Near the critical temperature, the simulation results are off from the experimental data.     

Computational investigation of order, structure, and dynamics in modified water models

Model liquids have been constructed to study the role of local structure in the anomalous properties of liquid water. The intermolecular potentials were modified by increasing the weight of the Lennard-Jones term relative to the electrostatic term in the SPC/E model for water. The resulting family of liquids varies from SPC/E water to a Lennard-Jones-like liquid. Properties were measured as a function of density and temperature. The local structure was described by two order parameters, one measuring the tetrahedral order and the other measuring the translational order. The translational order parameter was found to be large for both tetrahedral and Lennard-Jones liquids, but to go through a minimum as the potentials were modified, demonstrating that the two types of structure are incompatible. Just as in water several properties (e.g., the translational diffusion coefficient, entropy) exhibit anomalous density dependence as a result of the breakdown of local tetrahedrality, we observed nonmonotonic behavior of the translational diffusion constant and reorientational relaxation rate as the fluids were transformed from tetrahedral to Lennard-Jones-like. This is also an indication of the incompatibility between Lennard-Jones and water-like structure.     

Structure and dynamics of liquid water with different long-range interaction truncation and temperature control methods in molecular dynamics simulations

We have used molecular dynamics simulations to study the physical properties of modified TIP3P water model included in the CHARMM program, using four different methods-the Ewald summation technique, and three different spherical truncation methods-for the treatment of the long-range interactions. Both the structure and dynamics of the liquid water model were affected by the methods used to truncate the long-range interactions. For some of the methods artificial structuring of the model liquid was observed around the cutoff radius. The model liquid properties were also affected by the commonly applied temperature control methods. Four different methods for controlling the temperature of the system were studied, and the effects of these methods on the bulk properties for liquid water were analyzed. The system size was also found to change the dynamics of the model liquid water. Two control simulations with the SPC/E water model were carried out. The self-diffusion coefficient (D), the radial distribution function (g(OO)), the distance dependent Kirkwood G-factor [G(k)(r)] and the intermolecular potential energy (E(pot)) were determined from the different trajectories and compared with the experimental data.     

Formation of FeCl(2)/NaCl-nanoparticles in supercritical water investigated by molecular dynamics simulations: nucleation rates

Nucleation and growth of FeCl(2) in supercritical water containing NaCl at different state points between temperatures of 798 and 873 K and system densities between 0.24 and 0.14 g cm(-3) have been studied by molecular dynamics simulations. The number of NaCl ion pairs was chosen to simulate particle formation in seawater and brine of higher salinity. Rigid SPC/E water was used to model the water molecules while a combination of Coulomb and Lennard-Jones potentials was used for the ions. Two different methods for determination of nucleation rates are applied and their results compared. We find decreasing nucleation rates with both increasing temperature and decreasing system density. Our results are also compared to those we recently obtained in an investigation of pure FeCl(2) from supercritical water. We find both increasing nucleation rates and a decreasing size of the critical cluster with increasing amount of NaCl.     

Hydrophobic solvation of Gay-Berne particles in modified water models

The solvation of large hydrophobic solutes, modeled as repulsive and attractive Gay-Berne oblate ellipsoids, is characterized in several modified water liquids using the SPC/E model as the reference water fluid. We find that small amounts of attraction between the Gay-Berne particle and any model fluid result in wetting of the hydrophobic surface. However, significant differences are found among the modified and SPC/E water models and the critical distances in which they dewet the hydrophobic surfaces of pairs of repulsive Gay-Berne particles. We find that the dewetting trends for repulsive Gay-Berne particles in the various model liquids correlate directly with their surface tensions, the widths of the interfaces they form, and the openness of their network structure. The largest critical separations are found in liquids with the smallest surface tensions and the broadest interfaces as measured by the Egelstaff-Widom length.     

Soft sticky dipole-quadrupole-octupole potential energy function for liquid water: an approximate moment expansion

A new, efficient potential energy function for liquid water is presented here. The new model, which is referred here as the soft sticky dipole-quadrupole-octupole (SSDQO) model, describes a water molecule as a Lennard-Jones sphere with point dipole, quadrupole, and octupole moments. It is a single-point model and resembles the hard-sphere sticky dipole potential model for water by Bratko et al. [J. Chem. Phys. 83, 6367 (1985)] and the soft sticky dipole model by Ichiye and Liu [J. Phys. Chem. 100, 2723 (1996)] except now the sticky potential consists of an approximate moment expansion for the dimer interaction potential, which is much faster than the true moment expansion. The object here is to demonstrate that the SSDQO potential energy function can accurately mimic the potential energy function of a multipoint model using the moments of that model. First, the SSDQO potential energy function using the dipole, quadruple, and octupole moments from SPC/E, TIP3P, or TIP5P is shown to reproduce the dimer potential energy functions of the respective multipoint model. In addition, in Monte Carlo simulations of the pure liquid at room temperature, SSDQO reproduces radial distribution functions of the respective model. However, the Monte Carlo simulations using the SSDQO model are about three times faster than those using the three-point models and the long-range interactions decay faster for SSDQO (1/r(3) and faster) than for multipoint models (1/r). Moreover, the contribution of each moment to the energetics and other properties can be determined. Overall, the simplicity, efficiency, and accuracy of the SSDQO potential energy function make it potentially very useful for studies of aqueous solvation by computer simulations.     

Higher-order virial coefficients of water models

We use the Mayer sampling method, with both direct and overlap sampling, to calculate and compare classical virial coefficients up to B6 for various water models (SPC, SPC/E, MSPC/E, TIP3P, and TIP4P). The precision of the computed values ranges from 0.1% for B2 to an average of 25% for B6. When expressed in a form scaled by the critical properties, the values of the coefficients for SPC water are observed to greatly exceed the magnitude of corresponding coefficients for the simple Lennard-Jones model. We examine the coefficients in the context of the equation of state and the Joule-Thomson coefficient. Comparisons of these properties are made both to established molecular simulation data for each respective model and to real water. For all models, the virial series up to B5 describes the equation of state along the saturated vapor line better than the series that includes B6. At supercritical temperatures, however, the sixth-order series often describes pressure-volume-temperature behavior better than the fifth-order series. For example, the sixth-order virial equation of state for SPC/E water predicts the 673 K isotherm within 8% of published molecular simulation values up to a density of 9 mol/L (roughly half the critical density of SPC/E water).     

Development of new Cd2+ and Pb2+ Lennard-Jones parameters for liquid simulations

We present new Lennard-Jones parameters for Cd2+ and Pb2+ ion-water interactions and describe a general methodology to obtain these parameters for any ion. Our strategy is based on the adjustment of ion parameters to reproduce simultaneously experimental absolute hydration free energy and structural properties, namely, g(r) and coordination numbers, obtained from X-ray liquid scattering and quantum mechanical/molecular mechanical (QM/MM) calculations. The validation of the obtained parameters is made by the calculation of dynamical properties and comparing them with experimental values and theoretical results from the literature. The transferability of parameters is checked by the calculation of thermodynamic, structural, and dynamical properties cited above with four different water models. The results obtained for Cd2+ and Pb2+ show an overall agreement with reference values. The absolute hydration free energy calculated with the TIP3P, SPC/E, SPC, and TIP4P water models presents, respectively, percent differences of 3.8, 3.0, 4.3, and 7.2% for lead(II) and 9.8, 8.4, 10.2, and 14.1% for cadmium(II) when compared with experimental values. Ion-water mean distance and coordination numbers for the first coordination shell are in good agreement with experimental and QM/MM results for both ions. Cd2+ shows a lesser diffusion coefficient compared to that of Pb2+ despite its smaller atomic radius, indicating a more persistent first coordination shell for the cadmium(II) ion, a result confirmed with calculations of the mean residence time of water molecules in the first coordination shell.     

Dynamical and energetic properties of hydrogen and hydrogen-tetrahydrofuran clathrate hydrates

Classical equilibrium molecular dynamics (MD) simulations have been performed to investigate the dynamical and energetic properties in hydrogen and mixed hydrogen-tetrahydrofuran sII hydrates at 30 and 200 K and 0.05 kbar, and also at intermediate temperatures, using SPC/E and TIP4P-2005 water models. The potential model is found to have a large impact on overall density, with the TIP4P-2005 systems being on average 1% more dense than their SPC/E counterparts, due to the greater guest-host interaction energy. For the lightly-filled mixed H(2)-THF system, in which there is single H(2) occupation of the small cage (1s1l), we find that the largest contribution to the interaction energy of both types of guest is the van der Waals component with the surrounding water molecules in the constituent cavities. For the more densely-filled mixed H(2)-THF system, in which there is double H(2) occupation in the small cage (2s1l), we find that there is no dominant component (i.e., van der Waals or Coulombic) in the H(2) interaction energy with the rest of the system, but for the THF molecules, the dominant contribution is again the van der Waals interaction with the surrounding cage-water molecules; again, the Coulombic component increases in importance with increasing temperature. The lightly-filled pure H(2) hydrate (1s4l) system exhibits a similar pattern vis-à-vis the H(2) interaction energy as for the lightly-filled mixed H(2)-THF system, and for the more densely-filled pure H(2) system (2s4l), there is no dominant component of interaction energy, due to the multiple occupancy of the cavities. By consideration of Kubic harmonics, there is some evidence of preferential alignment of the THF molecules, particularly at 200 K; this was found to arise at higher temperatures due to transient hydrogen bonding of the oxygen atom in THF molecules with the surrounding cage-water molecules.     

Computation of methodology-independent single-ion solvation properties from molecular simulations. IV. Optimized Lennard-Jones interaction parameter sets for the alkali and halide ions in water

The raw single-ion solvation free energies computed from atomistic (explicit-solvent) simulations are extremely sensitive to the boundary conditions and treatment of electrostatic interactions used during these simulations. However, as shown recently [M. A. Kastenholz and P. H. Hu?nenberger, J. Chem. Phys. 124, 224501 (2006); M. M. Reif and P. H. Hu?nenberger, J. Chem. Phys. 134, 144103 (2010)], the application of appropriate correction terms permits to obtain methodology-independent results. The corrected values are then exclusively characteristic of the underlying molecular model including in particular the ion-solvent van der Waals interaction parameters, determining the effective ion size and the magnitude of its dispersion interactions. In the present study, the comparison of calculated (corrected) hydration free energies with experimental data (along with the consideration of ionic polarizabilities) is used to calibrate new sets of ion-solvent van der Waals (Lennard-Jones) interaction parameters for the alkali (Li(+), Na(+), K(+), Rb(+), Cs(+)) and halide (F(-), Cl(-), Br(-), I(-)) ions along with either the SPC or the SPC/E water models. The experimental dataset is defined by conventional single-ion hydration free energies [Tissandier et al., J. Phys. Chem. A 102, 7787 (1998); Fawcett, J. Phys. Chem. B 103, 11181] along with three plausible choices for the (experimentally elusive) value of the absolute (intrinsic) hydration free energy of the proton, namely, ΔG(hyd)(?)[H(+)] = -1100, -1075 or -1050 kJ mol(-1), resulting in three sets L, M, and H for the SPC water model and three sets L(E), M(E), and H(E) for the SPC/E water model (alternative sets can easily be interpolated to intermediate ΔG(hyd)(?)[H(+)] values). The residual sensitivity of the calculated (corrected) hydration free energies on the volume-pressure boundary conditions and on the effective ionic radius entering into the calculation of the correction terms is also evaluated and found to be very limited. Ultimately, it is expected that comparison with other experimental ionic properties (e.g., derivative single-ion solvation properties, as well as data concerning ionic crystals, melts, solutions at finite concentrations, or nonaqueous solutions) will permit to validate one specific set and thus, the associated ΔG(hyd)(?)[H(+)] value (atomistic consistency assumption). Preliminary results (first-peak positions in the ion-water radial distribution functions, partial molar volumes of ionic salts in water, and structural properties of ionic crystals) support a value of ΔG(hyd)(?)[H(+)] close to -1100 kJ·mol(-1).     

The importance of polarizability in the modeling of solubility: quantifying the effect of solute polarizability on the solubility of small nonpolar solutes in popular models of water

In recent work by Paschek [J. Chem. Phys. 120, 6674 (2004)] and others [see H. Docherty et al., J. Chem. Phys. 125, 074510 (2006) for a review] it has been suggested that, when coupled to a simple Lennard-Jones model for various small nonpolar solute molecules, the most common models of water (e.g., SPC/E and TIP4P) fail to reproduce quantitatively the solubility of small nonpolar solute molecules in water due in part to failing to account for polarization of the solute molecule. Given the importance of such systems as test-case prototype models of the solubility of proteins and biomolecules, in this work, we investigate the impact of using a polarizable solute model with the SPC/E, TIP3P, TIP4P, TIP4P-Ew, and TIP4P/2005 rigid water models. Specifically we consider Ne, Ar, Kr, Xe, and methane as solutes. In all cases we observe that the use of a polarizable solute improves agreement between experiment and simulations, with the best agreement seen for the largest solutes, Kr, CH(4), and Xe and the modern reparametrizations of the TIP4P model, i.e., the TIP4P-Ew and TIP4P/2005 models.     

Thermodynamic properties of water in the lattice gas model with сonsideration of the vibrational motions of molecules

A molecular model developed earlier for a polar fluid within the lattice gas model is supplemented by considering the vibrational motions of molecules using water as an example. A combination of point dipole and Lennard-Jones potentials from SPC parametrization is chosen as the force field model for the molecule. The main thermodynamic properties of liquid water (density, internal energy, and entropy) are studied as functions of temperature. There is qualitative agreement between the calculation results and the experimental data. Ways of refining the molecular theory are discussed.     

Determination of alkali and halide monovalent ion parameters for use in explicitly solvated biomolecular simulations

Alkali (Li(+), Na(+), K(+), Rb(+), and Cs(+)) and halide (F(-), Cl(-), Br(-), and I(-)) ions play an important role in many biological phenomena, roles that range from stabilization of biomolecular structure, to influence on biomolecular dynamics, to key physiological influence on homeostasis and signaling. To properly model ionic interaction and stability in atomistic simulations of biomolecular structure, dynamics, folding, catalysis, and function, an accurate model or representation of the monovalent ions is critically necessary. A good model needs to simultaneously reproduce many properties of ions, including their structure, dynamics, solvation, and moreover both the interactions of these ions with each other in the crystal and in solution and the interactions of ions with other molecules. At present, the best force fields for biomolecules employ a simple additive, nonpolarizable, and pairwise potential for atomic interaction. In this work, we describe our efforts to build better models of the monovalent ions within the pairwise Coulombic and 6-12 Lennard-Jones framework, where the models are tuned to balance crystal and solution properties in Ewald simulations with specific choices of well-known water models. Although it has been clearly demonstrated that truly accurate treatments of ions will require inclusion of nonadditivity and polarizability (particularly with the anions) and ultimately even a quantum mechanical treatment, our goal was to simply push the limits of the additive treatments to see if a balanced model could be created. The applied methodology is general and can be extended to other ions and to polarizable force-field models. Our starting point centered on observations from long simulations of biomolecules in salt solution with the AMBER force fields where salt crystals formed well below their solubility limit. The likely cause of the artifact in the AMBER parameters relates to the naive mixing of the Smith and Dang chloride parameters with AMBER-adapted Aqvist cation parameters. To provide a more appropriate balance, we reoptimized the parameters of the Lennard-Jones potential for the ions and specific choices of water models. To validate and optimize the parameters, we calculated hydration free energies of the solvated ions and also lattice energies (LE) and lattice constants (LC) of alkali halide salt crystals. This is the first effort that systematically scans across the Lennard-Jones space (well depth and radius) while balancing ion properties like LE and LC across all pair combinations of the alkali ions and halide ions. The optimization across the entire monovalent series avoids systematic deviations. The ion parameters developed, optimized, and characterized were targeted for use with some of the most commonly used rigid and nonpolarizable water models, specifically TIP3P, TIP4P EW, and SPC/E. In addition to well reproducing the solution and crystal properties, the new ion parameters well reproduce binding energies of the ions to water and the radii of the first hydration shells.     

Communication: Shifted forces in molecular dynamics

Simulations involving the Lennard-Jones potential usually employ a cutoff at r = 2.5σ. This communication investigates the possibility of reducing the cutoff. Two different cutoff implementations are compared, the standard shifted potential cutoff and the less commonly used shifted forces cutoff. The first has correct forces below the cutoff, whereas the shifted forces cutoff modifies Newton's equations at all distances. The latter is nevertheless superior; we find that for most purposes realistic simulations may be obtained using a shifted forces cutoff at r = 1.5σ, even though the pair force is here 30 times larger than at r = 2.5σ.     

Chemical potential perturbation: extension of the method to lattice sum treatment of intermolecular potentials

The recently developed chemical potential perturbation (CPP) method [S. G. Moore and D. R. Wheeler, J. Chem. Phys. 134, 114514 (2011)] is extended to the lattice (Ewald) sum treatment of intermolecular potentials. The CPP method predicts chemical potentials for a range of composition points using the local (position-dependent) pressure tensor of an inhomogeneous system. When computing the local pressure tensor, one can use the Irving-Kirkwood (IK) or Harasima (H) contours of distributing the pressure. We compare these two contours and show that for a planar interface, the homogeneous pressure and resulting chemical potential can be approximated with the CPP method using either the IK or the H contour, though with the lattice sum method the H contour has much greater computational efficiency. The proposed methods are validated by calculating the chemical potentials of the Lennard-Jones fluid and extended simple point-charge (SPC/E) water, and results show a high level of agreement with respective equations of state.     

A novel approach for designing simple point charge models for liquid water with three interaction sites

A simultaneous improvement of the diffusion and dielectric properties of the simple point charge (SPC) model for liquid water appears to be very difficult with conventional reparametrization of the commonly used Lennard-Jones and Coulomb interaction functions and without including a self-energy correction in the effective pair-potential as is done in the SPC/E model. Here, a different approach to circumvent this problem is presented. A short-range interaction term, which corrects the oxygen-oxygen energy at small distances by small amounts of energy, was introduced in the nonbonded interaction function. This additional force-field term allows to derive new parameter sets for SPC-like water models that yield better agreement with experimental data on liquid water. Based on previous investigations of the force-field parameter dependence of the water properties of SPC-like models, the necessary parameter changes to obtain a lower diffusion coefficient and a larger dielectric permittivity were specified and accordingly six new models were developed. They all represent an improvement over SPC in terms of structural and diffusional properties, four of them show better dielectric properties also. One model, SPC/S, has been characterized in more detail, and represents most properties of liquid water better than SPC while avoiding the larger discrepancies with experimental values regarding density, thermal compressibility, energy, and free energy of the SPC/E model. We conclude that the use of a simple, short-ranged additional oxygen-oxygen interaction term makes a simultaneous improvement of the diffusion coefficient and the dielectric properties of water feasible.     

Interfacial water at hydrophobic and hydrophilic surfaces: depletion versus adsorption

The structure of the water-solid interface for widely varying surface properties is investigated with Monte Carlo simulations using the SPC/E water model. Of particular interest is the relation between the wetting coefficient as a measure of the hydrophobicity of the substrate and the density depletion close to the solid surface. The substrates are modeled as rigid ordered lattices of sites that interact with water molecules through an orientation-independent Lennard-Jones potential of varying strength. Hydrophilic character is obtained by addition of polar hydroxyl groups on the substrate surface, and the influence of density, spatial distribution, and angular orientation of the polar groups on the interfacial water structure is studied.     

Molecular simulation of the hydration Gibbs energy of barbiturates

In the present work, molecular dynamics calculations of the Gibbs energy of hydration of 10 different substituted barbiturates in SPC/E water were performed using thermodynamic integration. Given that experimental determination of the Gibbs hydration energy for this class of compounds is currently unfeasible, computer simulations appear as the only alternative for the estimation of this important quantity. Several simulation parameters are discussed and optimized based on calculations for barbituric acid. It is concluded that accounting for electrostatic interactions with the Reaction-Field method can be up to two times faster than with Particle-Mesh-Ewald method, without loss of accuracy. Different number of solvent molecules and simulation lengths were also tested. Lennard-Jones and electrostatic contributions were scaled down to zero in an independent way. It is shown that the electrostatic contribution is dominant (representing approximately 90% of the total Gibbs energy of hydration) and that barbiturate intra-molecular interactions cannot be neglected. The importance of the electrostatic contribution is attributed to the formation of hydrogen bonds between the barbiturates and water, which play an important role in the solvation process. The influence of the different substituents and their contribution to the Gibbs energy of hydration was assessed. Finally, the Lennard-Jones contributions and the total hydration Gibbs energy can both be correlated against molecular weight or partition coefficient data for mono- and di-substituted barbiturates.