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.

The Microscopic Physical Cause for the Density Maximum of Liquid Water.

The existence of a density maximum at 277 K is probably the most prominent anomaly among the many very special thermodynamic properties of liquid water. While usually attributed to so-called hydrogen bonding, the microscopic physical cause of this prominent anomaly is still elusive. Here we show that the density anomaly is caused by those short-range electrostatic forces, which are generated by the quadrupole and higher moments of the charge distributions present in liquid-phase water molecules. This conclusion derives from 20 ns replica exchange molecular-dynamics simulations with closely related polarizable four-, five-, and six-point water models. As soon as the model complexity suffices to represent the higher electrostatic moments with sufficient accuracy, the density temperature profile n(T) calculated for T ∈ [250,320] K at the standard pressure 1 bar locks in to the experimental observation. The corresponding six-point model is, therefore, the most simple available cartoon for liquid-phase water molecules.

Polarizable water models from mixed computational and empirical optimization.

Here we suggest a mixed computational and empirical approach serving to optimize the parameters of complex and polarizable molecular mechanics (PMM) models for complicated liquids. The computational part of the parameter optimization relies on hybrid calculations combining density functional theory (DFT) for a solute molecule with a PMM treatment of its solvent environment at well-defined thermodynamic conditions. As an application we have developed PMM models for water featuring ν = 3, 4, and 5 points of force action, a Gaussian inducible dipole and a Buckingham potential at the oxygen, the experimental liquid phase geometry, the experimental gas phase polarizability α(exp)(g) = 1.47 ų, and, for ν = 4 and 5, the gas phase value µ(exp)(g) = 1.855 D for the static dipole moment. The widths of the Gaussian dipoles and, for ν = 4 and 5, also the electrostatic geometries of these so-called TLνP models are derived from self-consistent DFT/PMM calculations, and the parameters of the Buckingham potentials (and the static TL3P dipole moment) are estimated from molecular dynamics (MD) simulations. The high quality of the resulting models is demonstrated for the observables targeted during optimization (potential energy per molecule, pressure, radial distribution functions) and a series of predicted properties (quadrupole moments, density at constant pressure, dielectric constant, diffusivity, viscosity, compressibility, heat capacity) at certain standard conditions. Remaining deficiencies and possible ways for their removal are discussed.

Coupling density functional theory to polarizable force fields for efficient and accurate Hamiltonian molecular dynamics simulations.

Hybrid molecular dynamics (MD) simulations, in which the forces acting on the atoms are calculated by grid-based density functional theory (DFT) for a solute molecule and by a polarizable molecular mechanics (PMM) force field for a large solvent environment composed of several 10(3)-10(5) molecules, pose a challenge. A corresponding computational approach should guarantee energy conservation, exclude artificial distortions of the electron density at the interface between the DFT and PMM fragments, and should treat the long-range electrostatic interactions within the hybrid simulation system in a linearly scaling fashion. Here we describe a corresponding Hamiltonian DFT/(P)MM implementation, which accounts for inducible atomic dipoles of a PMM environment in a joint DFT/PMM self-consistency iteration. The long-range parts of the electrostatics are treated by hierarchically nested fast multipole expansions up to a maximum distance dictated by the minimum image convention of toroidal boundary conditions and, beyond that distance, by a reaction field approach such that the computation scales linearly with the number of PMM atoms. Short-range over-polarization artifacts are excluded by using Gaussian inducible dipoles throughout the system and Gaussian partial charges in the PMM region close to the DFT fragment. The Hamiltonian character, the stability, and efficiency of the implementation are investigated by hybrid DFT/PMM-MD simulations treating one molecule of the water dimer and of bulk water by DFT and the respective remainder by PMM.

Vibrational spectra of phosphate ions in aqueous solution probed by first-principles molecular dynamics.

We have carried out "first-principles" Born-Oppenheimer molecular dynamics (BOMD) simulations of the phosphate ions H2PO4⁻ and HPO4²â» in liquid water and have calculated their IR spectra by Fourier transform techniques from the trajectories. IR bands were assigned by a so-called "generalized normal coordinate analysis". The effects of including Hartree-Fock (HF) exchange into the density functional theory (DFT) computation of forces were studied by comparing results obtained with the well-known BP, BLYP, and B3LYP functionals. The neglect of dispersion in the functionals was empirically corrected. The inclusion of HF exchange turned out to yield dramatically improved and, thus, quite accurate descriptions of the IR spectra observed for H2PO4⁻ and HPO4²â» in aqueous solution. An analysis of earlier computational results (Klähn, M. et al. J. Phys. Chem. A 2004, 108, 6186-6194) on these vibrational spectra, which had been obtained in a hybrid setting combining a BP description of the respective phosphate with a simple molecular mechanics (MM) model of its aqueous environment, revealed three different sources of error, (i) the BP force field of the phosphates is much too soft and would have required a substantial scaling of frequencies, (ii) the oversimplified water force field entailed incorrect solvation structures and, thus, qualitatively wrong patterns of solvatochromic band shifts, and (iii) quantitative frequency computations additionally required the inclusion of HF exchange. Thus, the results of the B3LYP BOMD simulations do not only characterize physical properties like the IR spectra or the solvation structures of the phosphate systems but also provide clues for the future design of simplified but nevertheless reasonably accurate DFT/MM methods applicable to phosphates.

Optimizing the Accuracy and Efficiency of Fast Hierarchical Multipole Expansions for MD Simulations.

Based on p'th order Cartesian Taylor expansions of Coulomb interactions and on hierarchical decompositions of macromolecular simulation systems into hierarchies of nested, structure-adapted, and adaptively formed clusters of increasing size, fast multipole methods are constructed for rapid and accurate calculations of electrostatic interactions. These so-called SAMMp algorithms are formulated through totally symmetric and traceless tensors describing the multipole moments and the coefficients of local Taylor expansions. Simple recursions for the efficient evaluation and shifting of multipole moments are given. The required tensors are explicitly given up to order p = 4. The SAMMp algorithms are shown to guarantee the reaction principle. For systems with periodic boundaries, a reaction field (RF) correction is applied, which introduces at distances beyond the "minimum image convention" boundary a dielectric continuum surrounding each cluster at the top level of coarse graining. The correctness of the present SAMMp implementation is demonstrated by analyzing the scaling of the residuals and by checking the numerical accuracy of the reaction principle for a pair of distant molecular ions in vacuum. Molecular dynamics simulations of pure water and aqueous solutions containing artificial ions, which are enclosed by periodic boundaries, demonstrate the stability and low-noise behavior of SAMMp/RF.