Empirical Water Models Research Articles

Assessment of an Electrostatic Energy-Based Charge Model for Modeling the Electrostatic Interactions in Water Solvent.

The protein force field based on the restrained electrostatic potential (RESP) charges has limitations in accurately describing hydrogen bonding interactions in proteins. To address this issue, we propose an alternative approach called the electrostatic energy-based charges (EEC) model, which shows improved performance in describing electrostatic interactions (EIs) of hydrogen bonds in proteins. In this study, we further investigate the performance of the EEC model in modeling EIs in water solvent. Our findings demonstrate that the fixed EEC model can effectively reproduce the quantum mechanics/molecular mechanics (QM/MM)-calculated EIs between a water molecule and various water solvent environments. However, to achieve the same level of computational accuracy, the electrostatic potential (ESP) charge model needs to fluctuate according to the electrostatic environment. Our analysis indicates that the requirement for charge adjustments depends on the specific mathematical and physical representation of EIs as a function of the environment for deriving charges. By comparing with widely used empirical water models calibrated to reproduce experimental properties, we confirm that the performance of the EEC model in reproducing QM/MM EIs is similar to that of general purpose TIP4P-like water models such as TIP4P-Ew and TIP4P/2005. When comparing the computed 10,000 distinct EI values within the range of -40 to 0 kcal/mol with the QM/MM results calculated at the MP2/aug-cc-pVQZ/TIP3P level, we noticed that the mean unsigned error (MUE) for the EEC model is merely 0.487 kcal/mol, which is remarkably similar to the MUE values of the TIP4P-Ew (0.63 kcal/mol) and TIP4P/2005 (0.579 kcal/mol) models. However, both the RESP method and the TIP3P model exhibit a tendency to overestimate the EIs, as evidenced by their higher MUE values of 1.761 and 1.293 kcal/mol, respectively. EEC-based molecular dynamics simulations have demonstrated that, when combined with appropriate van der Waals parameters, the EEC model can closely reproduce oxygen-oxygen radial distribution function and density of water, showing a remarkable similarity to the well-established TIP4P-like empirical water models. Our results demonstrate that the EEC model has the potential to build force fields with comparable accuracy to more sophisticated empirical TIP4P-like water models.

A coarse-grained deep neural network model for liquid water

We introduce a coarse-grained deep neural network (CG-DNN) model for liquid water that utilizes 50 rotational and translational invariant coordinates and is trained exclusively against energies of ∼30 000 bulk water configurations. Our CG-DNN potential accurately predicts both the energies and the molecular forces of water, within 0.9 meV/molecule and 54 meV/Å of a reference (coarse-grained bond-order potential) model. The CG-DNN water model also provides good prediction of several structural, thermodynamic, and temperature dependent properties of liquid water, with values close to those obtained from the reference model. More importantly, CG-DNN captures the well-known density anomaly of liquid water observed in experiments. Our work lays the groundwork for a scheme where existing empirical water models can be utilized to develop a fully flexible neural network framework that can subsequently be trained against sparse data from high-fidelity albeit expensive beyond-DFT calculations.

Systematic Design-of-Experiments, factorial-design approaches for tuning simple empirical water models

ABSTRACT Using systematic Design-of-Experiments, factorial-design approaches, we identify how five variables (site charge, sigma, epsilon, OH bond length and HOH bond angle) affect the diffusivity and radial distribution function accuracy within the framework of the Simple Point Charge water model. From this (by a full-functional DoX-effects model), we glean response equations for the effects of these systematic variations on both estimation of site-site radial distribution functions, in addition to self-diffusivities, allowing scope to estimate more precisely optimal parameters for use in tailored force-fields for particular purposes (e.g. improved fidelity towards particular desired structural or dynamical properties). To this end, the diffusivity error was reduced from an initial value of 35.2–0.64%, whilst the radial distribution function was more difficult to improve, with a reduction in error of 2.83%.

Water models for interfacial water simulations

Interfacial water is of fundamental importance in many technological fields, such as biological processes, chemical reactions and lubrications. A prevalent way to study the structure and dynamics of interfacial water is carrying out molecular dynamics simulations with empirical potential water models. However, discrepant results have been reported due to their different charge geometries and target properties. Here we investigated the interfacial water structures on smooth surfaces of varying hydrophobicity at low temperatures by comprehensive molecular dynamics simulations with the prevailing water models. It is shown that the choice of the water model can significantly change the water structure on the hydrophilic surface, while has a minor effect on the contact angle on a hydrophobic surface. Furthermore, zero-dimensional ice pyramids and one-dimensional icicles were formed under the regulation of external charges injected to the substrate or imposed electric field, respectively. These results offer new insights into the water structures on different surfaces and reasonable choice of parameters in molecular simulations, and the development of water models.

