Improving pKa Predictions with Reparameterized Force Fields and Free Energy Calculations

Given the growing interest in designing targeted covalent inhibitors, methods for rapidly and accurately probing pKas─and, by extension, the reactivities─of target cysteines are highly desirable. Complementary to cysteine, histidine is similarly relevant due to its frequent presence in protein active sites and its unique ability to exist in two tautomeric states. Here, we demonstrate that nonequilibrium free energy calculations can accurately determine the pKa values of both residues, often outperforming conventional predictors. Importantly, we find that (1) increasing the van der Waals radius of cysteine's sulfur atom, (2) modifying the backbone charges of histidine, and (3) introducing effective polarization by downscaling the side chain partial charges of both residues can all significantly improve pKa prediction accuracy. Using the modified CHARMM36m force field on the full dataset reduces the prediction error from 2.12 ± 0.27 pK to 1.28 ± 0.15 pK and increases the correlation with experiment from 0.25 ± 0.09 to 0.58 ± 0.08. Similarly, using the modified Amber14SB force field decreases the error from 3.21 ± 0.29 pK to 1.69 ± 0.23 pK and improves the correlation from 0.15 ± 0.10 to 0.36 ± 0.10.

Evaluating the Chemical Reactivity of DFT-Simulated Liquid Water with Hydrated Electrons via the Dual Descriptor

Modeling the various properties of liquid water, particularly its reactivity, has been a longstanding challenge for simulation methods. Recently, ab initio simulations based on density functional theory (DFT) have come to the fore as tenable methods for calculating the properties and reactivity of water, with varying degrees of success for different exchange-correlation functionals. In particular, hybrid-GGA and meta-GGA functionals have been shown to reproduce many of the structural, dynamical, and energetic properties of water to a high degree of accuracy relative to their computational cost. Here, we show that the dual descriptor (DD) measure of nucleophilicity and electrophilicity, which is sometimes used to elucidate organic chemistry reaction mechanisms, can also be used to characterize the reactivity of DFT-simulated liquid water. The DD is especially apt for understanding the reactivity of excess electrons with water as its calculation explicitly involves adding and removing an excess electron from a reference system. We use the DD to explore the reactivity of water simulated using three different DFT functionals: the LDA functional (LDA), a hybrid-GGA functional (PBE0), and a hybrid meta-GGA functional (SCAN0). Using the DD, we show that the SCAN0 functional with the standard 25% Hartree-Fock exchange produces simulated liquid water with many regions that are far more reactive than either PBE0 or LDA. To understand the implications of these highly reactive regions, we then add a strong nucleophile in the form of an excess electron and find that although PBE0 and LDA predict stable hydrated electrons, the excess electron reacts nearly instantaneously with SCAN0 water via proton abstraction to form a hydrogen atom and hydroxide ion. We show that the DD provides the ability to not only predict whether or not liquid water will react with a hydrated electron but also which particular waters will be involved solely from analyzing pure water configurations generated with each functional. We rationalize this result in terms of the known trap-seeking behavior of injected hydrated electrons, which are able to find the most electronegative region in bulk water. These results highlight the utility of the dual descriptor as a fast and interpretable method for investigating condensed-phase reactivity with excess electrons.

Crumbling crystals: on the dissolution mechanism of NaCl in water

Dissolution of ionic salts in water is ubiquitous, particularly for NaCl. However, an atomistic scale understanding of the process remains elusive. Simulations lend themselves conveniently to studying dissolution since they provide the spatio-temporal resolution that can be difficult to obtain experimentally. Nevertheless, the complexity of various inter- and intra-molecular interactions require careful treatment and long time scale simulations, both of which are typically hindered by computational expense. Here, we use advances in machine learning potential methodology to resolve at an ab initio level of theory the dissolution mechanism of NaCl in water. The picture that emerges is that of a steady ion-wise unwrapping of the crystal preceding its rapid disintegration, reminiscent of crumbling. The onset of crumbling can be explained by a strong increase in the ratio of the surface area to volume of the crystal. Overall, dissolution comprises a series of highly dynamical microscopic sub-processes, resulting in an inherently stochastic mechanism. These atomistic level insights contribute to the general understanding of dissolution mechanisms in other crystals, and the methodology is primed for more complex systems of recent interest such as water/salt interfaces under flow and salt crystals under confinement.

