Absolute proton hydration free energy, surface potential of water, and redox potential of the hydrogen electrode from first principles: QM/MM MD free-energy simulations of sodium and potassium hydration

The Determination of Free Energy of Hydration of Water Ions from First Principles

We model the autoionization of water by determining the free energy of hydration of the major intermediate species of water ions. We represent the smallest ions─the hydroxide ion OH-, the hydronium ion H3O+, and the Zundel ion H5O2+─by bonded models and the more extended ionic structures by strong nonbonded interactions (e.g., the Eigen H9O4+ = H3O+ + 3(H2O) and the Stoyanov H13O6+ = H5O2+ + 4(H2O)). Our models are faithful to the precise QM energies and their components to within 1% or less. Using the calculated free energies and atomization energies, we compute the pKa of pure water from first principles as a consistency check and arrive at a value within 1.3 log units of the experimental one. From these calculations, we conclude that the hydronium ion, and its hydrated state, the Eigen cation, are the dominant species in the water autoionization process.

Reshaping the QM Region On-the-Fly: Adaptive-Shape QM/MM Dynamic Simulations of a Hydrated Proton in Bulk Water

Adaptive quantum mechanics/molecular mechanics (QM/MM) reclassifies on-the-fly a molecule or molecular fragment as QM or MM during dynamics simulations without abrupt changes in the energy or forces. Notably, the permuted adaptive-partitioning (PAP) algorithms have been applied to simulate a hydrated proton, with a mobile QM zone anchored at a pseudoatom called a proton indicator. The position of the proton indicator approximates the location of the delocalized excess proton, yielding a smooth trajectory of the proton diffusing via the Grotthuss mechanism in aqueous solutions. The mobile QM zone, which has been taken to be a sphere with a preset radius, follows the proton wherever it goes. Although the simulations are successful, the use of a spherical QM zone has one disadvantage: A large preset radius must be utilized to minimize the chance of missing water molecules that are important to proton translocation. A large radius leads to a large QM zone, which is computationally expensive. In this work, we report a new way to set up the QM zone, where one includes only the water molecules important to proton transfer. The importance of a given water molecule is quantified by its "weight" that depends on its relation to the reaction path of proton transfer. The weight varies smoothly, ensuring that a water molecule gradually appears in or disappears from the QM zone without abrupt changes, as required by the PAP method. Consequently, the shape of the QM zone evolves on-the-fly, keeping the QM zone as small as possible and as large as necessary. Test simulations demonstrate that the new algorithm significantly improves the computation efficiency while maintaining the proper descriptions of proton transfer in bulk water.

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.

Experimental Compilation and Computation of Hydration Free Energies for Ionic Solutes

Although charged solutes are common in many chemical systems, traditional solvation models perform poorly in calculating solvation energies of ions. One major obstacle is the scarcity of experimental data for solvated ions. In this study, we release an experiment-based aqueous ionic solvation energy data set, IonSolv-Aq, that contains hydration free energies for 118 anions and 155 cations, more than 2 times larger than the set of hydration free energies for singly charged ions contained in the 2012 Minnesota Solvation Database commonly used in benchmarking studies. We discuss sources of systematic uncertainty in the data set and use the data to examine the accuracy of popular implicit solvation models COSMO-RS and SMD for predicting solvation free energies of singly charged ionic solutes in water. Our results indicate that most SMD and COSMO-RS modeling errors for ionic solutes are systematic and correctable with empirical parameters. We discuss two systematic offsets: one across all ions and one that depends on the functional group of the ionization site. After correcting for these offsets, solvation energies of singly charged ions are predicted using COSMO-RS to 3.1 kcal mol-1 MAE against a challenging test set and 1.7 kcal mol-1 MAE (about 3% relative error) with a filtered test set. The performance of SMD is similar, with MAE against those same test sets of 2.7 and 1.7 kcal mol-1. These results underscore the importance of compiling larger experimental data sets to improve solvation model parametrization and fairly assess performance.

Adaptive-Partitioning Multilayer Dynamics Simulations: 2. Implementations of the Permuted and Interpolated Adaptive-Partitioning Gradients

Recently, an adaptive-partitioning multilayer Q1/Q2/MM method was proposed, where Q1 and Q2 denote, respectively, two distinct quantum-mechanical levels of theory and MM, the molecular-mechanical force fields. Such a multilayer model resembles the ONIOM (our own N-layered integrated molecular orbital and molecular mechanics) model by Morokuma and co-workers, but it is distinguished by on-the-fly reclassifying atoms to be Q1, Q2, or MM in dynamics simulations. To smoothly blend the levels of descriptions of the atoms, buffer zones are introduced between adjacent layers, and the energy is smoothly interpolated. In particular, the Q1/Q2 interaction energy was expressed in two different formalisms: permuted and interpolated adaptive-partitioning (PAP and IAP), respectively. While the PAP energy is based on a weighted many-body expansion, the IAP energy is derived via alchemical quantum calculations with interpolated Fock and overlap matrices. In this article, we examine in-depth the irregularities in the IAP energy near the boundary between the buffer and Q2 zones, which were found prominent in some calculations. These irregularities are due to basis-set linear dependencies, which can be effectively suppressed using a cutoff for the weighted atomic orbital coefficients. Furthermore, we derived and implemented the gradients for both PAP and IAP. Test calculations on a series of water cluster models show perfectly smooth gradients in PAP, while a minor discontinuity occurs in IAP gradients at the buffer/Q2 boundary. The energy and gradient discontinuities in IAP become smaller when moving the buffer/Q2 boundary further away from the Q1 center and when increasing the size of the basis sets used. Overall, those discontinuities are controllable, and possible ways to further diminish them are discussed.

