A systematic development of a polarizable potential of water

Accuracy limit of non-polarizable four-point water models: TIP4P/2005 vs OPC. Should water models reproduce the experimental dielectric constant?

The last generation of four center non-polarizable models of water can be divided into two groups: those reproducing the dielectric constant of water, as OPC, and those significantly underestimating its value, as TIP4P/2005. To evaluate the global performance of OPC and TIP4P/2005, we shall follow the test proposed by Vega and Abascal in 2011 evaluating about 40 properties to fairly address this comparison. The liquid-vapor and liquid-solid equilibria are computed, as well as the heat capacities, isothermal compressibilities, surface tensions, densities of different ice polymorphs, the density maximum, equations of state at high pressures, and transport properties. General aspects of the phase diagram are considered by comparing the ratios of different temperatures (namely, the temperature of maximum density, the melting temperature of hexagonal ice, and the critical temperature). The final scores are 7.2 for TIP4P/2005 and 6.3 for OPC. The results of this work strongly suggest that we have reached the limit of what can be achieved with non-polarizable models of water and that the attempt to reproduce the experimental dielectric constant deteriorates the global performance of the water force field. The reason is that the dielectric constant depends on two surfaces (potential energy and dipole moment surfaces), whereas in the absence of an electric field, all properties can be determined simply from just one surface (the potential energy surface). The consequences of the choice of the water model in the modeling of electrolytes in water are also discussed.

Accurate Free Energies of Aqueous Electrolyte Solutions from Molecular Simulations with Non-polarizable Force Fields

Non-polarizable force fields fail to accurately predict free energies of aqueous electrolytes without compromising the predictive ability for densities and transport properties. A new approach is presented in which (1) TIP4P/2005 water and scaled charge force fields are used to describe the interactions in the liquid phase and (2) an additional Effective Charge Surface (ECS) is used to compute free energies at zero additional computational expense. The ECS is obtained using a single temperature-independent charge scaling parameter per species. Thereby, the chemical potential of water and the free energies of hydration of various aqueous salts (e.g., NaCl and LiCl) are accurately described (deviations less than 5% from experiments), in sharp contrast to calculations where the ECS is omitted (deviations larger than 20%). This approach enables accurate predictions of free energies of aqueous electrolyte solutions using non-polarizable force fields, without compromising liquid-phase properties.

Comparing machine learning potentials for water: Kernel-based regression and Behler-Parrinello neural networks

In this paper, we investigate the performance of different machine learning potentials (MLPs) in predicting key thermodynamic properties of water using RPBE + D3. Specifically, we scrutinize kernel-based regression and high-dimensional neural networks trained on a highly accurate dataset consisting of about 1500 structures, as well as a smaller dataset, about half the size, obtained using only on-the-fly learning. This study reveals that despite minor differences between the MLPs, their agreement on observables such as the diffusion constant and pair-correlation functions is excellent, especially for the large training dataset. Variations in the predicted density isobars, albeit somewhat larger, are also acceptable, particularly given the errors inherent to approximate density functional theory. Overall, this study emphasizes the relevance of the database over the fitting method. Finally, this study underscores the limitations of root mean square errors and the need for comprehensive testing, advocating the use of multiple MLPs for enhanced certainty, particularly when simulating complex thermodynamic properties that may not be fully captured by simpler tests.

Neural Network Water Model Based on the MB-Pol Many-Body Potential

The MB-pol many-body potential accurately predicts many properties of water, including cluster, liquid phase, and vapor-liquid equilibrium properties, but its high computational cost can make applying it in large-scale simulations quite challenging. In order to address this limitation, we developed a "deep potential" neural network (DPMD) model based on the MB-pol potential for water. We find that a DPMD model trained on mostly liquid configurations yields a good description of the bulk liquid phase but severely underpredicts vapor-liquid coexistence densities. By contrast, adding cluster configurations to the neural network training set leads to a good agreement for the vapor coexistence densities. Liquid phase densities under supercooled conditions are also represented well, even though they were not included in the training set. These results confirm that neural network models can combine accuracy and transferability if sufficient attention is given to the construction of a representative training set for the target system.

Computation of Electrical Conductivities of Aqueous Electrolyte Solutions: Two Surfaces, One Property