Structures of Th4+ aqueous solutions: insights from AIMD and metadynamics simulations

Solution chemistry of actinide ions is critical to understanding the solvation behaviors and hydrolysis process. Using tetravalent thorium ion Th4+ as a representative example, we investigate the local structures and dynamic behaviors of hydrated Th4+ ions by ab initio molecular dynamics (AIMD) simulations using the recently developed norm-conserving pseudopotentials and basis sets optimized for actinides (J.-B. Lu et al., J. Chem. Theory Comput. 2021, 17, 3360-3371). AIMD simulations reveal two distinct solvation shells, with the first shell comprising 9 water molecules at approximately rTh-O = 2.50 Å and exhibiting a tricapped trigonal prism geometry. These conclusions are confirmed through metadynamics simulations and further structural analysis. AIMD simulations also show the slight effect of temperature and counterions on the structure of the solution. The structured solvation shells of the highly charged Th4+ ion with the specific geometry, distinct from the structure of liquid water, lead to corresponding structural changes in the hydrogen bond network in water. Additionally, beyond the solvent-shared ion pair (SIP) state observed in the unbiased AIMD simulations, the metadynamics simulations reconstruct a two-dimensional free energy surface that clearly indicates the potential stability of the contact ion pair (CIP) state in the system with Cl- as a counterion. The findings in this work provide insights into the solution chemistry of actinides and serve as a reference for studying other actinide solution systems.

Ab Initio Simulation of Liquid Water without Artificial High Temperature

Comprehending the structure and dynamics of water is crucial in various fields, such as water desalination, ion separation, electrocatalysis, and biochemical processes. While reported works show that the ab initio molecular dynamics (AIMD) can accurately portray water's structure, the artificial high temperature (AHT) from 120 to 30 K is needed to mimic the quantum nature of hydrogen-bond network from GGA, metaGGA to hybrid functionals. The AHT proves to be an inadequate approach for systems involving aqueous multiphase mixtures, such as water-solid interfaces and aqueous solutions. This is due to the activation of additional phonons in other phases, which can lead to an overestimation of the dynamics of nearby water molecules. In this work, we find that the regularized SCAN (rSCAN) functional effectively captures both the structure and dynamics of liquid water at ambient conditions without AHT. Moreover, rSCAN closely matches experimental results for the hydration structures of alkali, alkali earth, and halide ions. We anticipate that the versatile and accurate rSCAN functional will emerge as a key tool based on ab initio simulation for investigating chemical processes in aqueous environments.

To Pair or not to Pair? Machine-Learned Explicitly-Correlated Electronic Structure for NaCl in Water

The extent of ion pairing in solution is an important phenomenon to rationalize transport and thermodynamic properties of electrolytes. A fundamental measure of this pairing is the potential of mean force (PMF) between solvated ions. The relative stabilities of the paired and solvent shared states in the PMF and the barrier between them are highly sensitive to the underlying potential energy surface. However, direct application of accurate electronic structure methods is challenging, since long simulations are required. We develop wave function based machine learning potentials with the random phase approximation (RPA) and second order Møller-Plesset (MP2) perturbation theory for the prototypical system of Na and Cl ions in water. We show both methods in agreement, predicting the paired and solvent shared states to have similar energies (within 0.2 kcal/mol). We also provide the same benchmarks for different DFT functionals as well as insight into the PMF based on simple analyses of the interactions in the system.

Characterization of the Coordination and Solvation Dynamics of Solvated Systems─Implications for the Analysis of Molecular Interactions in Solutions and Pure H2O

The characterization of solvation shells of atoms, ions, and molecules in solution is essential to relate solvation properties to chemical phenomena such as complex formation and reactivity. Different definitions of the first-shell coordination sphere from simulation data can lead to potentially conflicting data on the structural properties and associated ligand exchange dynamics. The definition of a solvation shell is typically based on a given threshold distance determined from the respective solute-solvent pair distribution function g(r) (i.e., GC). Alternatively, a nearest neighbor (NN) assignment based on geometric properties of the coordination complex without the need for a predetermined cutoff criterion, such as the relative angular distance (RAD) or the modified Voronoi (MV) tessellation, can be applied. In this study, the effect of different NN algorithms on the coordination number and ligand exchange dynamics evaluated for a series of monatomic ions in aqueous solution, carbon dioxide in aqueous and dichloromethane solutions, and pure liquid water has been investigated. In the case of the monatomic ions, the RAD approach is superior in achieving a well separated definition of the first solvation layer. In contrast, the MV algorithm provides a better separation of the NNs from a molecular point of view, leading to better results in the case of solvated CO2. When analyzing the coordination environment in pure water, the cutoff-based GC framework was found to be the most reliable approach. By comparison of the number of ligand exchange reactions and the associated mean ligand residence times (MRTs) with the properties of the coordination number autocorrelation functions, it is shown that although the average coordination numbers are sensitive to the different definitions of the first solvation shell, highly consistent estimates for the associated MRT of the solvated system are obtained in the majority of cases.

