The quantum mechanics-based polarizable force field for water simulations

The heterogeneity of aqueous solutions: the current situation in the context of experiment and theory

The advancement of experimental methods has provided new information about the structure and structural fluctuations of water. Despite the appearance of numerous models, which aim to describe a wide range of thermodynamic and electrical characteristics of water, there is a deficit in systemic understanding of structuring in aqueous solutions. A particular challenge is the fact that even pure water is a heterogeneous, multicomponent system composed of molecular and supramolecular structures. The possibility of the existence of such structures and their nature are of fundamental importance for various fields of science. However, great difficulties arise in modeling relatively large supramolecular structures (e.g. extended hydration shells), where the bonds between molecules are characterized by low energy. Generally, such structures may be non-equilibrium but relatively long-lived. Evidently, the short times of water microstructure exchanges do not mean short lifetimes of macrostructures, just as the instability of individual parts does not mean the instability of the entire structure. To explain this paradox, we review the data from experimental and theoretical research. Today, only some of the experimental results on the lifetime of water structures have been confirmed by modeling, so there is not a complete theoretical picture of the structure of water yet. We propose a new hierarchical water macrostructure model to resolve the issue of the stability of water structures. In this model, the structure of water is presented as consisting of many hierarchically related levels (the stratification model). The stratification mechanism is associated with symmetry breaking at the formation of the next level, even with minimal changes in the properties of the previous level. Such a hierarchical relationship can determine the unique physico-chemical properties of water systems and, in the future, provide a complete description of them.

Quantum mechanics based non-bonded force field functions for use in molecular dynamics simulations of materials and systems: The nitrogen and oxygen columns

Accurate Force Fields (FFs) are essential for Molecular Dynamics (MD) simulations of the dynamics of realistic materials in terms of atomic-level interactions. The FF parameters of short-range valence interactions can be derived through Quantum Mechanical (QM) calculations on model systems practical for QM (<300 atoms). Similarly, the dynamic electrostatic interactions can be described with methods such as QEq or PQEq that allow charges and polarization to adjust dynamically. However, accurately extracting long-range van der Waals (vdW) interactions from QM calculations poses challenges due to the absence of a definitive method to distinguish between the different energetic components of electrostatics, polarization, vdW, hydrogen bonding, and valence interactions. To do this we use the Perdew-Burke-Ernzerhof flavor of Density Functional Theory, including empirical D3 vdW corrections, to predict the Equation of State for each element (keeping any covalent bonds fixed), from which we obtain the two-body vdW nonbond potential. Here, we extend these calculations to include non-bonded parameters for the N and O columns of the periodic table so that we now describe columns 15 (N), 16 (O), 17 (F), and 18 (Ne) of the periodic table. For these 20 elements, we find that the two-body vdW potentials can all be mapped to a single universal two-body curve, with just three scaling parameters: Re, De, and L. We refer to this as the Universal NonBond (UNB) potential. We expect this to be useful for new MD simulations and a helpful starting point to obtain UNB parameters for the remainder of the periodic table.

Prediction of 1H Singlet Relaxation via Intermolecular Dipolar Couplings Using the Molecular Dynamics Method

Dissolution dynamic nuclear polarization has been applied in various fields, including chemistry, biology, and medical science. To expand the scope of these applications, the nuclear singlet state, which is decoherence-free against dipolar relaxation between spin pairs, has been studied experimentally, theoretically, and numerically. The singlet state composed of proton spins is used in several applications, such as enhanced polarization preservation, molecular tagging to probe slow dynamic processes, and detection of ligand-protein complexes. In this study, we predict the lifetimes of the nuclear spin states composed of proton spin pairs using the molecular dynamics method and quantum chemistry simulations. We consider intramolecular dipolar, intermolecular dipolar between solvent and solute, chemical shift anisotropy, and spin-rotation interactions. In particular, the relaxation rate of intermolecular dipolar interactions is calculated using the molecular dynamics method for various solvents. The calculated values and the experimental values are of the same order of magnitude. Our program would provide insight into the molecular design of several NMR applications and would be helpful in predicting the nuclear spin relaxation time of synthetic molecules in advance.

Singular value decomposition analysis of the electron density changes occurring upon electrostatic polarization of water

In-depth elucidation of how molecules are electrically polarized would be one key factor for understanding the properties of those molecules under various thermodynamic and/or spatial conditions. Here this problem is tackled for the case of hydrogen-bonded water by conducting singular value decomposition of the electron density changes that occur upon electrostatic polarization. It is shown that all those electron density changes are approximately described as linear combinations of ten orthonormal basis “vectors”. One main component is the interatomic charge transfer through each OH bond, while some others are characterized as the atomic dipolar polarizations, meaning that both of these components are important for the electrostatic polarization of water. The interaction parameters that reasonably well reproduce the induced dipole moments are derived, which indicate the extent of mixing of the two components in electrostatic polarization.

