Systematic improvement of a classical molecular model of water

Thermophysical properties of supercritical water and bond flexibility

Molecular dynamics results are reported for the thermodynamic properties of supercritical water using examples of both rigid (TIP4P/2005) and flexible (TIP4P/2005f) transferable interaction potentials. Data are reported for pressure, isochoric and isobaric heat capacities, the thermal expansion coefficient, isothermal and adiabatic compressibilities, Joule-Thomson coefficient, speed of sound, self-diffusion coefficient, viscosities, and thermal conductivity. Many of these properties have unusual behavior in the supercritical phase such as maximum and minimum values. The effectiveness of bond flexibility on predicting these properties is determined by comparing the results to experimental data. The influence of the intermolecular potential on these properties is both variable and state point dependent. In the vicinity of the critical density, the rigid and flexible potentials yield very different values for the compressibilities, heat capacities, and thermal expansion coefficient, whereas the self-diffusion coefficient, viscosities, and thermal conductivities are much less potential dependent. Although the introduction of bond flexibility is a computationally expedient way to improve the accuracy of an intermolecular potential, it can be counterproductive in some cases and it is not an adequate replacement for incorporating the effects of polarization.

2D-Raman-THz spectroscopy: a sensitive test of polarizable water models

In a recent paper, the experimental 2D-Raman-THz response of liquid water at ambient conditions has been presented [J. Savolainen, S. Ahmed, and P. Hamm, Proc. Natl. Acad. Sci. U. S. A. 110, 20402 (2013)]. Here, all-atom molecular dynamics simulations are performed with the goal to reproduce the experimental results. To that end, the molecular response functions are calculated in a first step, and are then convoluted with the laser pulses in order to enable a direct comparison with the experimental results. The molecular dynamics simulation are performed with several different water models: TIP4P/2005, SWM4-NDP, and TL4P. As polarizability is essential to describe the 2D-Raman-THz response, the TIP4P/2005 water molecules are amended with either an isotropic or a anisotropic polarizability a posteriori after the molecular dynamics simulation. In contrast, SWM4-NDP and TL4P are intrinsically polarizable, and hence the 2D-Raman-THz response can be calculated in a self-consistent way, using the same force field as during the molecular dynamics simulation. It is found that the 2D-Raman-THz response depends extremely sensitively on details of the water model, and in particular on details of the description of polarizability. Despite the limited time resolution of the experiment, it could easily distinguish between various water models. Albeit not perfect, the overall best agreement with the experimental data is obtained for the TL4P water model.

Assessing the potential of atomistic molecular dynamics simulations to probe reversible protein-protein recognition and binding

Protein-protein recognition and binding are governed by diffusion, noncovalent forces and conformational flexibility, entangled in a way that only molecular dynamics simulations can dissect at high resolution. Here we exploited ubiquitin's noncovalent dimerization equilibrium to assess the potential of atomistic simulations to reproduce reversible protein-protein binding, by running submicrosecond simulations of systems with multiple copies of the protein at millimolar concentrations. The simulations essentially fail because they lead to aggregates, yet they reproduce some specificity in the binding interfaces as observed in known covalent and noncovalent ubiquitin dimers. Following similar observations in literature we hint at electrostatics and water descriptions as the main liable force field elements, and propose that their optimization should consider observables relevant to multi-protein systems and unfolded proteins. Within limitations, analysis of binding events suggests salient features of protein-protein recognition and binding, to be retested with improved force fields. Among them, that specific configurations of relative direction and orientation seem to trigger fast binding of two molecules, even over 50 Å distances; that conformational selection can take place within surface-to-surface distances of 10 to 40 Å i.e. well before actual intermolecular contact; and that establishment of contacts between molecules further locks their conformations and relative orientations.

Numerical study on the partitioning of the molecular polarizability into fluctuating charge and induced atomic dipole contributions