Density anomaly of water at negative pressures from first principles

Using molecular dynamics simulations based on ab initio trained high-dimensional neural network potentials, we study the equation of state of liquid water at negative pressures. From density isobars computed for various pressures down to p = −230 MPa we determine the line of density maxima for two potentials based on the BLYP and the RPBE functionals, respectively. In both cases, dispersion corrections are included to account for non-local long-range correlations that give rise to van der Waals forces. We have followed the density maximum down to negative pressures close to the spinodal instability. For both functionals, the temperature of maximum density increases with decreasing pressure under moderate stretching, but changes slope at MPa and MPa for BLYP and RPBE, respectively. Our calculations confirm qualitatively the retracing shape of the line of density maxima found for empirical water models, indicating that the spinodal line maintains a positive slope even at strongly negative pressures.

Communication: Role of explicit water models in the helix folding/unfolding processes

In the last years, it has become evident that computer simulations can assume a relevant role in modelling protein dynamical motions for their ability to provide a full atomistic image of the processes under investigation. The ability of the current protein force-fields in reproducing the correct thermodynamics and kinetics systems behaviour is thus an essential ingredient to improve our understanding of many relevant biological functionalities. In this work, employing the last developments of the metadynamics framework, we compare the ability of state-of-the-art all-atom empirical functions and water models to consistently reproduce the folding and unfolding of a helix turn motif in a model peptide. This theoretical study puts in evidence that the choice of the water models can influence the thermodynamic and the kinetics of the system under investigation, and for this reason cannot be considered trivial.

Assessing the accuracy of improved force-matched water models derived from Ab initio molecular dynamics simulations.

The accuracy of water models derived from ab initio molecular dynamics simulations by means on an improved force-matching scheme is assessed for various thermodynamic, transport, and structural properties. It is found that although the resulting force-matched water models are typically less accurate than fully empirical force fields in predicting thermodynamic properties, they are nevertheless much more accurate than generally appreciated in reproducing the structure of liquid water and in fact superseding most of the commonly used empirical water models. This development demonstrates the feasibility to routinely parametrize computationally efficient yet predictive potential energy functions based on accurate ab initio molecular dynamics simulations for a large variety of different systems. © 2016 Wiley Periodicals, Inc.

Assessing the Performance of the Diffusion Monte Carlo Method As Applied to the Water Monomer, Dimer, and Hexamer.

The diffusion Monte Carlo (DMC) method is applied to the water monomer, dimer, and hexamer using q-TIP4P/F, one of the most simple empirical water models with flexible monomers. The bias in the time step (Δτ) and population size (Nw) is investigated. For the binding energies, the bias in Δτ cancels nearly completely, whereas a noticeable bias in Nw remains. However, for the isotope shift (e.g, in the dimer binding energies between (H2O)2 and (D2O)2), the systematic errors in Nw do cancel. Consequently, very accurate results for the latter (within ∼0.01 kcal/mol) are obtained with moderate numerical effort (Nw ∼ 10(3)). For the water hexamer and its (D2O)6 isotopomer, the DMC results as a function of Nw are examined for the cage and prism isomers. For a given isomer, the issue of the walker population leaking out of the corresponding basin of attraction is addressed by using appropriate geometric constraints. The population size bias for the hexamer is more severe, and to maintain accuracy similar to that of the dimer, Nw must be increased by ∼2 orders of magnitude. Fortunately, when the energy difference between the cage and prism is taken, the biases cancel, thereby reducing the systematic errors to within ∼0.01 kcal/mol when using a population of Nw = 4.8 × 10(5) walkers. Consequently, a very accurate result for the isotope shift is also obtained. Notably, both the quantum and isotope effects for the prism-cage energy difference are small.

The curious case of the water hexamer: Cage vs. Prism

Small water clusters, such as the hexamer, provide a unique opportunity to advance the molecular-level understanding of water in all its phases. In particular, the water hexamer is the smallest cluster that possesses several nearly iso-energetic non-planar isomeric forms whose relative stability at low temperatures can be probed experimentally and investigated theoretically. Here, we report on the equilibrium populations of the isomers in the temperature range from 30K to 150K for both H2O and D2O as predicted by four different water potentials. The simulations, performed using path-integral molecular dynamics combined with the replica exchange method, highlight some deficiencies of empirical water models while providing support for the accuracy of more recent ab initio-based potentials. The theoretical predictions for the cage/prism isomeric equilibrium upon isotopic substitution suggest that rotational spectra measured for the deuterated cluster could deliver further insights on the ground-state properties of the water hexamer.