Charged species redistribution at electrochemical interfaces: a model system of the zirconium oxide/water interface

Quantifying the local distribution of charged defects in the solid state and charged ions in liquid solution near the oxide/liquid interface is key to understanding a range of important electrochemical processes, including oxygen reduction and evolution, corrosion and hydrogen evolution reactions. Based on a grand canonical approach relying on the electrochemical potential of individual charged species, a unified treatment of charged defects on the solid side and ions on the water side can be established. This approach is compatible with first-principles calculations where the formation free energy of individual charged species can be calculated and modulated by imposing certain electrochemical potential. Herein, we apply this framework to a system of monoclinic ZrO₂(1̄11)/water interface. The structure, defect chemistry and dynamical behavior of the electric double layer and space charge layer are analyzed with different pH values, water chemistry and doping elements in zirconium oxide. The model predicts ZrO₂ solubility in water and the point of zero charge consistent with the experimentally-measured values. We reveal the effect of dopant elements on the concentrations of oxygen and hydrogen species at the surface of the ZrO₂ passive layer in contact with water, uncovering an intrinsic trade-off between oxygen diffusion and hydrogen pickup during the corrosion of zirconium alloys. The solid/water interface model established here serves as the basis for modeling reaction and transport kinetics under doping and water chemistry effects.

Solvation of Manganese(III) Ion in Water and in Ammonia

In this work, we have studied the solvation of manganese(III) ion in water and in ammonia using three levels of theory: MP2, MN15, and ωB97XD associated with the aug-cc-pVDZ basis set. The studied systems are constituted of Mn³⁺(H₂O)₆ and Mn³⁺(NH₃)₆ in gas and solvent phases as well as Mn³⁺(H₂O)₁₈ and Mn³⁺(NH₃)₁₈ in the gas phase. Four aspects of the solvation of manganese(III) ion have been examined for the aforementioned systems at the three levels of theory. First, we started by locating the Jahn-Teller elongated and compressed configuration in Mn³⁺(H₂O)₆ and Mn³⁺(NH₃)₆. Second, we calculated the spin state energies and the spin state free energies for temperatures ranging from 50 to 400 K to look at possible spin crossover in the studied systems. Third, we carried out a quantum theory of atoms in molecules (QTAIM) analysis, and we determined the ionic radii of manganese(III) ion in water and in ammonia. Fourth, we calculated the solvation free energies and the solvation enthalpies of manganese(III) ion in water and in ammonia using the cluster continuum solvation model. For these four aspects of the solvation of manganese(III) ion, most of the reported properties are provided in this work for the first time. We particularly found that the calculated solvation enthalpy of the manganese(III) ion in water is in good agreement with an experimental estimate.

Tracking the Delocalized Proton in Concerted Proton Transfer in Bulk Water

A solvated proton in water is often characterized as a charge or structural defect, and it is important to track its evolution on-the-fly in certain dynamics simulations. Previously, we introduced the proton indicator, a pseudo-atom, whose position approximates the location of the excess proton modeled as a structural defect. The proton indicator generally yields a smooth trajectory of a hydrated proton diffusing in aqueous solutions, including in the events of stepwise proton transfer via the Grotthuss mechanism; however, the proton indicator did not perform well in the events of concerted proton transfer, for which it occasionally yielded large position displacements between two successive time steps. To overcome this hurdle, we develop a new algorithm of a proton indicator with greatly enhanced performance for concerted proton transfer in bulk water. A protocol is proposed to exhaustively explore the hydrogen-bonding network of the water wires over which the excess proton is delocalized and to properly account for the contributions of the water molecules in this network as the geometry evolves. The new proton indicator (called Indicator 2.0) is assessed in dynamics simulations of an excess proton in bulk water and in specially constructed model systems of more complex architectures. The results demonstrate that the new indicator yields a smooth trajectory in both stepwise and concerted proton transfers.

Electrokinetic, electrochemical, and electrostatic surface potentials of the pristine water liquid-vapor interface