In this work, we computed electrical conductivities under ambient conditions of aqueous NaCl and KCl solutions by using the Einstein-Helfand equation. Common force fields (charge q = ±1 e) do not reproduce the experimental values of electrical conductivities, viscosities, and diffusion coefficients. Recently, we proposed the idea of using different charges to describe the potential energy surface (PES) and the dipole moment surface (DMS). In this work, we implement this concept. The equilibrium trajectories required to evaluate electrical conductivities (within linear response theory) were obtained by using scaled charges (with the value q = ±0.75 e) to describe the PES. The potential parameters were those of the Madrid-Transport force field, which accurately describe viscosities and diffusion coefficients of these ionic solutions. However, integer charges were used to compute the conductivities (thus describing the DMS). The basic idea is that although the scaled charge describes the ion-water interaction better, the integer charge reflects the value of the charge that is transported due to the electric field. The agreement obtained with experiments is excellent, as for the first time electrical conductivities (and the other transport properties) of NaCl and KCl electrolyte solutions are described with high accuracy for the whole concentration range up to their solubility limit. Finally, we propose an easy way to obtain a rough estimate of the actual electrical conductivity of the potential model under consideration using the approximate Nernst-Einstein equation, which neglects correlations between different ions.

Self-Consistent Implementation of a Solvation Free Energy Framework to Predict the Salt Solubilities of Six Alkali Halides

To assess the salt solubilities of six alkali halides in aqueous systems, we proposed a thermodynamic cycle and an efficient molecular modeling methodology. The Gibbs free energy changes for vaporization, dissociation, and dissolution were calculated using the experimental data of ionic thermodynamic properties obtained from the NBS tables. Additionally, the Marcus' and Tissandier's solvation free energy data for Li+, Na+, K+, Cl-, and Br- ions were compared with the conventional solvation free energies by substituting into our self-consistent thermodynamic cycle. Furthermore, Tissandier's absolute solvation free energy data were used as the training set to refit the Lennard-Jones parameters of OPLS-AA force field for ions. To predict salt solubilities, an assumption of a pseudo-solvent was proposed to characterize the coupling work of a solute with its environment from infinitely diluted to saturated solutions, indicating that the Gibbs energy change of solvation process is a function of ionic strength. Following the self-consistency of the cycle, the newly derived formulas were used to determine the salt solubilities by interpolating the intersection of Gibbs free energy of dissolution and the zero free energy line. The refined ion parameters can also predict the structure and thermodynamic properties of aqueous electrolyte solutions, such as densities, pair correlation functions, hydration numbers, mean activity coefficients, vapor pressures, and the radial dependences of the net charge at 298.15 K and 1 bar. Our method can be used to characterize the solid-liquid equilibria of ions or charged particles in aqueous systems. Furthermore, for highly concentrated strong electrolyte systems, it is essential to introduce accurate water models and polarizable force fields.

Does the Sign of Charge Affect the Surface Affinity of Simple Ions?

The role the charge sign of simple ions plays in determining their surface affinity in aqueous solutions is investigated by computer simulation methods. For this purpose, the free surface of aqueous solutions of fictitious salts is simulated at finite concentration both with nonpolarizable point-charge and polarizable Gaussian-charge potential models. The salts consist of monovalent cations and anions that are, apart from the sign of their charge, identical to each other. In particular, we consider the small Na⁺ and the large I^- ions together with their charge-inverted counterparts. In an attempt to avoid the interference even between the behavior of cations and anions, we also simulate systems containing only one of the above ions, and determine the free energy profile of these ions across the liquid-vapor interface of water at infinite dilution by potential of mean force (PMF) calculations. The obtained results reveal that, in the case of small ions, the anion is hydrated considerably stronger than the cation due to the close approach of water H atoms, bearing a positive fractional charge. As a consequence, the surface affinity of a small anion is even smaller than that of its cationic counterpart. However, considering that small ions are effectively repelled from the water surface, the importance of this difference is negligible. Further, a change in the hydration energy trends of the two oppositely charged ions is observed with their increasing size. This change is largely attributed to the fact that, with increasing ion size, the factor of 2 in the magnitude of the fractional charge of the closely approaching water atoms (i.e., O around cations and H around anions) outweighs the closer approach of the H than the O atom in the hydration energy. Thus, for large ions, being already surface active themselves, the surface affinity of the anion is larger than that of its positively charged counterpart. Further, such a difference is seen even in the case when the sign of the surface potential favors the adsorption of cations.

Scaled charges for ions: An improvement but not the final word for modeling electrolytes in water

In this work, we discuss the use of scaled charges when developing force fields for NaCl in water. We shall develop force fields for Na⁺ and Cl^- using the following values for the scaled charge (in electron units): ±0.75, ±0.80, ±0.85, and ±0.92 along with the TIP4P/2005 model of water (for which previous force fields were proposed for q = ±0.85 and q = ±1). The properties considered in this work are densities, structural properties, transport properties, surface tension, freezing point depression, and maximum in density. All the developed models were able to describe quite well the experimental values of the densities. Structural properties were well described by models with charges equal to or larger than ±0.85, surface tension by the charge ±0.92, maximum in density by the charge ±0.85, and transport properties by the charge ±0.75. The use of a scaled charge of ±0.75 is able to reproduce with high accuracy the viscosities and diffusion coefficients of NaCl solutions for the first time. We have also considered the case of KCl in water, and the results obtained were fully consistent with those of NaCl. There is no value of the scaled charge able to reproduce all the properties considered in this work. Although certainly scaled charges are not the final word in the development of force fields for electrolytes in water, its use may have some practical advantages. Certain values of the scaled charge could be the best option when the interest is to describe certain experimental properties.