In order to carry out a detailed analysis of the molecular static polarizability, which is the response of the molecule to a uniform external electric field, the molecular polarizability was computed using the finite-difference method for 21 small molecules, using density functional theory. Within nine charge population schemes (Löwdin, Mulliken, Becke, Hirshfeld, CM5, Hirshfeld-I, NPA, CHELPG, MK-ESP) in common use, the charge fluctuation contribution is found to dominate the molecular polarizability, with its ratio ranging from 59.9% with the Hirshfeld or CM5 scheme to 96.2% with the Mulliken scheme. The Hirshfeld-I scheme is also used to compute the other contribution to the molecular polarizability coming from the induced atomic dipoles, and the atomic polarizabilities in eight small molecules and water pentamer are found to be highly anisotropic for most atoms. Overall, the results suggest that (a) more emphasis probably should be placed on the charge fluctuation terms in future polarizable force field development and (b) an anisotropic polarizability might be more suitable than an isotropic one in polarizable force fields based entirely or partially on the induced atomic dipoles.

An efficient algorithm for multipole energies and derivatives based on spherical harmonics and extensions to particle mesh Ewald

Next-generation molecular force fields deliver accurate descriptions of non-covalent interactions by employing more elaborate functional forms than their predecessors. Much work has been dedicated to improving the description of the electrostatic potential (ESP) generated by these force fields. A common approach to improving the ESP is by augmenting the point charges on each center with higher-order multipole moments. The resulting anisotropy greatly improves the directionality of the non-covalent bonding, with a concomitant increase in computational cost. In this work, we develop an efficient strategy for enumerating multipole interactions, by casting an efficient spherical harmonic based approach within a particle mesh Ewald (PME) framework. Although the derivation involves lengthy algebra, the final expressions are relatively compact, yielding an approach that can efficiently handle both finite and periodic systems without imposing any approximations beyond PME. Forces and torques are readily obtained, making our method well suited to modern molecular dynamics simulations.

On combining Thole's induced point dipole model with fixed charge distributions in molecular mechanics force fields

The Thole induced point dipole model is combined with three different point charge fitting methods, Merz-Kollman (MK), charges from electrostatic potentials using a grid (CHELPG), and restrained electrostatic potential (RESP), and two multipole algorithms, distributed multipole analysis (DMA) and Gaussian multipole model (GMM), which can be used to describe the electrostatic potential (ESP) around molecules in molecular mechanics force fields. This is done to study how the different methods perform when intramolecular polarizability contributions are self-consistently removed from the fitting done in the force field parametrization. It is demonstrated that the polarizable versions of the partial charge models provide a good compromise between accuracy and computational efficiency in describing the ESP of small organic molecules undergoing conformational changes. For the point charge models, the inclusion of polarizability reduced the the average root mean square error of ESP over the test set by 4-10%.

Water dispersion interactions strongly influence simulated structural properties of disordered protein states

Many proteins can be partially or completely disordered under physiological conditions. Structural characterization of these disordered states using experimental methods can be challenging, since they are composed of a structurally heterogeneous ensemble of conformations rather than a single dominant conformation. Molecular dynamics (MD) simulations should in principle provide an ideal tool for elucidating the composition and behavior of disordered states at an atomic level of detail. Unfortunately, MD simulations using current physics-based models tend to produce disordered-state ensembles that are structurally too compact relative to experiments. We find that the water models typically used in MD simulations significantly underestimate London dispersion interactions, and speculate that this may be a possible reason for these erroneous results. To test this hypothesis, we create a new water model, TIP4P-D, that approximately corrects for these deficiencies in modeling water dispersion interactions while maintaining compatibility with existing physics-based models. We show that simulations of solvated proteins using this new water model typically result in disordered states that are substantially more expanded and in better agreement with experiment. These results represent a significant step toward extending the range of applicability of MD simulations to include the study of (partially or fully) disordered protein states.

Ab initio water pair potential with flexible monomers

A potential energy surface for the water dimer with explicit dependence on monomer coordinates is presented. The surface was fitted to a set of previously published interaction energies computed on a grid of over a quarter million points in the 12-dimensional configurational space using symmetry-adapted perturbation theory and coupled-cluster methods. The present fit removes small errors in published fits, and its accuracy is critically evaluated. The minimum and saddle-point structures of the potential surface were found to be very close to predictions from direct ab initio optimizations. The computed second virial coefficients agreed well with experimental values. At low temperatures, the effects of monomer flexibility in the virial coefficients were found to be much smaller than the quantum effects.