Although conceptually simple, the air-water interface displays rich behavior and is subject to intense experimental and theoretical investigations. Different definitions of the electrostatic surface potential as well as different calculation methods, each relevant for distinct experimental scenarios, lead to widely varying potential magnitudes and sometimes even different signs. Based on quantum-chemical density-functional-theory molecular dynamics (DFT-MD) simulations, different surface potentials are evaluated and compared to force-field (FF) MD simulations. As well explained in the literature, the laterally averaged electrostatic surface potential, accessible to electron holography, is dominated by the trace of the water molecular quadrupole moment, and using DFT-MD amounts to +4.35 V inside the water phase, very different from results obtained with FF water models which yield negative values of the order of -0.4 to -0.6 V. Thus, when predicting potentials within water molecules, as relevant for photoelectron spectroscopy and non-linear interface-specific spectroscopy, DFT simulations should be used. The electrochemical surface potential, relevant for ion transfer reactions and ion surface adsorption, is much smaller, less than 200 mV in magnitude, and depends specifically on the ion radius. Charge transfer between interfacial water molecules leads to a sizable surface potential as well. However, when probing electrokinetics by explicitly applying a lateral electric field in DFT-MD simulations, the electrokinetic ζ-potential turns out to be negligible, in agreement with predictions using continuous hydrodynamic models. Thus, interfacial polarization charges from intermolecular charge transfer do not lead to significant electrokinetic mobility at the pristine vapor-liquid water interface, even assuming these transfer charges are mobile in an external electric field.

Physical and numerical aspects of sodium ion solvation free energies via the cluster-continuum model

Sodium cation solvation Gibbs free energies (ΔG_solv(Na⁺)) have been obtained in water, dimethylformamide, dimethyl sulfoxide, ethanol, acetone, acetonitrile, and methanol through the "monomer cycle" cluster-continuum approach where a solvent reference state is described by infinitely separated molecules. The following steps are vital for obtaining reliable ΔG_solv(Na⁺) values: (a) a meticulous conformational search involving dispersion corrected density functional theory (DFT-D) and the continuum solvation model (CSM); (b) gas-phase DFT-D geometry optimization followed by single-point (SP) domain-based local pair natural orbital coupled clusters including single, double, and partly triple excitation (DLPNO-CCSD(T)) calculations in conjunction with the complete basis set extrapolation; (c) advanced statistical thermodynamic treatment of the low harmonic frequencies (<100 cm^-1) to obtain the robust gas-phase Gibbs free energy correction; (d) gas-phase and dielectric continuum SP with non-electrostatic contributions included in the CSM; (e) an evaluation of the relative thermodynamic stability of the Na⁺(S)_n clusters to identify the number of explicit solvent molecules n to be considered. Our refined computational protocol is promising with a Pearson correlation coefficient between the predicted and experimental data, ρ, of 0.82, and the mean signed and mean unsigned errors of 0.3 and 1.4 kcal mol^-1, respectively.

Hydrogen Bond Thermodynamics in Aqueous Acid Solutions: A Combined DFT and Classical Force-Field Approach

The thermodynamics of hydrogen bonds in aqueous and acidic solutions significantly impacts the kinetics and thermodynamics of acid reaction chemistry. We utilize in this work a multiscale approach, combining density functional theory (DFT) with classical molecular dynamics (MD) to model hydrogen bond thermodynamics in an acidic solution. Using thermodynamic cycles, we split the solution phase free energy into its gas phase counterpart plus solvation free energies. We validate this DFT/MD approach by calculating the aqueous phase hydrogen bond free energy between two water molecules (H₂O-···-H₂O), the free energy to transform an H₃O⁺ cation into an H₅O₂⁺ cation, and the hydrogen bond free energy of protonated water clusters (H₃O⁺-···-H₂O and H₅O₂⁺-···-H₂O). The computed equilibrium hydrogen bond free energy of H₂O-···-H₂O is remarkably accurate, especially considering the large individual contributions to the thermodynamic cycle. Turning to cations, we find the ion to be more stable than H₃O⁺ by roughly 1-2 k_BT. This small free energy difference allows for thermal fluctuation between the two idealized motifs, consistent with spectroscopic and simulation studies. Lastly, hydrogen bonding free energies between either H⁺ cation and H₂O in solution were found to be stronger than between two H₂O, though much less so than in vacuum because of dielectric screening in solution. Altogether, our results suggest the DFT/MD approach is promising for application in modeling hydrogen bonding and proton transfer thermodynamics in condensed phases.

Measurements and Utilization of Consistent Gibbs Energies of Transfer of Single Ions: Towards a Unified Redox Potential Scale for All Solvents

Utilizing the “ideal” ionic liquid salt bridge to measure Gibbs energies of transfer of silver ions between the solvents water, acetonitrile, propylene carbonate and dimethylformamide results in a consistent data set with a precision of 0.6 kJ mol⁻¹ over 87 measurements in 10 half‐cells. This forms the basis for a coherent experimental thermodynamic framework of ion solvation chemistry. In addition, we define the solvent independent peabsH2O ‐ and the EabsH2O values that account for the electronating potential of any redox system similar to the pHabsH2O value of a medium that accounts for its protonating potential. This EabsH2O scale is thermodynamically well‐defined enabling a straightforward comparison of the redox potentials (reducities) of all media with respect to the aqueous redox potential scale, hence unifying all conventional solvents′ redox potential scales. Thus, using the Gibbs energy of transfer of the silver ion published herein, one can convert and unify all hitherto published redox potentials measured, for example, against ferrocene, to the EabsH2O scale.

Generalizing Continuum Solvation in Crystal to Nonaqueous Solvents: Implementation, Parametrization, and Application to Molecules and Surfaces