Dynamics of Aqueous Electrolyte Solutions: Challenges for Simulations

This Perspective article focuses on recent simulation work on the dynamics of aqueous electrolytes. It is well-established that full-charge, nonpolarizable models for water and ions generally predict solution dynamics that are too slow in comparison to experiments. Models with reduced (scaled) charges do better for solution diffusivities and viscosities but encounter issues describing other dynamic phenomena such as nucleation rates of crystals from solution. Polarizable models show promise, especially when appropriately parametrized, but may still miss important physical effects such as charge transfer. First-principles calculations are starting to emerge for these properties that are in principle able to capture polarization, charge transfer, and chemical transformations in solution. While direct ab initio simulations are still too slow for simulations of large systems over long time scales, machine-learning models trained on appropriate first-principles data show significant promise for accurate and transferable modeling of electrolyte solution dynamics.

Improving Condensed-Phase Water Dynamics with Explicit Nuclear Quantum Effects: The Polarizable Q-AMOEBA Force Field

We introduce a new parametrization of the AMOEBA polarizable force field for water denoted Q-AMOEBA, for use in simulations that explicitly account for nuclear quantum effects (NQEs). This study is made possible thanks to the recently introduced adaptive Quantum Thermal Bath (adQTB) simulation technique which computational cost is comparable to classical molecular dynamics. The flexible Q-AMOEBA model conserves the initial AMOEBA functional form, with an intermolecular potential including an atomic multipole description of electrostatic interactions (up to quadrupole), a polarization contribution based on the Thole interaction model and a buffered 14-7 potential to model van der Waals interactions. It has been obtained by using a ForceBalance fitting strategy including high-level quantum chemistry reference energies and selected condensed-phase properties targets. The final Q-AMOEBA model is shown to accurately reproduce both gas-phase and condensed-phase properties, notably improving the original AMOEBA water model. This development allows the fine study of NQEs on water liquid phase properties such as the average H-O-H angle compared to its gas-phase equilibrium value, isotope effects, and so on. Q-AMOEBA also provides improved infrared spectroscopy prediction capabilities compared to AMOEBA03. Overall, we show that the impact of NQEs depends on the underlying model functional form and on the associated strength of hydrogen bonds. Since adQTB simulations can be performed at near classical computational cost using the Tinker-HP package, Q-AMOEBA can be extended to organic molecules, proteins, and nucleic acids opening the possibility for the large-scale study of the importance of NQEs in biophysics.

Melting points of water models: Current situation

By using the direct coexistence method, we have calculated the melting points of ice I h at normal pressure for three recently proposed water models, namely, TIP3P-FB, TIP4P-FB, and TIP4P-D. We obtained T m = 216 K for TIP3P-FB, T m = 242 K for TIP4P-FB, and T m = 247 K for TIP4P-D. We revisited the melting point of TIP4P/2005 and TIP5P obtaining T m = 250 and 274 K, respectively. We summarize the current situation of the melting point of ice I h for a number of water models and conclude that no model is yet able to simultaneously reproduce the melting temperature of ice I h and the temperature of the maximum in density at room pressure. This probably points toward our both still incomplete knowledge of the potential energy surface of water and the necessity of incorporating nuclear quantum effects to describe both properties simultaneously.

Maximum in density of electrolyte solutions: Learning about ion-water interactions and testing the Madrid-2019 force field

In this work, we studied the effect of Li⁺, Na⁺, K⁺, Mg²⁺, and Ca²⁺ chlorides and sulfates on the temperature of maximum density (TMD) of aqueous solutions at room pressure. Experiments at 1 molal salt concentration were carried out to determine the TMD of these solutions. We also performed molecular dynamics simulations to estimate the TMD at 1 and 2 m with the Madrid-2019 force field, which uses the TIP4P/2005 water model and scaled charges for the ions, finding an excellent agreement between experiment and simulation. All the salts studied in this work shift the TMD of the solution to lower temperatures and flatten the density vs temperature curves (when compared to pure water) with increasing salt concentration. The shift in the TMD depends strongly on the nature of the electrolyte. In order to explore this dependence, we have evaluated the contribution of each ion to the shift in the TMD concluding that Na⁺, Ca²⁺, and SO₄ ^2- seem to induce the largest changes among the studied ions. The volume of the system has been analyzed for salts with the same anion and different cations. These curves provide insight into the effect of different ions upon the structure of water. We claim that the TMD of electrolyte solutions entails interesting physics regarding ion-water and water-water interactions and should, therefore, be considered as a test property when developing force fields for electrolytes. This matter has been rather unnoticed for almost a century now and we believe it is time to revisit it.