Why many semiempirical molecular orbital theories fail for liquid water and how to fix them

Water is an extremely important liquid for chemistry and the search for more accurate force fields for liquid water continues unabated. Neglect of diatomic differential overlap (NDDO) molecular orbital methods provide and intriguing generalization of classical force fields in this regard because they can account both for bond breaking and electronic polarization of molecules. However, we show that most standard NDDO methods fail for water because they give an incorrect description of hydrogen bonding, water's key structural feature. Using force matching, we design a reparameterized NDDO model and find that it qualitatively reproduces the experimental radial distribution function of water, as well as various monomer, dimer, and bulk properties that PM6 does not. This suggests that the apparent limitations of NDDO models are primarily due to poor parameterization and not to the NDDO approximations themselves. Finally, we identify the physical parameters that most influence the condensed phase properties. These results help to elucidate the chemistry that a semiempirical molecular orbital picture of water must capture. We conclude that properly parameterized NDDO models could be useful for simulations that require electronically detailed explicit solvent, including the calculation of redox potentials and simulation of charge transfer and photochemistry.

Revised Parameters for the AMOEBA Polarizable Atomic Multipole Water Model

A set of improved parameters for the AMOEBA polarizable atomic multipole water model is developed. An automated procedure, ForceBalance, is used to adjust model parameters to enforce agreement with ab initio-derived results for water clusters and experimental data for a variety of liquid phase properties across a broad temperature range. The values reported here for the new AMOEBA14 water model represent a substantial improvement over the previous AMOEBA03 model. The AMOEBA14 model accurately predicts the temperature of maximum density and qualitatively matches the experimental density curve across temperatures from 249 to 373 K. Excellent agreement is observed for the AMOEBA14 model in comparison to experimental properties as a function of temperature, including the second virial coefficient, enthalpy of vaporization, isothermal compressibility, thermal expansion coefficient, and dielectric constant. The viscosity, self-diffusion constant, and surface tension are also well reproduced. In comparison to high-level ab initio results for clusters of 2-20 water molecules, the AMOEBA14 model yields results similar to AMOEBA03 and the direct polarization iAMOEBA models. With advances in computing power, calibration data, and optimization techniques, we recommend the use of the AMOEBA14 water model for future studies employing a polarizable water model.

Combined covalent-electrostatic model of hydrogen bonding improves structure prediction with Rosetta

Interactions between polar atoms are challenging to model because at very short ranges they form hydrogen bonds (H-bonds) that are partially covalent in character and exhibit strong orientation preferences; at longer ranges the orientation preferences are lost, but significant electrostatic interactions between charged and partially charged atoms remain. To simultaneously model these two types of behavior, we refined an orientation dependent model of hydrogen bonds [Kortemme et al. J. Mol. Biol. 2003, 326, 1239] used by the molecular modeling program Rosetta and then combined it with a distance-dependent Coulomb model of electrostatics. The functional form of the H-bond potential is physically motivated and parameters are fit so that H-bond geometries that Rosetta generates closely resemble H-bond geometries in high-resolution crystal structures. The combined potentials improve performance in a variety of scientific benchmarks including decoy discrimination, side chain prediction, and native sequence recovery in protein design simulations and establishes a new standard energy function for Rosetta.

Systematic procedure to parametrize force fields for molecular fluids

A new strategy to develop force fields for molecular fluids is presented. The intermolecular parameters are fitted to reproduce experimental values of target properties at ambient conditions and also the critical temperature. The partial charges are chosen to match the dielectric constant. The Lennard-Jones parameters, εii and σii, are fitted to reproduce the surface tension at the vapor-liquid interface and the liquid density, respectively. The choice of those properties allows obtaining systematically the final parameters using a small number of simulations. It is shown that the use of surface tension as a target property is better than the choice of heat of vaporization. The method is applied to molecules, from all atoms to a coarse-grained level, such as pyridine, dichloromethane, methanol, and 1-ethyl-3-methylimidazolium tetrafluoroborate (EMIM-BF4) at different temperatures and pressures. The heat of vaporization, radial distribution functions, and self-diffusion coeficient are also calculated.