Consistent density functional theory-based description of ion hydration through density-corrected many-body representations

Delocalization error constrains the accuracy of density functional theory in describing molecular interactions in ion-water systems. Using Na+ and Cl- in water as model systems, we calculate the effects of delocalization error in the SCAN functional for describing ion-water and water-water interactions in hydrated ions, and demonstrate that density-corrected SCAN (DC-SCAN) predicts n-body and interaction energies with an accuracy approaching coupled cluster theory. The performance of DC-SCAN is size-consistent, maintaining an accurate description of molecular interactions well beyond the first solvation shell. Molecular dynamics simulations at ambient conditions with many-body MB-SCAN(DC) potentials, derived from the many-body expansion, predict the solvation structure of Na+ and Cl- in quantitative agreement with reference data, while simultaneously reproducing the structure of liquid water. Beyond rationalizing the accuracy of density-corrected models of ion hydration, our findings suggest that our unified density-corrected MB formalism holds great promise for efficient DFT-based simulations of condensed-phase systems with chemical accuracy.

Solvent Dynamics of Aqueous Halides before and after Photoionization

Electron transfer reactions can be strongly influenced by solvent dynamics. We study the photoionization of halides in water as a model system for such reactions. There are no internal nuclear degrees of freedom in the solute, allowing the dynamics of the solvent to be uniquely identified. We simulate the equilibrium solvent dynamics for Cl^-, Br^-, I^-, and their respective neutral atoms in water, comparing quantum mechanical/molecular mechanical (QM/MM) and classical molecular dynamics (MD) methods. On the basis of the obtained configurations, we calculate the extended X-ray absorption fine structure (EXAFS) spectra rigorously based on the MD snapshots and compare them in detail with other theoretical and experimental results available in the literature. We find our EXAFS spectra based on QM/MM MD simulations in good agreement with their experimental counterparts for the ions. Classical MD simulations for the ions lead to EXAFS spectra that agree equally well with the experiment when it comes to the oscillatory period of the signal, even though they differ from the QM/MM radial distribution functions extracted from the MD. The amplitude is, however, considerably overestimated. This suggests that to judge the reliability of theoretical simulation methods or to elucidate fine details of the atomistic dynamics of the solvent based on EXAFS spectra, the amplitude as well as the oscillatory period need to be considered. If simulations fail qualitatively, as does the classical MD for the aqueous neutral halogen atoms, the resulting EXAFS will also be strongly affected in both oscillatory period and amplitude. The good reliability of QM/MM-based EXAFS simulations, together with clear qualitative differences in the EXAFS spectra found between halides and their atomic counterparts, suggests that a combined theory and experimental EXAFS approach is suitable for elucidating the nonequilibrium solvent dynamics in the photoionization of halides and possibly also for electron transfer reactions in more complex systems.

Quantifying the Structure of Water and Hydrated Monovalent Ions by Density Functional Theory-Based Molecular Dynamics

The accurate description of the structures of water and hydrated ions is important in electrochemical desalination, ion separation, and supercapacitors. In this work, we present an ab initio atomistic simulation-based study to explore the structure of water and hydrated monovalent ions (Li⁺, Na⁺, K⁺, Rb⁺, F^-, and Cl^-) at ambient conditions using generalized gradient approximation (GGA)-based methods with and without van der Waals correction (PBE, PBE + D3, and revPBE + D3) and recently developed strongly constrained and appropriately normed (SCAN) meta-GGA. We find that both revPBE + D3 and SCAN can well capture the structure of bulk water with +30 K artificial high temperature in contrast to overstructuring water using PBE and PBE + D3. However, being the same as PBE + D3, revPBE + D3 overestimates the structure of the hydration shell, especially for monovalent cations. Surprisingly, SCAN can well match the experimental results of hydrated monovalent ions. Detailed structure analyzes of entropy reveal that the hydration shell under the level of PBE + D3 and revPBE + D3 is more disordered and looser than SCAN. The successful prediction of the flexible SCAN functional could facilitate the exploration of complex ionic processes in the aqueous phase, the interactions of hydrated ions with surfaces, and solvation states in nanopores at an accurate, efficient, predictive, and ab initio level.