Electron-induced fragmentation of water droplets: Simulation study

The transport of free electrons in a water environment is still poorly understood. We show that additional insight can be brought about by investigating fragmentation patterns of finite-size particles upon electron impact ionization. We have developed a composite protocol aiming to simulate fragmentation of water clusters by electrons with kinetic energies in the range of up to 100 eV. The ionization events for atomistically described molecular clusters are identified by a kinetic Monte Carlo procedure. We subsequently model the fragmentation with classical molecular dynamics simulations, calibrated by non-adiabatic quantum mechanics/molecular mechanics simulations of the ionization process. We consider one-electron ionizations, energy transfer via electronic excitation events, elastic scattering, and also the autoionization events through intermolecular Coulombic decay. The simulations reveal that larger water clusters are often ionized repeatedly, which is the cause of substantial fragmentation. After losing most of its energy, low-energy electrons further contribute to fragmentation by electronic excitations. The simultaneous measurement of cluster size distribution before and after the ionization represents a sensitive measure of the energy transferred into the system by an incident electron.

Freezing point depression of salt aqueous solutions using the Madrid-2019 model

Salt aqueous solutions are relevant in many fields, ranging from biological systems to seawater. Thus, the availability of a force-field that is able to reproduce the thermodynamic and dynamic behavior of salt aqueous solutions would be of great interest. Unfortunately, this has been proven challenging, and most of the existing force-fields fail to reproduce much of their behavior. In particular, the diffusion of water or the salt solubility are often not well reproduced by most of the existing force-fields. Recently, the Madrid-2019 model was proposed, and it was shown that this force-field, which uses the TIP4P/2005 model for water and non-integer charges for the ions, provides a good description of a large number of properties, including the solution densities, viscosities, and the diffusion of water. In this work, we assess the performance of this force-field on the evaluation of the freezing point depression. Although the freezing point depression is a colligative property that at low salt concentrations depends solely on properties of pure water, a good model for the electrolytes is needed to accurately predict the freezing point depression at moderate and high salt concentrations. The coexistence line between ice and several salt aqueous solutions (NaCl, KCl, LiCl, MgCl₂, and Li₂SO₄) up to the eutectic point is estimated from direct coexistence molecular dynamics simulations. Our results show that this force-field reproduces fairly well the experimentally measured freezing point depression with respect to pure water freezing for all the salts and at all the compositions considered.

Comparison of molecular dynamics simulations of water with neutron and X-ray scattering experiments

The atomistic structure and dynamics obtained from molecular dynamics (MD) simulations with the example of TIP3P (rigid and flexible) and TIP4P/2005 (rigid) water is compared to neutron and X-ray scattering data at ambient conditions. Neutron and X-ray diffractograms are calculated from the simulations for four isotopic substitutions as well as the incoherent intermediate scattering function for neutrons. The resulting curves are compared to each other and to published experimental data. Differences between simulated and measured intermediate scattering functions are quantified by fitting an analytic model to the computed values. The sensitivity of the scattering curves to the parameters of the MD simulations is demonstrated on the example of two parameters, bond length and angle.

Investigating the Accuracy of Water Models through the Van Hove Correlation Function

We present molecular-simulation-based calculations of the Van Hove correlation function (VHF) of water using multiple modeling approaches: classical molecular dynamics with simple three-site nonpolarizable models, with a polarizable model, and with a reactive force field; density functional tight-binding molecular dynamics; and ab initio molecular dynamics. Due to the many orders of magnitude difference in the computational cost of these approaches, we investigate how small and short the simulations can be while still yielding sufficiently accurate and interpretable results for the VHF. We investigate the accuracy of the different models by comparing them to recently published inelastic X-ray scattering measurements of the VHF. We find that all of the models exhibit qualitative agreement with the experiments, and in some models and for some properties, the agreement is quantitative. This work lays the foundation for future simulation approaches to calculating the VHF for aqueous solutions in bulk and under nanoconfinement.

Charge State Dependence of Amino Acid Propensity at Water Surface: Mechanisms Elucidated by Molecular Dynamics Simulations