We present an extension of a generalized finite-difference Poisson-Boltzmann (FDPB) continuum solvation model based on a self-consistent reaction field treatment to nonaqueous solvents. Implementation and reparametrization of the cavitation, dispersion, and structural (CDS) effects nonelectrostatic model are presented in CRYSTAL, with applications to both finite and infinite periodic systems. For neutral finite systems, computed errors with respect to available experimental data on free energies of solvation of 2523 solutes in 91 solvents, as well as 144 transfer energies from water to 14 organic solvents are on par with the reference SM12 solvation model for which the CDS parameters have been developed. Calculations performed on a TiO₂ anatase surface and compared to VASPsol data revealed an overall very good agreement of computed solvation energies, surface energies, as well as band structure changes upon solvation in three different solvents, validating the general applicability of the reparametrized FDPB approach to neutral nonperiodic and periodic solutes in aqueous and nonaqueous solvents. For ionic species, while the reparametrized CDS model led to large errors on free energies of solvation of anions, addition of a corrective term based on Abraham's acidity of the solvent significantly improved the accuracy of the proposed continuum solvation model, leading to errors on aqueous pK_a of a test set of 83 solutes divided by a factor of 4 compared to the reference solvation model based on density (SMD). Overall, therefore, these encouraging results demonstrate that the generalized FDPB continuum solvation model can be applied to a broad range of solutes in various solvents, ranging from finite neutral or charged solutes to extended periodic surfaces.

Adaptive-Partitioning Multilayer Dynamics Simulations: 1. On-the-Fly Switch between Two Quantum Levels of Theory

We propose to generalize the previously developed two-layer permuted adaptive-partitioning quantum-mechanics/molecular-mechanics (QM/MM), which reclassifies atoms as QM or MM on-the-fly in dynamics simulations, to multilayer adaptive-partitioning algorithms that enable multiple levels of theory. In this work, we formulate two new algorithms that smoothly interpolate the energy between two QM (Q1 and Q2) levels of theory. The first "permuted adaptive-partitioning" scheme is based on the weighted many-body expansion of the potential, as in the adaptive-partitioning QM/MM. Unconventional and potentially more efficient, the second "interpolated adaptive-partitioning" method employs alchemical QM calculations with Q1/Q2-mixed basis sets, Fock matrices, and overlap matrices. To our knowledge, this is the first time that such alchemical calculations are performed in QM, although they are routinely done in MM. Test calculations on water-cluster models show that both new algorithms indeed yield smooth energy curves when water molecules shift between Q1 and Q2.

Toward a First-Principles Framework for Predicting Collective Properties of Electrolytes

Given the universal importance of electrolyte solutions, it is natural to expect that we have a nearly complete understanding of the fundamental properties of these solutions (e.g., the chemical potential) and that we can therefore explain, predict, and control the phenomena occurring in them. In fact, reality falls short of these expectations. But, recent advances in the simulation and modeling of electrolyte solutions indicate that it should soon be possible to make progress toward these goals. In this Account, we will discuss the use of first-principles interaction potentials based in quantum mechanics (QM) to enhance our understanding of electrolyte solutions. Specifically, we will focus on the use of quantum density functional theory (DFT) combined with molecular dynamics simulation (DFT-MD) as the foundation for our approach. The overarching concept is to understand and accurately reproduce the balance between local or short-ranged (SR) structural details and long-range (LR) correlations, allowing the prediction of the thermodynamics of both single ions in solution as well as the collective interactions characterized by activity/osmotic coefficients. In doing so, relevant collective motions and driving forces characterized by chemical potentials can be determined.In this Account, we will make the case that understanding electrolyte solutions requires a faithful QM representation of the SR nature of the ion-ion, ion-water, and water-water interactions. However, the number of molecules that is required for collective behavior makes the direct application of high-level QM methods that contain the best SR physics untenable, making methods that balance accuracy and efficiency a practical goal. Alternatives such as continuum solvent models (CSMs) and empirically based classical molecular dynamics have been extensively employed to resolve this problem but without yet overcoming the fundamental issue of SR accuracy. We will demonstrate that accurately describing the SR interaction is imperative for predicting both intrinsic properties, namely, at infinite dilution, and collective properties of electrolyte solutions.DFT has played an important role in our understanding of condensed phase systems, e.g., bulk liquid water, the air-water interface, ions in bulk, and at the air-water interface. This approach holds huge promise to provide benchmark calculations of electrolyte solution properties that will allow for the development and improvement of more efficient methods, as well as an enhanced understanding of fundamental phenomena. However, the standard protocol using the generalized gradient approximation with van der Waals (vdW) correction requires improvement in order to achieve a high level of quantitative accuracy. Simply simulating with higher level DFT functionals may not be the best route considering the significant computational cost. Alternative methods of incorporating information from higher levels of QM should be explored; e.g., using force matching techniques on small clusters, where high level benchmark calculations are possible, to develop ideal correction terms to the DFT functional is a promising possibility. We argue that DFT with statistical mechanics is becoming an increasingly useful framework enabling the prediction of collective electrolyte properties.

Uncovering Differences in Hydration Free Energies and Structures for Model Compound Mimics of Charged Side Chains of Amino Acids