Bayesian energy landscape tilting: towards concordant models of molecular ensembles

Predicting biological structure has remained challenging for systems such as disordered proteins that take on myriad conformations. Hybrid simulation/experiment strategies have been undermined by difficulties in evaluating errors from computational model inaccuracies and data uncertainties. Building on recent proposals from maximum entropy theory and nonequilibrium thermodynamics, we address these issues through a Bayesian energy landscape tilting (BELT) scheme for computing Bayesian hyperensembles over conformational ensembles. BELT uses Markov chain Monte Carlo to directly sample maximum-entropy conformational ensembles consistent with a set of input experimental observables. To test this framework, we apply BELT to model trialanine, starting from disagreeing simulations with the force fields ff96, ff99, ff99sbnmr-ildn, CHARMM27, and OPLS-AA. BELT incorporation of limited chemical shift and (3)J measurements gives convergent values of the peptide's α, β, and PPII conformational populations in all cases. As a test of predictive power, all five BELT hyperensembles recover set-aside measurements not used in the fitting and report accurate errors, even when starting from highly inaccurate simulations. BELT's principled framework thus enables practical predictions for complex biomolecular systems from discordant simulations and sparse data.

Building Water Models: A Different Approach

Simplified classical water models are currently an indispensable component in practical atomistic simulations. Yet, despite several decades of intense research, these models are still far from perfect. Presented here is an alternative approach to constructing widely used point charge water models. In contrast to the conventional approach, we do not impose any geometry constraints on the model other than the symmetry. Instead, we optimize the distribution of point charges to best describe the "electrostatics" of the water molecule. The resulting "optimal" 3-charge, 4-point rigid water model (OPC) reproduces a comprehensive set of bulk properties significantly more accurately than commonly used rigid models: average error relative to experiment is 0.76%. Close agreement with experiment holds over a wide range of temperatures. The improvements in the proposed model extend beyond bulk properties: compared to common rigid models, predicted hydration free energies of small molecules using OPC are uniformly closer to experiment, with root-mean-square error <1 kcal/mol.

Predictions for water clusters from a first-principles two- and three-body force field

A new rigid-monomer three-body potential has been developed for water by fitting it to more than 70 thousand trimer interaction energies computed ab initio using coupled-cluster methods and augmented triple-zeta-quality basis sets. This potential was used together with a modified form of a previously developed two-body potential and with a polarization model of four- and higher-body interactions to predict the energetics of the water trimer, hexamer, and 24-mer. Despite using the rigid-monomer approximation, these predictions agree better with flexible-monomer benchmarks than published results obtained with flexible-monomer force fields. An unexpected finding of our work is that simple polarization models predict four-body interactions to within a few percent, whereas for three-body interactions these models are known to have errors on the order of 50%.

Charge Anisotropy: Where Atomic Multipoles Matter Most

Specific intermolecular interactions are largely guided by electrostatics. However, the most common model for electrostatic interactions-atomic point charges-fails to reproduce anisotropic charge distributions, such as lone pairs and sigma holes. Although this has long been known, point charges are still widely used in chemical modeling and reasoning. In this contribution, we analyze the deficiency of atomic point charges in reproducing the electrostatic potential (ESP) around molecules and find that multipole moments up to quadrupoles can, with a much lower error than point charges, reproduce the relevant ESP for all cases. Mapping the surface to the closest atom allows to compare ESP errors between atom types and to identify cases with the most urgent need for atomic multipoles. Our analysis shows that almost all heteroatoms require multipoles to correctly describe their charge distribution, with the most serious cases being nitrogen, sulfur, and halogens. Comparison with small molecule crystallography data studies supports our findings and emphasizes the need for incorporating anisotropic charge descriptions in chemical models. The scheme introduced here can be used to identify anisotropic binding preferences for atom types where there is too little coverage in crystal structure databases.

Recent Advances in Polarizable Force Fields for Macromolecules: Microsecond Simulations of Proteins Using the Classical Drude Oscillator Model