The Ion-Dipole Correction of the 3DRISM Solvation Model to Accurately Compute Water Distributions around Negatively Charged Biomolecules

The 3D reference interaction site model (3DRISM) provides an efficient grid-based solvation model to compute the structural and thermodynamic properties of biomolecules in aqueous solutions. However, it remains challenging for existing 3DRISM methods to correctly predict water distributions around negatively charged solute molecules. In this paper, we first show that this challenge is mainly due to the orientation of water molecules in the first solvation shell of the negatively charged solute molecules. To properly consider this orientational preference, position-dependent two-body intramolecular correlations of solvent need to be included in the 3DRISM theory, but direct evaluations of these position-dependent two-body intramolecular correlations remain numerically intractable. To address this challenge, we introduce the Ion-Dipole Correction (IDC) to the 3DRISM theory, in which we incorporate the orientation preference of water molecules via an additional solute-solvent interaction term (i.e., the ion-dipole interaction) while keeping the formulism of the 3DRISM equation unchanged. We prove that this newly introduced IDC term is equivalent to an effective direct correlation function which can effectively consider the orientation effect that arises from position dependent two-body correlations. We first quantitatively validate our 3DRISM-IDC theory combined with the PSE3 closure on Cl^-, [ClO]^- (a two-site anion), and [NO₂]^- (a three-site anion). For all three anions, we show that our 3DRISM-IDC theory significantly outperforms the 3DRISM theory in accurately predicting the solvation structures in comparison to MD simulations, including RDFs and 3D water distributions. Furthermore, we have also demonstrated that the 3DRISM-IDC can improve the accuracy of hydration free-energy calculation for Cl^-. We further demonstrate that our 3DRISM-IDC theory yields significant improvements over the 3DRISM theory when applied to compute the solvation structures for various negatively charged solute molecules, including adenosine triphosphate (ATP), a short peptide containing 19 residues, a DNA hairpin containing 24 nucleotides, and a riboswitch RNA molecule with 77 nucleotides. We expect that our 3DRISM-IDC-PSE3 solvation model holds great promise to be widely applied to study solvation properties for nucleic acids and other biomolecules containing negatively charged functional groups.

Extended X-ray absorption fine structure spectroscopy measurements and ab initio molecular dynamics simulations reveal the hydration structure of the radium(II) ion

Radium is refocused from the viewpoint of an environmental pollutant and cancer therapy using alpha particles, where it mainly exists as a hydrated ion. We investigated the radium hydration structure and the dynamics of water molecules by extended X-ray absorption fine structure (EXAFS) spectroscopy and ab initio molecular dynamics (AIMD) simulation. The EXAFS experiment showed that the coordination number and average distance between radium ion and the oxygen atoms in the first hydration shell are 9.2 ± 1.9 and 2.87 ± 0.06 Å, respectively. They are consistent with those obtained from the AIMD simulations, 8.4 and 2.88 Å. The AIMD simulations also revealed that the water molecules in the first hydration shell of radium are less structured and more mobile than those of barium, which is an analogous element of radium. Our results indicate that radium can be more labile than barium in terms of interactions with water.

Hydration in aqueous NaCl

Atomistic details about the hydration of ions in aqueous solutions are still debated due to the disordered and statistical nature of the hydration process. However, many processes from biology, physical chemistry to materials sciences rely on the complex interplay between solute and solvent. Oxygen K-edge X-ray excitation spectra provide a sensitive probe of the local atomic and electronic surrounding of the excited sites. We used ab initio molecular dynamics simulations together with extensive spectrum calculations to relate the features found in experimental oxygen K-edge spectra of a concentration series of aqueous NaCl with the induced structural changes upon solvation of the salt and distill the spectral fingerprints of the first hydration shells around the Na⁺- and Cl^--ions. By this combined experimental and theoretical approach, we find the strongest spectral changes to indeed result from the first hydration shells of both ions and relate the observed shift of spectral weight from the post- to the main-edge to the origin of the post-edge as a shape resonance.