Free energies of hydration are of fundamental interest for modeling and understanding conformational and phase equilibria of macromolecular solutes in aqueous phases. Of particular relevance to systems such as intrinsically disordered proteins are the free energies of hydration and hydration structures of model compounds that mimic charged side chains of Arg, Lys, Asp, and Glu. Here, we deploy a Thermodynamic Cycle-based Proton Dissociation (TCPD) approach in conjunction with data from direct measurements to obtain estimates for the free energies of hydration for model compounds that mimic the side chains of Arg+, Lys+, Asp-, and Glu-. Irrespective of the choice made for the hydration free energy of the proton, the TCPD approach reveals clear trends regarding the free energies of hydration for Arg+, Lys+, Asp-, and Glu-. These trends include asymmetries between the hydration free energies of acidic (Asp- and Glu-) and basic (Arg+ and Lys+) residues. Further, the TCPD analysis, which relies on a combination of experimental data, shows that the free energy of hydration of Arg+ is less favorable than that of Lys+. We sought a physical explanation for the TCPD-derived trends in free energies of hydration. To this end, we performed temperature-dependent calculations of free energies of hydration and analyzed hydration structures from simulations that use the polarizable Atomic Multipole Optimized Energetics for Biomolecular Applications (AMOEBA) force field and water model. At 298 K, the AMOEBA model generates estimates of free energies of hydration that are consistent with TCPD values with a free energy of hydration for the proton of ca. -259 kcal/mol. Analysis of temperature-dependent simulations leads to a structural explanation for the observed differences in free energies of hydration of ionizable residues and reveals that the heat capacity of hydration is positive for Arg+ and Lys+ and negative for Asp- and Glu-.

Analysis of the Ordering Effects in Anthraquinone Thin Films and Its Potential Application for Sodium Ion Batteries

The ordering effects in anthraquinone (AQ) stacking forced by thin-film application and its influence on dimer solubility and current collector adhesion are investigated. The structural characteristics of AQ and its chemical environment are found to have a substantial influence on its electrochemical performance. Computational investigation for different charged states of AQ on a carbon substrate obtained via basin hopping global minimization provides important insights into the physicochemical thin-film properties. The results reveal the ideal stacking configurations of the individual AQ-carrier systems and show ordering effects in a periodic supercell environment. The latter reveals the transition from intermolecular hydrogen bonding toward the formation of salt bridges between the reduced AQ units and a stabilizing effect upon the dimerlike rearrangement, while the strong surface-molecular interactions in the thin-film geometries are found to be crucial for the formed dimers to remain electronically active. Both characteristics, the improved current collector adhesion and the stabilization due to dimerization, are mutual benefits of thin-film electrodes over powder-based systems. This hypothesis has been further investigated for its potential application in sodium ion batteries. Our results show that AQ thin-film electrodes exhibit significantly better specific capacities (233 vs 87 mAh g-1 in the first cycle), Coulombic efficiencies, and long-term cycling performance (80 vs 4 mAh g-1 after 100 cycles) over the AQ powder electrodes. By augmenting the experimental findings via computational investigations, we are able to suggest design strategies that may foster the performance of industrially desirable powder-based electrode materials.

On the Accuracy of the Direct Method to Calculate pKa from Electronic Structure Calculations

The direct method (HA_(soln) ⇌ A_(soln)^– + H_(soln)⁺) for calculating pK_a of monoprotic acids is as efficient as thermodynamic cycles. A selective adjustment of proton free energy in solution was used with experimental pK_a data. The procedure was analyzed at different levels of theory. The solvent was described by the solvation model density (SMD) model, including or not explicit water molecules, and three training sets were tested. The best performance under any condition was obtained by the G4CEP method with a mean absolute error close to 0.5 units of pK_a and an uncertainty around ±1 unit of pK_a for any training set including or excluding explicit solvent molecules. PM6 and AM1 performed very well with average absolute errors below 0.75 units of pK_a but with uncertainties up to ±2 units of pK_a, using only the SMD solvent model. Density functional theory (DFT) results were highly dependent on the basis functions and explicit water molecules. The best performance was observed for the local spin density approximation (LSDA) functional in almost all calculations and under certain conditions, as high as those obtained by G4CEP. Basis set complexity and explicit solvent molecules were important factors to control DFT calculations. The training set molecules should consider the diversity of compounds.

Acetonitrile Transition Metal Interfaces from First Principles

Acetonitrile is among the most commonly used nonaqueous solvents in catalysis and electrochemistry. We study its interfaces with multiple facets of the metals Ag, Cu, Pt, and Rh using density functional theory calculations; the structures reported shed new light on experimental observations and underscore the importance of solvent-solvent interactions at high coverage. We investigate the relationship of potential of zero charge (PZC) to metal work function, reporting results in agreement with experimental measurements. We develop a model to explain the effects of solvent chemisorption and orientation on the PZC to within a mean absolute deviation of 0.08-0.12 V for all facets studied. Our electrostatic field dependent phase diagram agrees with spectroscopic observations and sheds new light on electrostatic field effects. This work provides new insight into experimental observations on acetonitrile metal interfaces and provides guidance for future studies of acetonitrile and other nonaqueous solvent interfaces with transition metals.

Absolute ion hydration free energy scale and the surface potential of water via quantum simulation