In this Perspective, we summarize recent efforts to include the explicit treatment of induced electronic polarization in biomolecular force fields. Methods used to treat polarizability, including the induced dipole, fluctuating charge, and classical Drude oscillator models, are presented, including recent advances in force fields using those methods. This is followed by recent results obtained with the Drude model, including microsecond molecular dynamics (MD) simulations of multiple proteins in explicit solvent. Results show significant variability of backbone and side-chain dipole moments as a function of environment, including significant changes during individual simulations. Dipole moments of water in the vicinity of the proteins reveal small but systematic changes, with the direction of the changes dependent on the environment. Analyses of the full proteins show that the polarizable Drude model leads to larger values of the dielectric constant of the protein interior, especially in the case of hydrophobic regions. These results indicate that the inclusion of explicit electronic polarizability leads to significant differences in the physical forces affecting the structure and dynamics of proteins, which can be investigated in a computationally tractable fashion in the context of the Drude model.

Molecular simulation of water and hydration effects in different environments: challenges and developments for DFTB based models

We discuss the description of water and hydration effects that employs an approximate density functional theory, DFTB3, in either a full QM or QM/MM framework. The goal is to explore, with the current formulation of DFTB3, the performance of this method for treating water in different chemical environments, the magnitude and nature of changes required to improve its performance, and factors that dictate its applicability to reactions in the condensed phase in a QM/MM framework. A relatively minor change (on the scale of kBT) in the O-H repulsive potential is observed to substantially improve the structural properties of bulk water under ambient conditions; modest improvements are also seen in dynamic properties of bulk water. This simple change also improves the description of protonated water clusters, a solvated proton, and to a more limited degree, a solvated hydroxide. By comparing results from DFTB3 models that differ in the description of water, we confirm that proton transfer energetics are adequately described by the standard DFTB3/3OB model for meaningful mechanistic analyses. For QM/MM applications, a robust parametrization of QM-MM interactions requires an explicit consideration of condensed phase properties, for which an efficient sampling technique was developed recently and is reviewed here. The discussions help make clear the value and limitations of DFTB3 based simulations, as well as the developments needed to further improve the accuracy and transferability of the methodology.

Some practical approaches to treating electrostatic polarization of proteins

Conspectus Electrostatic interaction plays a significant role in determining many properties of biomolecules, which exist and function in aqueous solution, a highly polar environment. For example, proteins are composed of amino acids with charged, polar, and nonpolar side chains and their specific electrostatic properties are fundamental to the structure and function of proteins. An important issue that arises in computational study of biomolecular interaction and dynamics based on classical force field is lack of polarization. Polarization is a phenomenon in which the charge distribution of an isolated molecule will be distorted when interacting with another molecule or presented in an external electric field. The distortion of charge distribution is intended to lower the overall energy of the molecular system, which is counter balanced by the increased internal energy of individual molecules due to the distorted charge distributions. The amount of the charge redistribution, which characterizes the polarizability of a molecule, is determined by the level of the charge distortion. Polarization is inherently quantum mechanical, and therefore classical force fields with fixed atomic charges are incapable of capturing this important effect. As a result, simulation studies based on popular force fields, AMBER, CHARMM, etc., lack the polarization effect, which is a widely known deficiency in most computational studies of biomolecules today. Many efforts have been devoted to remedy this deficiency, such as adding additional movable charge on the atom, allowing atomic charges to fluctuate, or including induced multipoles. Although various successes have been achieved and progress at various levels has been reported over the past decades, the issue of lacking polarization in force field based simulations is far from over. For example, some of these methods do not always give converged results, and other methods require huge computational cost. This Account reviews recent work on developing polarized and polarizable force fields based on fragment quantum mechanical calculations for proteins. The methods described here are based on quantum mechanical calculations of proteins in solution, but with a different level of rigor and different computational efficiency for the molecular dynamics applications. In the general approach, a fragment quantum mechanical calculation for protein with implicit solvation is carried out to derive a polarized protein-specific charge (PPC) for any given protein structure. The PPC correctly reflects the polarization state of the protein in a given conformation, and it can also be dynamically changed as the protein changes conformation in dynamics simulations. Another approach that is computationally more efficient is the effective polarizable bond method in which only polar bonds or groups can be polarized and their polarizabilities are predetermined from quantum mechanical calculations of these groups in external electric fields. Both methods can be employed for applications in various situations by taking advantage of their unique features.