The main features of the electron density changes that occur upon electrostatic polarization of water are elucidated by conducting singular value decomposition analysis of those changes.

Implicit Solvation Methods for Catalysis at Electrified Interfaces

Implicit solvation is an effective, highly coarse-grained approach in atomic-scale simulations to account for a surrounding liquid electrolyte on the level of a continuous polarizable medium. Originating in molecular chemistry with finite solutes, implicit solvation techniques are now increasingly used in the context of first-principles modeling of electrochemistry and electrocatalysis at extended (often metallic) electrodes. The prevalent ansatz to model the latter electrodes and the reactive surface chemistry at them through slabs in periodic boundary condition supercells brings its specific challenges. Foremost this concerns the difficulty of describing the entire double layer forming at the electrified solid–liquid interface (SLI) within supercell sizes tractable by commonly employed density functional theory (DFT). We review liquid solvation methodology from this specific application angle, highlighting in particular its use in the widespread ab initio thermodynamics approach to surface catalysis. Notably, implicit solvation can be employed to mimic a polarization of the electrode’s electronic density under the applied potential and the concomitant capacitive charging of the entire double layer beyond the limitations of the employed DFT supercell. Most critical for continuing advances of this effective methodology for the SLI context is the lack of pertinent (experimental or high-level theoretical) reference data needed for parametrization.

Two Faces of the Two-Phase Thermodynamic Model

The quantum harmonic model and the two-phase thermodynamic method (2PT) are widely used to obtain quantum-corrected properties such as isobaric heat capacities or molar entropies. 2PT heat capacities were calculated inconsistently in the literature. For water, the classical heat capacity was also considered, but for organic liquids, it was omitted. We reanalyzed the performance of different quantum corrections on the heat capacities of common organic solvents against experimental data. We have pointed out serious flaws in previous 2PT studies. The vibrational density of states was calculated incorrectly causing a 39% relative error in diffusion coefficients and 45% error in the 2PT heat capacities. The wrong conversion of isobaric and isochoric heat capacities also caused about 40% error but in the other direction. We have introduced the concept of anharmonic correction (AC), which is simply the deviation of the classical heat capacity from that of the harmonic oscillator model. This anharmonic contribution is around +30 to 40 J/(mol K) for water depending on the water model and −8 to −10 J/(mol K) for hydrocarbons and halocarbons. AC is unrealistically large, +40 J/(K mol) for alcohols and amines, indicating some deficiency of the OPLS force field. The accuracy of the computations was also assessed with the determination of the self-diffusion coefficients.

DL_FFLUX: A Parallel, Quantum Chemical Topology Force Field

DL_FFLUX is a force field based on quantum chemical topology that can perform molecular dynamics for flexible molecules endowed with polarizable atomic multipole moments (up to hexadecapole). Using the machine learning method kriging (aka Gaussian process regression), DL_FFLUX has access to atomic properties (energy, charge, dipole moment, etc.) with quantum mechanical accuracy. Newly optimized and parallelized using domain decomposition Message Passing Interface (MPI), DL_FFLUX is now able to deliver this rigorous methodology at scale while still in reasonable time frames. DL_FFLUX is delivered as an add-on to the widely distributed molecular dynamics code DL_POLY 4.08. For the systems studied here (10³–10⁵ atoms), DL_FFLUX is shown to add minimal computational cost to the standard DL_POLY package. In fact, the optimization of the electrostatics in DL_FFLUX means that, when high-rank multipole moments are enabled, DL_FFLUX is up to 1.25× faster than standard DL_POLY. The parallel DL_FFLUX preserves the quality of the scaling of MPI implementation in standard DL_POLY. For the first time, it is feasible to use the full capability of DL_FFLUX to study systems that are large enough to be of real-world interest. For example, a fully flexible, high-rank polarized (up to and including quadrupole moments) 1 ns simulation of a system of 10 125 atoms (3375 water molecules) takes 30 h (wall time) on 18 cores.

Entropic Stabilization of Water at Graphitic Interfaces