Atmospheric aerosols contain a variety of compounds, among them free amino acids and salt ions. The pH of the aerosol droplets depends on their origin and environment. Consequently, compounds like free amino acids found in the droplets will be at different charge states, since these states to a great extent depend on the surrounding pH condition. In droplets of marine origin, amino acids are believed to drive salt ions to the water surface and a pH-dependent amino acid surface propensity will, therefore, indirectly affect many processes in atmospheric chemistry and physics such as for instance cloud condensation. To understand the surface propensity of glycine, valine, and phenylalanine at acidic, neutral, and basic pH, we used molecular dynamics (MD) simulations to investigate them at three different charge states in water. Their respective surface propensities were obtained by the means of a potential of mean force (PMF) in an umbrella sampling approach. Glycine was found to have no preference for the surface, while both valine and phenylalanine showed high propensities. Among the charge states of the surface-enriched ones, the cation, representing the amino acids at low pH, was found to have the highest affinity. Free energy decomposition revealed that the driving forces depend strongly on the nature of the amino acid and its charge state. In phenylalanine, the main factor was found to be a substantial entropy gain, likely related to the side chain, whereas in valine, hydrogen bonding to the functional groups leads to favorable energies and, in turn, affects the surface propensity. A significant gain in water-water enthalpy was seen for both valine and phenylalanine.

Vapor-liquid equilibrium of water with the MB-pol many-body potential

Among the many existing molecular models of water, the MB-pol many-body potential has emerged as a remarkably accurate model, capable of reproducing thermodynamic, structural, and dynamic properties across water's solid, liquid, and vapor phases. In this work, we assessed the performance of MB-pol with respect to an important set of properties related to vapor-liquid coexistence and interfacial behavior. Through direct coexistence classical molecular dynamics simulations at temperatures of 400 K < T < 600 K, we calculated properties such as equilibrium coexistence densities, vapor-liquid interfacial tension, vapor pressure, and enthalpy of vaporization and compared the MB-pol results to experimental data. We also compared rigid vs fully flexible variants of the MB-pol model and evaluated system size effects for the properties studied. We found that the MB-pol model predictions are in good agreement with experimental data, even for temperatures approaching the vapor-liquid critical point; this agreement was largely insensitive to system sizes or the rigid vs flexible treatment of the intramolecular degrees of freedom. These results attest to the chemical accuracy of MB-pol and its high degree of transferability, thus enabling MB-pol's application across a large swath of water's phase diagram.

Single-Particle Dynamics at the Intrinsic Surface of Aqueous Alkali Halide Solutions

The distribution of ions in the proximity of the liquid-vapor interface of their aqueous solution has been the subject of an intense debate during the last decade. The effects of ionic polarizability have been one of its salient aspects. Much less has been said about the corresponding dynamical properties, which are substantially unexplored. Here, we investigate the single-particle dynamics at the liquid-vapor interface of several alkali halide solutions, using molecular dynamics simulations with polarizable and nonpolarizable force fields and intrinsic surface analysis. We analyze the diffusion coefficient, residence time, and velocity autocorrelation function of water and ions and investigate how these properties depend on the molecular layer where they reside. While anions are found in the first molecular layer for relatively long times, cations are only making quick excursions into it, thanks to thermal fluctuations. The in-layer residence time of ions and their molar fraction in the layer turned out to be linearly dependent on each other. We interpret this unexpected result using a simple two-state model. In addition, we found that, unlike water and other neat molecular liquids that show a different diffusion mechanism at the surface than in the bulk of their liquid phase, ions do not enjoy enhanced mobility in the surface layer of their aqueous solution. This result indicates that ions in the surface layer are shielded by their nearest water neighbors from being exposed to the vapor phase as much as possible. Such positions are available for the ions at the negatively curved troughs of the molecularly rugged liquid surface.

Linker residues regulate the activity and stability of hexokinase 2, a promising anticancer target

Hexokinase (HK) catalyzes the first step in glucose metabolism, making it an exciting target for the inhibition of tumor initiation and progression due to their elevated glucose metabolism. The upregulation of hexokinase-2 (HK2) in many cancers and its limited expression in normal tissues make it a particularly attractive target for the selective inhibition of cancer growth and the eradication of tumors with limited side effects. The design of such safe and effective anticancer therapeutics requires the development of HK2-specific inhibitors that will not interfere with other HK isozymes. As HK2 is unique among HKs in having a catalytically active N-terminal domain (NTD), we have focused our attention on this region. We previously found that NTD activity is affected by the size of the linker helix-α₁₃ that connects the N- and C-terminal domains of HK2. Three nonactive site residues (D447, S449, and K451) at the beginning of the linker helix-α₁₃ have been found to regulate the NTD activity of HK2. Mutation of these residues led to increased dynamics, as shown via hydrogen deuterium exchange analysis and molecular dynamic simulations. D447A contributed the most to the enhanced dynamics of the NTD, with reduced calorimetric enthalpy of HK2. Similar residues exist in the C-terminal domain (CTD) but are unnecessary for HK1 and HK2 activity. Thus, we postulate these residues serve as a regulatory site for HK2 and may provide new directions for the design of anticancer therapeutics that reduce the rate of glycolysis in cancer through specific inhibition of HK2.