Perspective: Markov models for long-timescale biomolecular dynamics

Molecular dynamics simulations have the potential to provide atomic-level detail and insight to important questions in chemical physics that cannot be observed in typical experiments. However, simply generating a long trajectory is insufficient, as researchers must be able to transform the data in a simulation trajectory into specific scientific insights. Although this analysis step has often been taken for granted, it deserves further attention as large-scale simulations become increasingly routine. In this perspective, we discuss the application of Markov models to the analysis of large-scale biomolecular simulations. We draw attention to recent improvements in the construction of these models as well as several important open issues. In addition, we highlight recent theoretical advances that pave the way for a new generation of models of molecular kinetics.

Optimizing potentials for a liquid mixture: a new force field for a tert-butanol and water solution

A technology for optimization of potential parameters from condensed-phase simulations (POP) is discussed and illustrated. It is based on direct calculations of the derivatives of macroscopic observables with respect to the potential parameters. The derivatives are used in a local minimization scheme, comparing simulated and experimental data. In particular, we show that the Newton trust region protocol allows for more accurate and robust optimization. We apply the newly developed technology to study the liquid mixture of tert-butanol and water. We are able to obtain, after four iterations, the correct phase behavior and accurately predict the value of the Kirkwood Buff (KB) integrals. We further illustrate that a potential that is determined solely by KB information, or the pair correlation function, is not necessarily unique.

Efficient and Minimal Method to Bias Molecular Simulations with Experimental Data

A primary goal in molecular simulations is to modify the potential energy of a system so that properties of the simulation match experimental data. This is traditionally done through iterative cycles of simulation and reparameterization. An alternative approach is to bias the potential energy so that the system matches experimental data. This can be done while minimally changing the underlying free energy of the molecular simulation. Current minimal biasing methods require replicas, which can lead to unphysical dynamics and introduces new complexity: the choice of replica number and their properties. Here, we describe a new method, called experiment directed simulation that does not require replicas, converges rapidly, can match many data simultaneously, and minimally modifies the potential. The experiment directed simulation method is demonstrated on model systems and a three-component electrolyte simulation. The theory used to derive the method also provides insight into how changing a molecular force-field impacts the expected value of observables in simulation.

Building Force Fields: An Automatic, Systematic, and Reproducible Approach

The development of accurate molecular mechanics force fields is a significant challenge that must be addressed for the continued success of molecular simulation. We developed the ForceBalance method to automatically derive accurate force field parameters using flexible combinations of experimental and theoretical reference data. The method is demonstrated in the parametrization of two rigid water models, yielding new parameter sets (TIP3P-FB and TIP4P-FB) that accurately describe many physical properties of water.

A fixed-charge model for alcohol polarization in the condensed phase, and its role in small molecule hydration

We present a simple optimization strategy for incorporating experimental dielectric response information on neat liquids in classical molecular models of alcohol. Using this strategy, we determine simple and transferable hydroxyl modulation rules that, when applied to an existing molecular parameter set, result in a newly dielectric corrected (DC) parameter set. We applied these rules to the general Amber force field (GAFF) to form an initial set of GAFF-DC parameters, and we found this to lead to significant improvement in the calculated dielectric constant and hydration free energy values for a wide variety of small molecule alcohol models. Tests of the GAFF-DC parameters in the SAMPL4 blind prediction event for hydration show these changes improve agreement with experiment. Surprisingly, these simple modifications also outperform detailed quantum mechanical electric field calculations using a self-consistent reaction field environment coupling term. This work provides a potential benchmark for future developments in methods for representing condensed-phase environments in electronic structure calculations.

Modeling organochlorine compounds and the σ-hole effect using a polarizable multipole force field