On the absolute thermodynamics of water from computer simulations: A comparison of first-principles molecular dynamics, reactive and empirical force fields

We present the absolute enthalpy, entropy, heat capacity, and free energy of liquid water at ambient conditions calculated by the two-phase thermodynamic method applied to ab initio, reactive and classical molecular dynamics simulations. We find that the absolute entropy and heat capacity of liquid water from ab initio molecular dynamics (AIMD) is underestimated, but falls within the range of the flexible empirical as well as the reactive force fields. The origin of the low absolute entropy of liquid water from AIMD simulations is due to an underestimation of the translational entropy by 20% and the rotational entropy by 40% compared to the TIP3P classical water model, consistent with previous studies that reports low diffusivity and increased ordering of liquid water from AIMD simulations. Classical MD simulations with rigid water models tend to be in better agreement with experiment (in particular TIP3P yielding the best agreement), although the TIP4P-ice water model, the only empirical force field that reproduces the experimental melting temperature, has the lowest entropy, perhaps expectedly. This reiterates the limitations of existing empirical water models in simultaneously capturing the thermodynamics of solid and liquid phases. We find that the quantum corrections to heat capacity of water can be as large as 60%. Although certain water models are computed to yield good absolute free energies of water compared to experiments, they are often due to the fortuitous enthalpy-entropy cancellation, but not necessarily due to the correct descriptions of enthalpy and entropy separately.

A Unified Framework for Symbiosis of Evolutionary Mechanisms with Application to Water Clusters Potential Model Design