The thermodynamic stability of water next to graphitic surfaces is of fundamental interest, as it underlies several natural phenomena and important industrial processes. It is commonly assumed that water wets graphite more than graphene due to increased, favorable van der Waals interactions between the interfacial water molecules with multiple carbon sheets. Here, we employed extensive computer simulations and analysis of the molecular correlation functions to show that the interfacial water thermodynamics is in fact dominated by surface entropy. We show that on graphite, destabilization of the interfacial hydrogen bond network leads to an overcompensating increase in population of low frequency translational and librational modes, which is ultimately responsible for the increased interfacial stability compared to graphene. The spectroscopic signature of this effect is an enhancement of the modes near 100 and 300 cm^-1. This subtle interplay between entropy and surface binding may have important consideration for interpretations of various phenomena, including the hydrophobic effect.

The phase diagram of carbon dioxide from correlation functions and a many-body potential

The phase stability and equilibria of carbon dioxide are investigated from 125-325 K and 1-10 000 atm using extensive molecular dynamics (MD) simulations and the Two-Phase Thermodynamics (2PT) method. We devise a direct approach for calculating phase diagrams, in general, by considering the separate chemical potentials of the isolated phase at specific points on the P-T diagram. The unique ability of 2PT to accurately and efficiently approximate the entropy and Gibbs energy of liquids allows for assignment of phase boundaries from relatively short (∼100 ps) MD simulations. We validate our approach by calculating the critical properties of the flexible elementary physical model 2, showing good agreement with previous results. We show, however, that the incorrect description of the short-range Pauli force and the lack of molecular charge polarization lead to deviations from experiments at high pressures. We, thus, develop a many-body, fluctuating charge model for CO₂, termed CO₂-Fq, from high level quantum mechanics (QM) calculations that accurately capture the condensed phase vibrational properties of the solid (including the Fermi resonance at 1378 cm^-1) as well as the diffusional properties of the liquid, leading to overall excellent agreement with experiments over the entire phase diagram. This work provides an efficient computational approach for determining phase diagrams of arbitrary systems and underscores the critical role of QM charge reorganization physics in molecular phase stability.

Deconstructing the Local Intermolecular Ordering and Dynamics of Liquid Chloroform and Bromoform

Local intermolecular structure and dynamics of the polar molecular liquids chloroform and bromoform are studied by molecular dynamics simulation. Structural distribution functions, including 1- and 2-D pair correlations and dipole contour plots allow direct comparison and show agreement with recent analyses of diffraction experiments. Studies of the haloforms' reorientational dynamics and longevity of structural features resulting from intermolecular interaction extend previous work toward deeper understanding of the factors controlling these features. Analyses of ensemble average structures and dynamical properties isolate mass, electrostatics, and steric packing as driving forces or contributing factors for the observed ordering and dynamics.

Signatures of a liquid-liquid transition in an ab initio deep neural network model for water

The possible existence of a metastable liquid-liquid transition (LLT) and a corresponding liquid-liquid critical point (LLCP) in supercooled liquid water remains a topic of much debate. An LLT has been rigorously proved in three empirically parametrized molecular models of water, and evidence consistent with an LLT has been reported for several other such models. In contrast, experimental proof of this phenomenon has been elusive due to rapid ice nucleation under deeply supercooled conditions. In this work, we combined density functional theory (DFT), machine learning, and molecular simulations to shed additional light on the possible existence of an LLT in water. We trained a deep neural network (DNN) model to represent the ab initio potential energy surface of water from DFT calculations using the Strongly Constrained and Appropriately Normed (SCAN) functional. We then used advanced sampling simulations in the multithermal-multibaric ensemble to efficiently explore the thermophysical properties of the DNN model. The simulation results are consistent with the existence of an LLCP, although they do not constitute a rigorous proof thereof. We fit the simulation data to a two-state equation of state to provide an estimate of the LLCP's location. These combined results-obtained from a purely first-principles approach with no empirical parameters-are strongly suggestive of the existence of an LLT, bolstering the hypothesis that water can separate into two distinct liquid forms.

Ten Facets, One Force Field: The GAL19 Force Field for Water-Noble Metal Interfaces