Surface Affinity of Alkali and Halide Ions in Their Aqueous Solution: Insight from Intrinsic Density Analysis

The surface tension of all aqueous alkali halide solutions is higher than that of pure water. According to the Gibbs adsorption equation, this indicates a net depletion of these ions in the interfacial region. However, simulations and experiments show that large, soft ions, such as I^-, can accumulate at the liquid/vapor interface. The presence of a loose hydration shell is usually considered to be the reason for this behavior. In this work, we perform computer simulations to characterize the liquid-vapor interface of aqueous alkali chloride and sodium halide solutions systematically, considering all ions from Li⁺ to Cs⁺ and from F^- to I^-. Using computational methods for the removal of surface fluctuations, we analyze the structure of the interface at a dramatically enhanced resolution, showing that the positive excess originates in the very first molecular layer and that the next 3-4 layers account for the net negative excess. With the help of a fictitious system with charge-inverted ion pairs, we also show that it is not possible to rationalize the surface affinity of ions in solutions in terms of the properties of anions and cations separately. Moreover, the surface excess is generally dominated by the smaller of the two ions.

Pressure in Molecular Simulations with Scaled Charges. 1. Ionic Systems

Charge scaling, rationalized as MDEC (molecular dynamics in electronic continuum) or ECC (electronic continuum correction), has become a widely used simple approach to how to avoid self-consistent induced dipoles yet approximately take into account the effects of electronic polarizability. It has been assumed that the continuum permittivity does not depend on density; in turn, pressure is calculated by standard formulas. In this work, we elaborate a complementary approximation of density-independent molecular polarizability and derive formulas for pressure corrections within the MDEC framework; real behavior lies between these two extremes. The pressure corrections for test ionic systems are huge and negative, leading to sizable densities in constant-pressure MDEC simulations. A comparison of MDEC results with equivalent polarizable systems gives a good pressure match for a crystal but very low MDEC pressures for ionic liquids. These results witness about the importance of a correct density dependence not only of continuum permittivity in MDEC simulations but also of polarizability in polarizable simulations.

Support Vector Regression-Based Monte Carlo Simulation of Flexible Water Clusters

Molecular simulations based on classical force fields are computationally efficient but lack accuracy due to the empirical formulation of non-bonded interactions. Quantum mechanical (QM) methods, albeit accurate, have inhibitory computational costs for large molecules and clusters. Hence, to overcome the bottleneck, machine learning (ML)-based methods have been employed in the recent years. We had earlier reported a combined scheme of many-body expansion (MBE) and ML to predict the interaction energies of rigid water clusters. In this work, we proceed toward building a flexible water model using the ML-MBE scheme. This ML-MBE scheme has an error of <1% for interaction energy prediction in comparison to the parent QM method for flexible water decamers. Machine learning-based Monte Carlo simulations (MLMC) are performed with this water model, and the structural properties of these configurations are compared with those obtained from ab initio molecular dynamics (AIMD) and the TIP3P classical force field. The radial distribution functions, tetrahedral order parameters, and number of hydrogen bonds in AIMD and MLMC have a similar qualitative and quantitative trend, whereas the classical force fields show a significant deviation.

Polarizable TIP7P Water Model with Perturbation Charges Evaluated from ABEEM

A polarizable version of the rigid seven-site (TIP7P) water model with the atom-bond electronegativity equalization method (ABEEM) is proposed. The model uses direct polarization, where an isolated water monomer in the equilibrium geometry is assumed as a reference state and the polarization of the monomer arises from interacting with other molecules as a perturbation of the reference state. The charge on each site of the monomer splits into reference charge and perturbation charge. The perturbation charge arises only because of other reference charges. The interaction of the perturbation charge with other perturbation charges is replaced using polarization scaling to enhance the interaction of perturbation charge with the reference charges of the sites from other molecules. The perturbation charges are updated by evaluating explicit expressions once. This direct polarization is time-reversible because the charge update is independent of the charges in previous simulation steps. A Slater-type damping function moderates the short-range electrostatics to treat charge diffusion. The Ewald method corrects the long-range electrostatics both in the nuclei movement and in electronegativity equalization to diminish the size effect. The water model is parameterized by fitting the ab initio results of water clusters and the experimental results of water monomers and thermodynamic properties for liquid water. Owing to polarizability, the model performs better than the TIP7P model in terms of vaporization enthalpy, isothermal compressibility, and shear viscosity of the liquid phase. It performs better at the melting point of ice but slightly worse under critical conditions than the TIP7P model. Direct polarization has a low time complexity of O(N) and is a good choice for ABEEM to improve its computational efficiency.