With a goal of determining an absolute free energy scale for ion hydration, quasi-chemical theory and ab initio quantum mechanical simulations are employed to obtain an accurate value for the bulk hydration free energy of the Na⁺ ion. The free energy is partitioned into three parts: 1) the inner-shell or chemical contribution that includes direct interactions of the ion with nearby waters, 2) the packing free energy that is the work to produce a cavity of size λ in water, and 3) the long-range contribution that involves all interactions outside the inner shell. The interfacial potential contribution to the free energy resides in the long-range term. By averaging cation and anion data for that contribution, cumulant terms of all odd orders in the electrostatic potential are removed. The computed total is then the bulk hydration free energy. Comparison with the experimentally derived real hydration free energy produces an effective surface potential of water in the range -0.4 to -0.5 V. The result is consistent with a variety of experiments concerning acid-base chemistry, ion distributions near hydrophobic interfaces, and electric fields near the surface of water droplets.

Cu2+ in liquid ammonia-The impact of solvent flexibility and electron correlation in ab initio quantum mechanical charge field molecular dynamics

The impact of solvent flexibility and electron correlation on the simulation results of Cu²⁺ in liquid ammonia has been investigated via an ab initio quantum mechanical charge field molecular dynamics (QMCF MD) simulation approach. To achieve this, three different simulation systems were considered in this study, namely Cu²⁺ in rigid and flexible ammonia at Hartree-Fock (HF) level of theory, as well as resolution of identity second order Møller-Plesset (MP2) perturbation theory in the rigid body case. In all cases, a stable octahedral [Cu(NH₃ )₆ ]²⁺ complex subject to dynamic Jahn-Teller distortions without the occurrence of ligand exchange was observed. The Cu²⁺  - NH₃ distance in the first shell agrees well with the experimental and other theoretical data. In all three cases, the structural data shows that the rigid-body ammonia model in conjunction with the HF level of theory provides accurate data for the first solvation shell, while at the same time, the computational demand and thus the achievable simulation time are much more beneficial. The vibrational analysis of the Cu²⁺  - NH₃ interaction yields similar force constants in the three investigated systems indicating that there is no distinct difference on the dynamical properties of the first solvation shell. In addition to the QMCF MD simulations, a number of natural bond orbital (NBO) analyses were carried out, confirming the strong electrostatic character of the Cu²⁺  - NH₃ interaction.

Computing pKa Shifts Using Traditional Molecular Dynamics: Example of Acid-Base Indicator Dyes in Organized Solutions

A compound's acidity constant (K_a) in a given medium determines its protonation state and, thus, its behavior and physicochemical properties. Therefore, it is among the key characteristics considered during the design of new compounds for the needs of advanced technology, medicine, and biological research, a notable example being pH sensors. The computational prediction of K_a for weak acids and bases in homogeneous solvents is presently rather well developed. However, it is not the case for more complex media, such as microheterogeneous solutions. The constant-pH molecular dynamics (MD) method is a notable contribution to the solution of the problem, but it is not commonly used. Here, we develop an approach for predicting K_a changes of weak small-molecule acids upon transfer from water to colloid solutions by means of traditional classical molecular dynamics. The approach is based on free energy (ΔG) computations and requires limited experiment data input during calibration. It was successfully tested on a series of pH-sensitive acid-base indicator dyes in micellar solutions of surfactants. The difficulty of finite-size effects affecting ΔG computation between states with different total charges is taken into account by evaluating relevant corrections; their impact on the results is discussed, and it is found non-negligible (0.1-0.4 pK_a units). A marked bias is found in the ΔG values of acid deprotonation, as computed from MD, which is apparently caused by force-field issues. It is hypothesized to affect the constant-pH MD and reaction ensemble MD methods as well. Consequently, for these methods, a preliminary calibration is suggested.

Embedded Mean-Field Theory for Solution-Phase Transition-Metal Polyolefin Catalysis

Decreasing the wall-clock time of quantum mechanics/molecular mechanics (QM/MM) calculations without sacrificing accuracy is a crucial prerequisite for widespread simulation of solution-phase dynamical processes. In this work, we demonstrate the use of embedded mean-field theory (EMFT) as the QM engine in QM/MM molecular dynamics (MD) simulations to examine polyolefin catalysts in solution. We show that employing EMFT in this mode preserves the accuracy of hybrid-functional DFT in the QM region, while providing up to 20-fold reductions in the cost per SCF cycle, thereby increasing the accessible simulation time-scales. We find that EMFT reproduces DFT-computed binding energies and optimized bond lengths to within chemical accuracy, as well as consistently ranking conformer stability. Furthermore, solution-phase EMFT/MM simulations provide insight into the interaction strength of strongly coordinating and bulky counterions.

Prediction of Solvation Free Energies of Ionic Solutes in Neutral Solvents