Understanding the structure of the water/metal interfaces plays an important role in many areas ranging from surface chemistry to environmental processes. The size, required phase-space sampling, and the slow diffusion of molecules at the water/metal interfaces motivate the development of accurate force fields. We develop and parametrize GAL19, a novel force field, to describe the interaction of water with two facets (111 and 100) of five metals (Pt, Pd, Au, Ag, Cu). To increase transferability compared to its predecessor GAL17, the water-metal interaction is described as a sum of pairwise terms. The interaction energy has three contributions: (i) physisorption is described via a Tang and Toennies potential, (ii) chemisorption and surface corrugation rely on an attractive Gaussian term, and (iii) the angular dependence is explicitly included as a truncated Fourier series. Thirteen parameters are used for each metal surface and were fitted on 250 water adsorption energies computed at the PBE+dDsC level. The performance of GAL19 was evaluated on a set of more than 600 DFT adsorption energies for each surface, leading to an average root-mean-square deviation of only 1 kcal/mol, correctly reproducing the adsorption trends: strong on Pt and Pd but weaker on Ag, Au, and Cu. This force field was then used to simulate the water/metal interface for all ten surfaces for 1 ns. Structural analyses reveal similar tendencies for all surfaces: a first, dense water layer that is mostly adsorbed on the metal top sites and a second layer up to around 6 Å, which is less structured. On Pt and Pd, the first layer is strongly organized with water lying flat on the surface. The pairwise additive functional form allows one to simulate the water adsorption on alloys, which is demonstrated at the example of Ag/Cu and Au/Pt alloys. The water/Ag-Cu interface is predicted to be disordered with water mostly adsorbed on Cu which should exacerbate the Ag reactivity. On the contrary, incorporating Pt into Au materials leads to a structuring of the water interface. Our promising results make GAL19 an ideal candidate to get representative sampling of complex metal/water interfaces as a first step toward accurate estimation of free energies of reactions in solution at the metal interface.

Nuclear Quantum Effects from the Analysis of Smoothed Trajectories: Pilot Study for Water

Nuclear quantum effects have significant contributions to thermodynamic quantities and structural properties; furthermore, very expensive methods are necessary for their accurate computation. In most calculations, these effects, for instance, zero-point energies, are simply neglected or only taken into account within the quantum harmonic oscillator approximation. Herein, we present a new method, Generalized Smoothed Trajectory Analysis, to determine nuclear quantum effects from molecular dynamics simulations. The broad applicability is demonstrated with the examples of a harmonic oscillator and different states of water. Ab initio molecular dynamics simulations have been performed for ideal gas up to the temperature of 5000 K. Classical molecular dynamics have been carried out for hexagonal ice, liquid water, and vapor at atmospheric pressure. With respect to the experimental heat capacity, our method outperforms previous calculations in the literature in a wide temperature range at lower computational cost than other alternatives. Dynamic and structural nuclear quantum effects of water are also discussed.

Anomalies in Supercooled Water at ∼230 K Arise from a 1D Polymer to 2D Network Topological Transformation

Puzzling anomalous properties of water are drastically enhanced in the supercooled region. However, the nature of these anomalies is not known. We report here molecular dynamics simulations using the RexPoN force field from 298 to 200 K along the 1 atm density curve. At 298 K, there are 2.1 strong hydrogen bonds (SHBs), leading to a dynamic branched one-dimensional (1D) polymer. Water remains 1D down to 240 K, but at and below 230 K, the number of SHBs becomes 3.0, leading to a two-dimensional (2D) network that persists to 200 K. We propose that this 1D-to-2D topological transition accounts for the anomalous properties of supercooled water. Near 230 K, the power spectra show dramatic increases in the angular vibrational frequency modes, while the diffusivity decreases dramatically, both arising from the 1D-to-2D transformation. This transition is not first order because free energy changes uniformly but fluctuations in the entropy near 230 K suggest a possible second-order transition.

Liquid water is a dynamic polydisperse branched polymer

We developed the RexPoN force field for water based entirely on quantum mechanics. It predicts the properties of water extremely accurately, with T _melt = 273.3 K (273.15 K) and properties at 298 K: ΔH_vap = 10.36 kcal/mol (10.52), density = 0.9965 g/cm³ (0.9965), entropy = 68.4 J/mol/K (69.9), and dielectric constant = 76.1 (78.4), where experimental values are in parentheses. Upon heating from 0.0 K (ice) to 273.0 K (still ice), the average number of strong hydrogen bonds (SHBs, r_OO ≤ 2.93 Å) decreases from 4.0 to 3.3, but upon melting at 273.5 K, the number of SHBs drops suddenly to 2.3, decreasing slowly to 2.1 at 298 K and 1.6 at 400 K. The lifetime of the SHBs is 90.3 fs at 298 K, increasing monotonically for lower temperature. These SHBs connect to form multibranched polymer chains (151 H₂O per chain at 298 K), where branch points have 3 SHBs and termination points have 1 SHB. This dynamic fluctuating branched polymer view of water provides a dramatically modified paradigm for understanding the properties of water. It may explain the 20-nm angular correlation lengths at 298 K and the critical point at 227 K in supercooled water. Indeed, the 15% jump in the SHB lifetime at 227 K suggests that the supercooled critical point may correspond to a phase transition temperature of the dynamic polymer structure. This paradigm for water could have a significant impact on the properties for protein, DNA, and other materials in aqueous media.