The charge distribution of halogen atoms on organochlorine compounds can be highly anisotropic and even display a so-called σ-hole, which leads to strong halogen bonds with electron donors. In this paper, we have systematically investigated a series of chloromethanes with one to four chloro substituents using a polarizable multipole-based molecular mechanics model. The atomic multipoles accurately reproduced the ab initio electrostatic potential around chloromethanes, including CCl₄, which has a prominent σ-hole on the Cl atom. The van der Waals parameters for Cl were fitted to the experimental density and heat of vaporization. The calculated hydration free energy, solvent reaction fields, and interaction energies of several homo- and heterodimer of chloromethanes are in good agreement with experimental and ab initio data. This study suggests that sophisticated electrostatic models, such as polarizable atomic multipoles, are needed for accurate description of electrostatics in organochlorine compounds and halogen bonds, although further improvement is necessary for better transferability.

Polarizable six-point water models from computational and empirical optimization

Tröster et al. (J. Phys. Chem B 2013, 117, 9486-9500) recently suggested a mixed computational and empirical approach to the optimization of polarizable molecular mechanics (PMM) water models. In the empirical part the parameters of Buckingham potentials are optimized by PMM molecular dynamics (MD) simulations. The computational part applies hybrid calculations, which combine the quantum mechanical description of a H2O molecule by density functional theory (DFT) with a PMM model of its liquid phase environment generated by MD. While the static dipole moments and polarizabilities of the PMM water models are fixed at the experimental gas phase values, the DFT/PMM calculations are employed to optimize the remaining electrostatic properties. These properties cover the width of a Gaussian inducible dipole positioned at the oxygen and the locations of massless negative charge points within the molecule (the positive charges are attached to the hydrogens). The authors considered the cases of one and two negative charges rendering the PMM four- and five-point models TL4P and TL5P. Here we extend their approach to three negative charges, thus suggesting the PMM six-point model TL6P. As compared to the predecessors and to other PMM models, which also exhibit partial charges at fixed positions, TL6P turned out to predict all studied properties of liquid water at p0 = 1 bar and T0 = 300 K with a remarkable accuracy. These properties cover, for instance, the diffusion constant, viscosity, isobaric heat capacity, isothermal compressibility, dielectric constant, density, and the isobaric thermal expansion coefficient. This success concurrently provides a microscopic physical explanation of corresponding shortcomings of previous models. It uniquely assigns the failures of previous models to substantial inaccuracies in the description of the higher electrostatic multipole moments of liquid phase water molecules. Resulting favorable properties concerning the transferability to other temperatures and conditions like the melting of ice are also discussed.

Hydrogen atom in water from ambient to high temperatures

The aqueous hydrogen atom is studied with molecular dynamics simulations from ambient temperature to near the critical point. The radial distribution functions find a hydrogen atom coordination number of about 13 water molecules at 300 K to about 4 water molecules at 646 K. The radial and angular distribution functions indicate that first-shell water molecules tend to orient to maximize hydrogen bonding interactions with other water molecules. These orientational tendencies diminish with temperature. The calculated diffusion coefficient agrees very well with experimental results known near ambient temperatures. It fits a simple activation model to about 575 K, above which the diffusion becomes much faster than predicted by the fit. To temperatures of at least 500 K there is evidence for caging on a time scale of about 1 ps, but the evidence disappears at very high temperatures. Values of the aqueous hydrogen hyperfine coupling constant are obtained by averaging the results of density functional calculations on clusters extracted from the simulations. The hyperfine coupling calculations do not agree well with experiment for reasons that are not understood now, pointing to the need for further research on this problem.

Advanced potential energy surfaces for condensed phase simulation

Computational modeling at the atomistic and mesoscopic levels has undergone dramatic development in the past 10 years to meet the challenge of adequately accounting for the many-body nature of intermolecular interactions. At the heart of this challenge is the ability to identify the strengths and specific limitations of pairwise-additive interactions, to improve classical models to explicitly account for many-body effects, and consequently to enhance their ability to describe a wider range of reference data and build confidence in their predictive capacity. However, the corresponding computational cost of these advanced classical models increases significantly enough that statistical convergence of condensed phase observables becomes more difficult to achieve. Here we review a hierarchy of potential energy surface models used in molecular simulations for systems with many degrees of freedom that best meet the trade-off between accuracy and computational speed in order to define a sweet spot for a given scientific problem of interest.