A FFLUX Water Model: Flexible, Polarizable and with a Multipolar Description of Electrostatics

Key to progress in molecular simulation is the development of advanced models that go beyond the limitations of traditional force fields that employ a fixed, point charge-based description of electrostatics. Taking water as an example system, the FFLUX framework is shown capable of producing models that are flexible, polarizable and have a multipolar description of the electrostatics. The kriging machine-learning methods used in FFLUX are able to reproduce the intramolecular potential energy surface and multipole moments of a single water molecule with chemical accuracy using as few as 50 training configurations. Molecular dynamics simulations of water clusters (25-216 molecules) using the new FFLUX model reveal that incorporating charge-quadrupole, dipole-dipole, and quadrupole-charge interactions into the description of the electrostatics results in significant changes to the intermolecular structuring of the water molecules. © 2019 The Authors. Journal of Computational Chemistry published by Wiley Periodicals, Inc.

Translational and Rotational Diffusion in Liquid Water at Very High Pressure: A Simulation Study

Translational and rotational diffusion coefficients of liquid water have been computed from molecular dynamics simulation with a recent polarizable potential at 298, 400, and 550 K at very high pressure. At 298 K, the model reproduces the initial increase and the occurrence of a maximum for the translational and rotational diffusion coefficients when the pressure increases. At 400 and 550 K, translational and rotational diffusion coefficients are found to monotonically decrease when pressure increases in the gigapascal range, with the translational coefficient decreasing faster than the rotational one. At 400 K, such an evolution of the rotational diffusion coefficient contrasts with quasielastic neutron scattering results predicting a near independence of the rotational diffusion with a pressure increase above ≃0.5 GPa. An interpretation is proposed to explain this discrepancy. The pressure dependence of the translation-rotation coupling is analyzed. The anisotropy of rotational diffusion is investigated by computing the rotational diffusion tensor in a molecular system of axes and the reorientational correlation times of rank 1 and rank 2 of the inertia axes and of the OH bond vector. Deviation of the simulation data with respect to the predictions of the isotropic Debye model of rotational diffusion are quantified and can be used to estimate experimental rotational diffusion coefficients from experimental reorientational correlation times.

A force field of Li+, Na+, K+, Mg2+, Ca2+, Cl-, and SO4 2- in aqueous solution based on the TIP4P/2005 water model and scaled charges for the ions

In this work, a force field for several ions in water is proposed. In particular, we consider the cations Li⁺, Na⁺, K⁺, Mg²⁺, and Ca²⁺ and the anions Cl^- and SO₄ ^2-. These ions were selected as they appear in the composition of seawater, and they are also found in biological systems. The force field proposed (denoted as Madrid-2019) is nonpolarizable, and both water molecules and sulfate anions are rigid. For water, we use the TIP4P/2005 model. The main idea behind this work is to further explore the possibility of using scaled charges for describing ionic solutions. Monovalent and divalent ions are modeled using charges of 0.85 and 1.7, respectively (in electron units). The model allows a very accurate description of the densities of the solutions up to high concentrations. It also gives good predictions of viscosities up to 3 m concentrations. Calculated structural properties are also in reasonable agreement with the experiment. We have checked that no crystallization occurred in the simulations at concentrations similar to the solubility limit. A test for ternary mixtures shows that the force field provides excellent performance at an affordable computer cost. In summary, the use of scaled charges, which could be regarded as an effective and simple way of accounting for polarization (at least to a certain extend), improves the overall description of ionic systems in water. However, for purely ionic systems, scaled charges will not adequately describe neither the solid nor the melt.

Glass polymorphism and liquid-liquid phase transition in aqueous solutions: experiments and computer simulations

One of the most intriguing anomalies of water is its ability to exist as distinct amorphous ice forms (glass polymorphism or polyamorphism). This resonates well with the possible first-order liquid-liquid phase transition (LLPT) in the supercooled state, where ice is the stable phase. In this Perspective, we review experiments and computer simulations that search for LLPT and polyamorphism in aqueous solutions containing salts and alcohols. Most studies on ionic solutes are devoted to NaCl and LiCl; studies on alcohols have mainly focused on glycerol. Less attention has been paid to protein solutions and hydrophobic solutes, even though they reveal promising avenues. While all solutions show polyamorphism and an LLPT only in dilute, sub-eutectic mixtures, there are differences regarding the nature of the transition. Isocompositional transitions for varying mole fractions are observed in alcohol but not in ionic solutions. This is because water can surround alcohol molecules either in a low- or high-density configuration whereas for ionic solutes, the water ion hydration shell is forced into high-density structures. Consequently, the polyamorphic transition and the LLPT are prevented near the ions, but take place in patches of water within the solutions. We highlight discrepancies and different interpretations within the experimental community as well as the key challenges that need consideration when comparing experiments and simulations. We point out where reinterpretation of past studies helps to draw a unified, consistent picture. In addition to the literature review, we provide original experimental results. A list of eleven open questions that need further consideration is identified.