Automated Workflow for Computation of Redox Potentials, Acidity Constants and Solvation Free Energies Accelerated by Machine Learning

Fast evolution of modern society stimulates intense development of new materials with novel functionalities in energy and environmental applications. Due to rapid progress of computer science, computational design of materials with target properties has recently attracted lots of interests. Accurate and efficient calculation of fundamental thermodynamic properties, including redox potentials, acidity constants, and solvation free energies, is of great importance for selection and design of desirable materials. Free energy calculation based on ab initio molecular dynamics (AIMD) can predict these properties with high accuracy at complex environments, however being impeded by high computational costs. To address this issue, this work develops an automated scheme that combines iterative training of machine learning potentials (MLPs) and free energy calculation, and deomonstrates that these thermodynamic properties can be computed by ML accelerated MD with ab initio accuracy and much longer time scale at cheaper costs, improving poor statistics and convergence of numerical integration by AIMD. Our automated scheme lays the foundation for computational chemistry-assisted materials design.

Core and Valence Level Photoelectron Spectroscopy of Nanosolvated KCl

The solvation of alkali and halide ions in the aqueous environment has been a subject of intense experimental and theoretical research with multidisciplinary interests; yet, a comprehensive molecular-level understanding has still not been obtained. In recent years, electron spectroscopy has been increasingly applied to study the electronic and structural properties of aqueous ions with implications, especially in atmospheric chemistry. In this work, we report core and valence level (Cl 2p, Cl 3p, and K 3p) photoelectron spectra of the common alkali halide, KCl, doped in gas-phase water clusters in the size range of a few hundred water molecules. The results indicate that the electronic structure of these nanosolutions shows a distinct character from that observed at the liquid-vapor interface in liquid microjets and ambient pressure setups. Insights are provided into the unique solvation properties of ions in a nanoaqueous environment, emerging properties of bulk electrolyte solutions with growing cluster size, and sensitivity of the electronic structure to varying solvation configurations.

Role of water model on ion dissociation at ambient conditions

We study ion pair dissociation in water at ambient conditions using a combination of classical and ab initio approaches. The goal of this study is to disentangle the sources of discrepancy observed in computed potentials of mean force. In particular, we aim to understand why some models favor the stability of solvent-separated ion pairs vs contact ion pairs. We found that some observed differences can be explained by non-converged simulation parameters. However, we also unveil that for some models, small changes in the solution density can have significant effects on modifying the equilibrium balance between the two configurations. We conclude that the thermodynamic stability of contact and solvent-separated ion pairs is very sensitive to the dielectric properties of the underlying simulation model. In general, classical models are very robust in providing a similar estimation of the contact ion pair stability, while this is much more variable in density functional theory-based models. The barrier to transition from the solvent-separated to contact ion pair is fundamentally dependent on the balance between electrostatic potential energy and entropy. This reflects the importance of water intra- and inter-molecular polarizability in obtaining an accurate description of the screened ion-ion interactions.

Failure of molecular dynamics to provide appropriate structures for quantum mechanical description of the aqueous chloride ion charge-transfer-to-solvent ultraviolet spectrum

The lowest band in the charge-transfer-to-solvent ultraviolet absorption spectrum of aqueous chloride ion is studied by experiment and computation. Interestingly, the experiments indicate that at concentrations up to at least 0.25 M, where calculations indicate ion pairing to be significant, there is no notable effect of ionic strength on the spectrum. The experimental spectra are fitted to aid comparison with computations. Classical molecular dynamic simulations are carried out on dilute aqueous Cl-, Na+, and NaCl, producing radial distribution functions in reasonable agreement with experiment and, for NaCl, clear evidence of ion pairing. Clusters are extracted from the simulations for quantum mechanical excited state calculations. Accurate ab initio coupled-cluster benchmark calculations on a small number of representative clusters are carried out and used to identify and validate an efficient protocol based on time-dependent density functional theory. The latter is used to carry out quantum mechanical calculations on thousands of clusters. The resulting computed spectrum is in excellent agreement with experiment for the peak position, with little influence from ion pairing, but is in qualitative disagreement on the width, being only about half as wide. It is concluded that simulation by classical molecular dynamics fails to provide an adequate variety of structures to explain the experimental CTTS spectrum of aqueous Cl-.