This article presents a theoretic model for facilitating the emergence of productive search profiles transpiring from the symbiosis of gene (stochastic variation) and meme (lifetime learning) working in synergy. The evolvability measure of the symbiotic search profiles for each individual is quantified by means of statistical learning on distinct sample vectors encountered along the search. The most productive search profile inferred for an individual, as defined by evolvability measure, is subsequently used to work on it, leading to the self-configuration of solvers that acclimatizes to suit the given problem of interest. Empirical studies on representative problems are presented to reflect the characteristics of symbiotic evolution. Assessment made against several recent state-of-the-art evolutionary and adaptive search algorithms highlighted the efficacy of the theoretic formalism of evolutionary mechanisms in symbiosis for autonomic search. As the design of computationally cheap advanced empirical water models for the understanding of enigmatic properties of water remains an important and unsolved problem, the article presents an illustration of symbiotic evolution for the design of (H <sub xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">2</sub> O) <sub xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">n</sub> or water clusters potential model.

Accurate Predictions of Water Cluster Formation, (H2O)n=2−10

An efficient mixed molecular dynamics/quantum mechanics model has been applied to the water cluster system. The use of the MP2 method and correlation consistent basis sets, with appropriate correction for BSSE, allows for the accurate calculation of electronic and free energies for the formation of clusters of 2-10 water molecules. This approach reveals new low energy conformers for (H(2)O)(n=7,9,10). The water heptamer conformers comprise five different structural motifs ranging from a three-dimensional prism to a quasi-planar book structure. A prism-like structure is favored energetically at low temperatures, but a chair-like structure is the global Gibbs free energy minimum past 200 K. The water nonamers exhibit less complexity with all the low energy structures shaped like a prism. The decamer has 30 conformers that are within 2 kcal/mol of the Gibbs free energy minimum structure at 298 K. These structures are categorized into four conformer classes, and a pentagonal prism is the most stable structure from 0 to 320 K. Results can be used as benchmark values for empirical water models and density functionals, and the method can be applied to larger water clusters.

Compatibility of Quantum Chemical Methods and Empirical (MM) Water Models in Quantum Mechanics/Molecular Mechanics Liquid Water Simulations

High-level (MP2, HF, and BLYP with the aug-cc-pVDZ basis set) quantum mechanics/molecular mechanics (QM/MM) Monte Carlo free energy simulations of liquid water are used here to test the compatibility of various QM methods with four standard empirical “molecular mechanics” (MM) water models. Consistency of QM methods with water models is of particular importance, given the aqueous environment of many of the systems of interest for QM/MM modeling (e.g., biological systems). The results show that treating a single water molecule using a QM method in bulk TIP3P can induce solvent structuring consistent with experiment. The results also show that the TIP4P model is the most suitable water model of those tested for such QM/MM simulations, while the TIP5P model is not well suited. The findings have important implications for future QM/MM method development and applications. They indicate that the choice of MM models should be made carefully for consistency and compatibility in QM/MM simulations.

Competing quantum effects in the dynamics of a flexible water model

Numerous studies have identified large quantum mechanical effects in the dynamics of liquid water. In this paper, we suggest that these effects may have been overestimated due to the use of rigid water models and flexible models in which the intramolecular interactions were described using simple harmonic functions. To demonstrate this, we introduce a new simple point charge model for liquid water, q-TIP4P/F, in which the O-H stretches are described by Morse-type functions. We have parametrized this model to give the correct liquid structure, diffusion coefficient, and infrared absorption frequencies in quantum (path integral-based) simulations. The model also reproduces the experimental temperature variation of the liquid density and affords reasonable agreement with the experimental melting temperature of hexagonal ice at atmospheric pressure. By comparing classical and quantum simulations of the liquid, we find that quantum mechanical fluctuations increase the rates of translational diffusion and orientational relaxation in our model by a factor of around 1.15. This effect is much smaller than that observed in all previous simulations of empirical water models, which have found a quantum effect of at least 1.4 regardless of the quantum simulation method or the water model employed. The small quantum effect in our model is a result of two competing phenomena. Intermolecular zero point energy and tunneling effects destabilize the hydrogen-bonding network, leading to a less viscous liquid with a larger diffusion coefficient. However, this is offset by intramolecular zero point motion, which changes the average water monomer geometry resulting in a larger dipole moment, stronger intermolecular interactions, and a slower diffusion. We end by suggesting, on the basis of simulations of other potential energy models, that the small quantum effect we find in the diffusion coefficient is associated with the ability of our model to produce a single broad O-H stretching band in the infrared absorption spectrum.

Structure of Aqueous Proline via Parallel Tempering Molecular Dynamics and Neutron Diffraction

The structure of aqueous L-proline amino acid has been the subject of much debate centering on the validity of various proposed models, differing widely in the extent to which local and long-range correlations are present. Here, aqueous proline is investigated by atomistic, replica exchange molecular dynamics simulations, and the results are compared to neutron diffraction and small angle neutron scattering (SANS) data, which have been reported recently (McLain, S.; Soper, A.; Terry, A.; Watts, A. J. Phys. Chem. B 2007, 111, 4568). Comparisons between neutron experiments and simulation are made via the static structure factor S(Q) which is measured and computed from several systems with different H/D isotopic compositions at a concentration of 1:20 molar ratio. Several different empirical water models (TIP3P, TIP4P, and SPC/E) in conjunction with the CHARMM22 force field are investigated. Agreement between experiment and simulation is reasonably good across the entire Q range although there are significant model-dependent variations in some cases. In general, agreement is improved slightly upon application of approximate quantum corrections obtained from gas-phase path integral simulations. Dimers and short oligomeric chains formed by hydrogen bonds (frequently bifurcated) coexist with apolar (hydrophobic) contacts. These emerge as the dominant local motifs in the mixture. Evidence for long-range association is more equivocal: No long-range structures form spontaneously in the MD simulations, and no obvious low-Q signature is seen in the SANS data. Moreover, associations introduced artificially to replicate a long-standing proposed mesoscale structure for proline correlations as an initial condition are annealed out by parallel tempering MD simulations. However, some small residual aggregates do remain, implying a greater degree of long-range order than is apparent in the SANS data.

“Tetrahedrality” and the Relationship between Collective Structure and Radial Distribution Functions in Liquid Water

Molecular dynamics simulations of pure liquid water under ambient conditions using four common empirical water models have been analyzed to determine how well the oxygen-oxygen radial distribution function, g(OO)(r), used as the sole criterion of congruence with experiment, captures variations in the actual anisotropic collective structuring for these models. The largest systematic deviations from tetrahedrality were found to be due to deformations of the angle between the two closet hydrogen bond donor neighbors, but for intrinsic geometric reasons, these were found to contribute less to g(OO)(r) than deformations of the angles between one hydrogen bond donor neighbor and one hydrogen bond acceptor neighbor. Relying exclusively on a qualitative characterization of the second peak in g(OO)(r) seems to overemphasize the differences between the structuring in some of these models.

DIRECT INTERACTION ENERGY: A COMPUTATIONAL QUANTITY FOR PARAMETERIZATION OF CONDENSED-PHASE FORCE FIELDS AND ITS APPLICATION TO HYDROGEN BONDING

Using N-Methylacetamide (NMA) dimer and NMA–water as model complexes, the solvent effect on the protein inter- N – H ⋯ O =C and intra- N – H ⋯ OH 2, and C = O ⋯ H 2 O hydrogen bonding have been studied by the polarizable continuum model (PCM) ab initio calculations in the four media (vacuum, ether, nitromethane and water). In contrast to the empirical approaches, we suggested using the direction interaction energies (DE) to consider the solvent polarization, which can be derived from PCM ab initio calculations. The DEs of the model compounds in solvents are larger than their in vacuo binding energies, which reflect the solvent polarization effect. As the solvents become increasingly polar, the binding free energies decrease while DEs increase. The increasing DE is consistent with the increasing hydrogen bond length. Considering the protein environment, the DEs of NMA-NMA dimer in ether, 9.14 and 9.41 kcal/mol for NMADI and NMADII, are recommended for the intra N – H ⋯ O =C hydrogen bonding. The DEs of NMA–water complex in water, -5.47 (NMAWI) and -5.41 kcal/mol (NMAWI'), -8.44 (NMAWII) and -8.68 kcal/mol (NMAWII'), respectively, are suggested for the inter- N – H ⋯ OH 2 and C – O ⋯ H 2 O hydrogen bonding of proteins. Using the same approach, we have also computed the DE of water dimer in liquid water. The computed DE of water dimer (-5.63 kcal/mol) is larger than the in vacuo water dimerization energy (-5.14 kcal/mol) and in reasonable agreement with the dimerization energies (ranging from –6.0 to 6.8 kcal/mol) of polarization-included empirical water models.

Ab initio x-ray scattering of liquid water

A direct calculation of the coherent x-ray scattering spectrum of liquid water under ambient conditions is presented using ab initio density functional theory. The experimental data are compared with the calculated x-ray scattering spectra retrieved from the trajectories of two Car–Parrinello molecular dynamics runs of about 10 ps with 32 and 64 water molecules in the simulation cell, respectively. Furthermore, the x-ray spectra obtained from molecular dynamics runs of 20 ps for 128 water molecules using empirical water models are presented. All calculations were performed with the QUICKSTEP code, which includes an all-electron implementation of the Gaussian augmented plane wave method. It is shown that by this approach one can assess how well different models describe the structural properties of water.

A Molecular Dynamics Study of Tryptophan in Water

Molecular dynamics simulations of a single tryptophan amino acid in water were performed with the CHARMM molecular mechanics program using the modified TIP3P, the refined SPC, and the original SPC/E water models. All these empirical water models are similar in nature, but small differences give significant differences in their properties for liquid water. Nanosecond molecular dynamics simulations were carried out in the NVT, NVE, or NPT ensembles using a cubic simulation cell furnished with periodic boundary conditions. The calculations of long-range interactions were performed with two different methods, the atom-based spherical cutoff method with force shifting, where the nonbonded interactions were smoothly shifted to zero at the cutoff distance, and the particle-mesh Ewald technique. Transition dipole autocorrelation functions and reorientation times of tryptophan were calculated and compared to experimental values. Tryptophan hydration was affected by the water model and the treatment of long-range electrostatic interactions. The reorientation time of tryptophan depended strongly on the water model used, and the viscosity of the model liquid was found to correlate with the tryptophan translational and rotational dynamics in water. Statistical accuracy and errors in the analyses are discussed.

Structure and Dynamics of the TIP3P, SPC, and SPC/E Water Models at 298 K

Molecular dynamics simulations of five water models, the TIP3P (original and modified), SPC (original and refined), and SPC/E (original), were performed using the CHARMM molecular mechanics program. All simulations were carried out in the microcanonical NVE ensemble, using 901 water molecules in a cubic simulation cell furnished with periodic boundary conditions at 298 K. The SHAKE algorithm was used to keep water molecules rigid. Nanosecond trajectories were calculated with all water models for high statistical accuracy. The characteristic self-diffusion coefficients D and radial distribution functions, gOO, gOH, and gHH for all five water models were determined and compared to experimental data. The effects of velocity rescaling on the self-diffusion coefficient D were examined. All these empirical water models used in this study are similar by having three interaction sites, but the small differences in their pair potentials composed of Lennard-Jones (LJ) and Coulombic terms give significant differences in the calculated self-diffusion coefficients, and in the height of the second peak of the radial distribution function gOO.