The prediction of solvation free energies is essential for a variety of applications. Solvation free energies of neutral systems can be predicted quite accurately. The accuracy of predictions for solvation free energies of ionic solutes dissolved in neutral solvents, however, has been reported to be worse by at least 1 order of magnitude. In this study, the performance of three approaches for solvation free energy prediction of several hundred ions dissolved in neutral solvents is evaluated. The applied methods are COSMO-RS, cluster continuum model (CCM) together with COSMO-RS, and COSMO-RS-ES. It is emphasized that the reference data for model evaluation are subject to large uncertainties stemming from the impossibility to measure the so-called elusive absolute free energies of solvation of a single ion. Consequently, such uncertainty must be considered during the evaluation of prediction methods. Therefore, a straightforward approach to account for the underlying uncertainty is applied here. Hereby, it is revealed that the true performance of the method is better than what is often reported. The average absolute deviation (AAD) of COSMO-RS is calculated to be 2.3 kcal mol^-1, while applying the CCM and COSMO-RS-ES each results in AADs of 2.0 kcal mol^-1. This accuracy allows for qualitative assessment of solvation free energy-dependent quantities, such as reaction rate constants.

Hybrid Solvation Model with First Solvation Shell for Calculation of Solvation Free Energy

We present a hybrid solvation model with first solvation shell to calculate solvation free energies. This hybrid model combines the quantum mechanics and molecular mechanics methods with the analytical expression based on the Born solvation model to calculate solvation free energies. Based on calculated free energies of solvation and reaction profiles in gas phase, we set up a unified scheme to predict reaction profiles in solution. The predicted solvation free energies and reaction barriers are compared with experimental results for twenty bimolecular nucleophilic substitution reactions. These comparisons show that our hybrid solvation model can predict reliable solvation free energies and reaction barriers for chemical reactions of small molecules in aqueous solution.

First-principles modeling of chemistry in mixed solvents: Where to go from here?

Mixed solvents (i.e., binary or higher order mixtures of ionic or nonionic liquids) play crucial roles in chemical syntheses, separations, and electrochemical devices because they can be tuned for specific reactions and applications. Apart from fully explicit solvation treatments that can be difficult to parameterize or computationally expensive, there is currently no well-established first-principles regimen for reliably modeling atomic-scale chemistry in mixed solvent environments. We offer our perspective on how this process could be achieved in the near future as mixed solvent systems become more explored using theoretical and computational chemistry. We first outline what makes mixed solvent systems far more complex compared to single-component solvents. An overview of current and promising techniques for modeling mixed solvent environments is provided. We focus on so-called hybrid solvation treatments such as the conductor-like screening model for real solvents and the reference interaction site model, which are far less computationally demanding than explicit simulations. We also propose that cluster-continuum approaches rooted in physically rigorous quasi-chemical theory provide a robust, yet practical, route for studying chemical processes in mixed solvents.

Mean Inner Potential of Liquid Water

Improving our experimental and theoretical knowledge of electric potentials at liquid-solid boundaries is essential to achieve a deeper understanding of the driving forces behind interfacial processes. Electron holography has proved successful in probing solid-solid interfaces but requires knowledge of the materials' mean inner potential (MIP, V_{0}), which is a fundamental bulk material property. Combining off-axis electron holography with liquid phase transmission electron microscopy (LPTEM), we provide the first quantitative MIP determination of liquid water V_{0}=+4.48±0.19  V. This value is larger than most theoretical predictions, and to explain the disagreement we assess the dominant factors needed in quantum simulations of liquid water. A precise MIP lays the foundations for nanoscale holographic potential measurements in liquids, and provides a benchmark to improve quantum mechanical descriptions of aqueous systems and their interfaces in, e.g., electrochemistry, solvation processes, and spectroscopy.

Electromechanics of the liquid water vapour interface

Two collective properties distinguishing the thin liquid water vapour interface from the bulk liquid are the anisotropy of the pressure tensor giving rise to surface tension and the orientational alignment of the molecules leading to a finite dipolar surface potential. Both properties can be regarded as capillary phenomena and are likely to be coupled. We have investigated this coupling by determining the response of the tangential component of the surface tension to the application of an electric field normal to the surface using finite field molecular dynamics simulations. We find an upside down parabola with a maximum shifted away from zero field. Comparing the molecular dynamics results to a phenomenological electromechanical model we relate the zero field derivative of the tangential part of the surface tension to the electrostatic potential generated by the spontaneous interface polarization. When interpreted with this model our simulations also indicate that Kelvin forces due to electric field gradients at a polarized interface play an important role in the effective dielectric response.

Large-Sized Ammonia Clusters and Solvation Energies of the Proton in Ammonia

The absolute solvation energies (free energies and enthalpies) of the proton in ammonia are used to compute the pK_a of species embedded in ammonia. They are also used to compute the solvation energies of other ions in ammonia. Despite their importance, it is not possible to determine experimentally the solvation energies of the proton in a given solvent. We propose in this work a direct approach to compute the solvation energies of the proton in ammonia from large-sized neutral and protonated ammonia clusters. To undertake this investigation, we performed a geometry optimization of neutral and protonated ammonia 30-mer, 40-mer, and 50 mer to locate stable structures. These structures have been fully optimized at both APFD/6-31++g(d,p) and M06-2X/6-31++g(d,p) levels of theory. An infrared spectroscopic study of these structures has been provided to assess the reliability of our investigation. Using these structures, we have computed the absolute solvation free energy and the absolute solvation enthalpy of the proton in ammonia. It comes out that the absolute solvation free energy of the proton in ammonia is calculated to be -1192 kJ mol^-1 , whereas the absolute solvation enthalpy is evaluated to be -1214 kJ mol^-1 . © 2019 Wiley Periodicals, Inc.