Molecular Dynamics Simulations of Ionic Liquids and Electrolytes Using Polarizable Force Fields

Many applications in chemistry, biology, and energy storage/conversion research rely on molecular simulations to provide fundamental insight into structural and transport properties of materials with high ionic concentrations. Whether the system is comprised entirely of ions, like ionic liquids, or is a mixture of a polar solvent with a salt, e.g., liquid electrolytes for battery applications, the presence of ions in these materials results in strong local electric fields polarizing solvent molecules and large ions. To predict properties of such systems from molecular simulations often requires either explicit or mean-field inclusion of the influence of polarization on electrostatic interactions. In this manuscript, we review the pros and cons of different treatments of polarization ranging from the mean-field approaches to the most popular explicit polarization models in molecular dynamics simulations of ionic materials. For each method, we discuss their advantages and disadvantages and emphasize key assumptions as well as their adjustable parameters. Strategies for the development of polarizable models are presented with a specific focus on extracting atomic polarizabilities. Finally, we compare simulations using polarizable and nonpolarizable models for several classes of ionic systems, discussing the underlying physics that each approach includes or ignores, implications for implementation and computational efficiency, and the accuracy of properties predicted by these methods compared to experiments.

New scaling relations to compute atom-in-material polarizabilities and dispersion coefficients: part 1. Theory and accuracy

Polarizabilities and London dispersion forces are important to many chemical processes. Force fields for classical atomistic simulations can be constructed using atom-in-material polarizabilities and C n (n = 6, 8, 9, 10…) dispersion coefficients. This article addresses the key question of how to efficiently assign these parameters to constituent atoms in a material so that properties of the whole material are better reproduced. We develop a new set of scaling laws and computational algorithms (called MCLF) to do this in an accurate and computationally efficient manner across diverse material types. We introduce a conduction limit upper bound and m-scaling to describe the different behaviors of surface and buried atoms. We validate MCLF by comparing results to high-level benchmarks for isolated neutral and charged atoms, diverse diatomic molecules, various polyatomic molecules (e.g., polyacenes, fullerenes, and small organic and inorganic molecules), and dense solids (including metallic, covalent, and ionic). We also present results for the HIV reverse transcriptase enzyme complexed with an inhibitor molecule. MCLF provides the non-directionally screened polarizabilities required to construct force fields, the directionally-screened static polarizability tensor components and eigenvalues, and environmentally screened C6 coefficients. Overall, MCLF has improved accuracy compared to the TS-SCS method. For TS-SCS, we compared charge partitioning methods and show DDEC6 partitioning yields more accurate results than Hirshfeld partitioning. MCLF also gives approximations for C8, C9, and C10 dispersion coefficients and quantum Drude oscillator parameters. This method should find widespread applications to parameterize classical force fields and density functional theory (DFT) + dispersion methods.

Molecular polarizability in open ensemble simulations of aqueous nanoconfinements under electric field

Molecular polarization at aqueous interfaces involves fast degrees of freedom that are often averaged-out in atomistic-modeling approaches. The resulting effective interactions depend on a specific environment, making explicit account of molecular polarizability particularly important in solutions with pronounced anisotropic perturbations, including solid/liquid interfaces and external fields. Our work concerns polarizability effects in nanoscale confinements under electric field, open to an unperturbed bulk environment. We model aqueous molecules and ions in hydrophobic pores using the Gaussian-charge-on-spring BK3-AH representation. This involves nontrivial methodology developments in expanded ensemble Monte Carlo simulations for open systems with long-ranged multibody interactions and necessitates further improvements for efficient modeling of polarizable ions. Structural differences between fixed-charge and polarizable models were captured in molecular dynamics simulations for a set of closed systems. Our open ensemble results with the BK3 model in neat-aqueous systems capture the ∼10% reduction of molecular dipoles within the surface layer near the hydrophobic pore walls in analogy to reported quantum mechanical calculations at water/vapor interfaces. The polarizability affects the interfacial dielectric behavior and weakens the electric-field dependence of water absorption at pragmatically relevant porosities. We observe moderate changes in thermodynamic properties and atom and charged-site spatial distributions; the Gaussian distribution of mobile charges on water and ions in the polarizable model shifts the density amplitudes and blurs the charge-layering effects associated with increased ion absorption. The use of polarizable force field indicates an enhanced response of interfacial ion distributions to applied electric field, a feature potentially important for in silico modeling of electric double layer capacitors.