Quantifying the hydration structure of sodium and potassium ions: taking additional steps on Jacob's Ladder

The ability to reproduce the experimental structure of water around the sodium and potassium ions is a key test of the quality of interaction potentials due to the central importance of these ions in a wide range of important phenomena.

Solvation energies of ions with ensemble cluster-continuum approach

An alternative cluster-continuum approach for the calculation of solvation free energies of ions.

A Calculation Model of the General Theory of Interaction Potentials for Stoichiometric Lanthanide Type Crystals: Applications to the Cs2KLnCl6 System

This article aims to develop a generalized model calculation model to be applicable to the general theory of interaction potentials with reference to the stoichiometric elpasolite type crystals. In this study, we have chosen to report both a theoretical model and a calculation strategy to undertake semi empirical calculations of thermodynamic properties, such as reticular energies and heats of formation for the series of systems such as: Cs₂KLnCl₆. We have also carried out quite a number of calculations for a variety of systems such as: Cs₂NaLnF₆, Cs₂NaLnCl₆, Cs₂NaLnBr₆, Rb₂NaLnF₆and Cs₂KLnF₆ in the Fm3m space group since we aim to check the strengths and weaknesses of our model calculations. We have analyzed a substantial number of approximate theoretical models and have carried a formidable amount of computing simulations to estimate the reticular energies and the corresponding heat of formation for these type of crystal using a semi empirical model. We made use of the thermodynamic cycles of Born-Haber so as to get a broad view with reference to the accuracy of our semi empirical theoretical models. The problem itself is quite challenging since we have focused our attention upon trivalent lanthanide ions \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L{n}^{+3}$$\end{document}Ln+3 in the first inner transition series of the chemical elements: (Ce, Pr, Nd, Pm, Sm, Eu, Gd, Tb, Dy, Ho, Er, Tm, Yb, Lu). There are a significant amount of outstanding research works published in the literature with reference to structural analysis, one photon spectroscopy, vibrionic intensity model calculations and generalized models to deal with these kind of complex crystals. The calculated energy values associated with these observables seems to be most reasonable, and these follow the expected trends, as may be expected on both theoretical and experimental grounds. Both, the advantages and disadvantages of the current model calculations, have been tested against other previous calculations performed for this type of complex systems. It is of a paramount importance, the results obtained and reported in this article with regards to convergence tests as well as some master equations derived to account for the various contributions to the total energy. The Born-Mayer-Buckingham potential is carefully examined with reference to these lanthanide type crystals Cs₂KLnCl₆. Finally but not at last, the most likely sources for improvement are carefully discussed in this work. We strongly believe that there is enough room for improvement and have therefore initiated a new research program of activities tackling systems of well-known optical and structural properties.

An effective partial charge model for bulk and surface properties of cubic ZrO2, Y2O3 and yttrium-stabilised zirconia

In this work a newly parametrised Coulomb plus Buckingham potential formulation for cubic ZrO₂, Y₂O₃ and yttrium-stabilised zirconia (YSZ) is presented. The density and pair distributions obtained for neat ZrO₂ and Y₂O₃ under ambient conditions are in excellent agreement with experimental data, while the vibrational power spectra are highly similar compared to those obtained via ab initio molecular dynamics simulations at the PBEsol level. In addition, it is shown that the use of effective partial charges has several advantages compared to interaction potentials employing the oxidation states in the evaluation of the coulombic interactions: (i) the diffusion coefficient and the associated activation energy of oxygen ions evaluated for YSZn (n = 4 to 12) display the best agreement with experimental data; (ii) no unphysical reorganisation of the interface and the bulk are observed in simulations of the (110) and (111) surfaces of cubic ZrO₂ and Y₂O₃, while due to the strong coulombic contributions in the case of the tested full-charge models a pronounced restructuring of the interface and the bulk is observed in the ZrO₂ case, and (iii) the use of effective partial charges ensures compatibility with existing solvent models and force-fields for the treatment of molecular compounds.

The Jahn-Teller effect in mixed aqueous solution: the solvation of Cu2+ in 18.6% aqueous ammonia obtained from ab initio quantum mechanical charge field molecular dynamics

The solvation structure and dynamics of Cu2+ in 18.6 % aqueous ammonia have been investigated using an ab initio quantum mechanical charge field molecular dynamics (QMCF MD) simulation approach at the Hartree–Fock (HF) level of theory applying the LANL2DZ ECP and Dunning DZP basis sets for Cu2+, ammonia and water, respectively. During a simulation time of 20 ps, only NH3 molecules are observed within the first solvation shell of Cu2+, resulting in the formation of an octahedral [Cu(NH3)6]2+ complex. While no exchange of these ligands with the second solvation shell are observed along the simulation, the monitoring of the associated N-Ntrans distances highlight the dynamics of the associated Jahn-Teller distortions, showing on average 2 elongated axial (2.19 Å) and 4 equatorial Cu–N bonds (2.39 Å). The observed structural properties are found in excellent agreement with experimental studies. In addition, an NBO analysis was carried out, confirming the strong electrostatic character of the Cu2+–NH3 interaction.