Direct Absolute pKa Predictions and Proton Transfer Mechanisms of Small Molecules in Aqueous Solution by QM/MM-MD

Variations in proton transfer pathways and energetics on pristine and defect-rich quartz surfaces in water: Insights into the bimodal acidities of quartz

HYPOTHESIS

Understanding the mechanisms of proton transfer on quartz surfaces in water is critical for a range of processes in geochemical, environmental, and materials sciences. The wide range of surface acidities (>9 pKa units) found on the ubiquitous mineral quartz is caused by the structural variations of surface silanol groups. Molecular scale simulations provide essential tools for elucidating the origin of site-specific surface acidities.

SIMULATIONS

We used density-functional tight-binding-based molecular dynamics combined with rare-event metadynamics simulations to probe the mechanisms of deprotonation reactions from ten representative surface silanol groups found on both pristine and defect-rich quartz (101) surfaces with Si vacancies.

FINDINGS

The results show that deprotonation is a highly dynamic process where both the surface hydroxyls and bridging oxygen atoms serve as the proton acceptors, in addition to water. Deprotonation of embedded silanols through intrasurface proton transfer exhibited lower pKa values with less H-bond participation and higher energy barriers, suggesting a new mechanism to explain the bimodal acidity observed on quartz surface. Defect sites, recently shown to comprise a significant portion of the quartz (101) surface, diversify the coordination and local H-bonding environments of the surface silanols, changing both the deprotonation pathways and energetics, leading to a wider range of pKa values (2.4 to 11.5) than that observed on pristine quartz surface (10.4 and 12.1).

Identifying the acidic or basic behavior of surface water: a QM/MM-MD study

Controversies on the water surface were theoretically addressed with the help of large scale quantum mechanical molecular dynamics (QMMD) simulations on water surface model systems with and without excess hydroniums and hydroxides. It was revealed that the thermodynamic surface structures of these ions strongly depend on their location and dipole orientation. Fast hydronium diffusion by proton transfer establishes a wider kinetic depth distribution (∼6 Å) than that predicted by its thermodynamic affinity for the water surface, while slow hydroxide is shallowly trapped below the outermost molecular layer (3-4 Å). In addition, the anisotropic orientation of surface water dipole can generate a substantial magnitude of surface potential, which extends to a depth of a few molecular layers. With these distinctively different surface properties of two ions and water molecules, the seemingly contradictory observations of acidic and negatively charged water surfaces may be successfully explained. That is, the negative surface charge of neutral water mostly stems from intrinsic water properties such as water dipole orientation and electron density spillage at the surface, rather than surface OH- ions. The enhanced acidity of the water surface can be attributed in large part to the kinetic depth profile of ion density in addition to static thermodynamic origin. Furthermore, the different depth profiles of the two ions may differently affect the surface-sensitive spectroscopic observations.

New insights into the degradation mechanism of ibuprofen in the UV/H2O2 process: role of natural dissolved matter in hydrogen transfer reactions

Ibuprofen (IBU), a widely used antipyretic and analgesic, has been frequently detected in various natural water systems. Advanced oxidation processes (AOPs) are effective ways to remove pollutants from water. The degradation of IBU under UV/H2O2 conditions in the presence of various kinds of natural dissolved matter was investigated using density functional theory (DFT). The eco-toxicological properties were predicted based on a quantitative structure-activity relationship (QSAR) model. The calculated results showed that two H-abstraction reactions occurring at the side chain are predominant pathways in the initial reaction. H2O, NH3, CH3OH, C2H5OH, HCOOH and CH3COOH can catalyze the H transfer in the degradation process through decreasing the energy barriers and the catalysis effects follow the order of NH3 > alcohols > acids > H2O. The catalysis effects differ under acid or alkaline conditions. The overall rate coefficient of the reaction of IBU with ˙OH is calculated to be 5.04 × 109 M-1 s-1 at 298 K. IBU has harmful effects on aquatic organisms and human beings and the degradation process cannot significantly reduce its toxicity. Among all products, 2-(4-formylphenyl)propanoic acid, which is more toxic than IBU, is the most toxic with acute and chronic toxicity, developmental toxicity, mutagenicity, genotoxic carcinogenicity and irritation/corrosivity to skin. The findings in this work provide new insights into the degradation of IBU and can help to assess its environmental risks.

Accurate Computation of Aqueous pKas of Biologically Relevant Organic Acids: Overcoming the Challenges Posed by Multiple Conformers, Tautomeric Equilibria, and Disparate Functional Groups with the Fully Black-Box pK-Yay Method

The use of static electronic structure calculations to compute solution-phase pKas offers a great advantage in that a macroscopic bulk property could be computed via microscopic computations involving very few molecules. There are various sources of errors in the quantum chemical calculations though. Overcoming these errors to accurately compute pKas of a plethora of acids is an active area of research in physical chemistry pursued by both computational as well as experimental chemists. We recently developed the pK-Yay method in our attempt to accurately compute aqueous pKas of strong and weak acids. The method is fully black-box, computationally inexpensive, and is very easy for even a nonexpert to use. However, the method was thus far tested on very few molecules (only 16 in all). Herein, in order to assess the future applicability of pK-Yay, we study the effect of multiple conformers, the presence of tautomers under equilibrium, and the impact of a wide variety of functional groups (derivatives of acetic acid with substituents at various positions, dicarboxylic acids, aromatic carboxylic acids, amines and amides, phenols and thiols, and fluorine bearing organic acids). Starting with more than 1000 conformers and tautomers, this study establishes that overall errors of ∼ 1.0 pKa units are routinely obtained for a majority of the molecules. Larger errors are noted in cases where multiple charges, intramolecular hydrogen bonding, and several ionizable functional groups are simultaneously present. An important conclusion to emerge from this work is that, the computed pKas are insensitive (difference <0.5) to whether we consider multiple conformers/tautomers or only choose the most stable conformer/tautomer. Further, pK-Yay captures the stereoelectronic effects arising due to differing axial vs equatorial pattern, and is useful to predict the dominant acid-base equilibrium in a system featuring several equilibria. Overall, pK-Yay may be employed in several chemical applications featuring organic molecules and biomonomers.

Speciation of La3+-Cl- Complexes in Hydrothermal Fluids from Deep Potential Molecular Dynamics

The stability of rare earth element (REE) complexes plays a crucial role in quantitatively assessing their hydrothermal migration and transformation. However, reliable data are lacking under high-temperature hydrothermal conditions, which hampers our understanding of the association behavior of REE. Here a deep learning potential model for the LaCl3-H2O system in hydrothermal fluids is developed based on the first-principles density functional theory calculations. The model accurately predicts the radial distribution functions compared to ab initio molecular dynamics (AIMD) simulations. Furthermore, species of La-Cl complexes, the dissociation pathway of the La-Cl complexes dissociation process, and the potential of mean forces and corresponding association constants (logK) for LaCln3-n (n = 1-4) are extensively investigated under a wide range of temperatures and pressures. Empirical density models for logK calculation are fitted with these data and can accurately predict logK data from both experimental results and AIMD simulations. The distribution of La-Cl species is also evaluated across a wide range of temperatures, pressures, and initial chloride concentration conditions. The results show that La-Cl complexes are prone to forming in a low-density solution, and the number of bonded Cl- ions increases with rising temperature. In contrast, in a high-density solution, La3+ dominates and becomes the more prevalent species.

Predicting Solvent Effects on SN2 Reaction Rates: Comparison of QM/MM, Implicit, and MM Explicit Solvent Models

Solvents are one of the key variables in the optimization of a synthesis yield or properties of a synthesis product. In this paper, contemporary solvent models are applied to predict the rates of S_N2 reactions in a range of aqueous and non-aqueous solvents. High-level CCSD(T)/CBS//M06-2X/6-31+G(d) gas phase energies were combined with solvation free energies from SMD, SM12, and ADF-COSMO-RS continuum solvent models, as well as molecular mechanics (MM) explicit solvent models with different atomic charge schemes to predict the rate constants of three S_N2 reactions in eight protic and aprotic solvents. It is revealed that the prediction of rate constants in organic solvents is not necessarily less challenging than in water and popular solvent models struggle to predict their rate constants to within 3 log units of experimental values. Among the continuum solvent models, the ADF-COSMO-RS model performed the best in predicting absolute rate contants while the SM12 model was best at predicting relative rate constants with an average accuracy of about 1.5 and 0.8 log units, respectively. The use of computationally more demanding MM explicit solvent models did not translate to improvements in absolute rate constants but was quite effective at predicting relative rate constants due to systematic error cancellation. Free energy barriers obtained from umbrella sampling with explicit solvent QM/MM simulations led to excellent agreement with experimental values, provided that a validated level of theory is used to treat the QM region.

Accelerated Deep Learning Dynamics for Atomic Layer Deposition of Al(Me)3 and Water on OH/Si(111)

Knowledge of the detailed mechanism behind the atomic layer deposition (ALD) can greatly facilitate the optimization of the manufacturing process. Computational modeling can potentially foster the understanding; however, the presently available capabilities of the accurate ab initio computational techniques preclude their application to modeling surface processes occurring on a long time scale, such as ALD. Although the situation can be greatly improved using machine learning (ML), this technique requires an enormous amount of data for training datasets. Here, we propose an iterative protocol for optimizing ML training datasets and apply ML-assisted ab initio calculations to model surface reactions occurring during the Al(Me)₃/H₂O ALD process on the OH-terminated Si (111) surface. The protocol uses a recently developed low-dimensional projection technique (TDUS), greatly reducing the amount of information required to achieve high accuracy (ca. 1 kcal/mol or less) of the developed ML models. The resulting free energy landscapes reveal fine details of various aspects of the target ALD process, such as the surface proton transfer, zwitterionic surface configurations, elimination-addition/addition-elimination, and S_N2 reactions as well as the role of the surface entropic and temperature effects. Simulations of adsorption dynamics predict that the maximum physisorption rate of ca. 70% is achieved at the incidence velocity u_rms of the reactants in the range of 15-20 Å/ps. Hence, the proposed protocol furnishes a very effective tool to study complex chemical reaction dynamics at a much reduced computational cost.

Advanced Theory and Simulation to Guide the Development of CO2 Capture Solvents

Increasing atmospheric concentrations of greenhouse gases due to industrial activity have led to concerning levels of global warming. Reducing carbon dioxide (CO₂) emissions, one of the main contributors to the greenhouse effect, is key to mitigating further warming and its negative effects on the planet. CO₂ capture solvent systems are currently the only available technology deployable at scales commensurate with industrial processes. Nonetheless, designing these solvents for a given application is a daunting task requiring the optimization of both thermodynamic and transport properties. Here, we discuss the use of atomic scale modeling for computing reaction energetics and transport properties of these chemically complex solvents. Theoretical studies have shown that in many cases, one is dealing with a rich ensemble of chemical species in a coupled equilibrium that is often difficult to characterize and quantify by experiment alone. As a result, solvent design is a balancing act between multiple parameters which have optimal zones of effectiveness depending on the operating conditions of the application. Simulation of reaction mechanisms has shown that CO₂ binding and proton transfer reactions create chemical equilibrium between multiple species and that the agglomeration of resulting ions and zwitterions can have profound effects on bulk solvent properties such as viscosity. This is balanced against the solvent systems needing to perform different functions (e.g., CO₂ uptake and release) depending on the thermodynamic conditions (e.g., temperature and pressure swings). The latter constraint imposes a "Goldilocks" range of effective parameters, such as binding enthalpy and pK _a, which need to be tuned at the molecular level. The resulting picture is that solvent development requires an integrated approach where theory and simulation can provide the necessary ingredients to balance competing factors.

Connecting Gas-Phase Computational Chemistry to Condensed Phase Kinetic Modeling: The State-of-the-Art

In recent decades, quantum chemical calculations (QCC) have increased in accuracy, not only providing the ranking of chemical reactivities and energy barriers (e.g., for optimal selectivities) but also delivering more reliable equilibrium and (intrinsic/chemical) rate coefficients. This increased reliability of kinetic parameters is relevant to support the predictive character of kinetic modeling studies that are addressing actual concentration changes during chemical processes, taking into account competitive reactions and mixing heterogeneities. In the present contribution, guidelines are formulated on how to bridge the fields of computational chemistry and chemical kinetics. It is explained how condensed phase systems can be described based on conventional gas phase computational chemistry calculations. Case studies are included on polymerization kinetics, considering free and controlled radical polymerization, ionic polymerization, and polymer degradation. It is also illustrated how QCC can be directly linked to material properties.

pK-Yay: A Black-Box Method Using Density Functional Theory and Implicit Solvation Models to Compute Aqueous pKa Values of Weak and Strong Acids

The concept of pK_a is very important in chemistry and biology. Over the past two decades, electronic structure calculations have made huge strides toward becoming reliable counterparts to experiments in determining solution-phase pK_a values. However, invariably, the computation of pK_a values involves the use of methods whose error bars are intrinsically larger (definitely >2 kcal/mol with density functionals) than the accuracy desired for the estimation of pK_a values (<1 or 1.5 pK_a units). This scenario presents an ample scope for innovation in developing systematic error cancellation methods even today. In this work, we have developed the pK-Yay method. It is a user-friendly black-box method used to compute aqueous pK_a values of strong and weak acids employing routinely used and computationally inexpensive density functionals and implicit solvation models. It does not require the use of any explicit solvent molecule or modifying any other parameter in an electronic structure program. As part of evaluating the method, a comprehensive test set of 26 weak and strong organic and mineral acids covering 35 pK_a units (20 to -15) was assembled. The detailed description of the method, its performance with different functionals (ωB97X-D performed best for organic acids with a mean absolute error (MAE) of 0.8 pK_a units, and B2-PLYP performed best for strong acids with an MAE of 1.6 pK_a units), the strengths and limitations in the present version, and a future scope to improve the accuracy and reduce its empirical nature are presented herein.

The Impact of Electron Correlation on Describing QM/MM Interactions in the Attendant Molecular Dynamics Simulations of CO in Myoglobin

The impact of the dispersion and electron correlation effects on describing quantum mechanics/molecular mechanics (QM/MM) interactions in QM/MM molecular dynamics (MD) simulations was explored by performing a series of up to 2 ns QM/MM MD simulations on the B states of the myoglobin-carbon monoxide (MbCO) system. The results indicate that both dispersion and electron correlations play significant roles in the simulation of the ratios of two B states (B₁/B₂), which suggests that the inclusion of the electron correlation effects is essential for accurately modeling the interactions between QM and MM subsystems. We found that the QM/MM interaction energies between the CO and the surroundings statistically present a linear correlation with the electric fields along the CO bond. This indicates that QM/MM interactions can be described by a simple physical model of a dipole with constant moment under the action of the electric fields. The treatment provides us with an accurate and effective approach to account for the electron correlation effects in QM/MM MD simulations.

Low-dimensional projection approach for efficient sampling of molecular recognition and polymer aggregation

The one-dimensional projection (ODP) approach is extended to two-dimensional umbrella sampling (TDUS) and is applied to three different complex systems in combination with a reactive force field (ReaxFF). TDUS is capable of showing detailed features of the free-energy surface (FES) of the double-proton transfer of the acetic acid dimer. It also revealed the direct relationship between the types of hydrogen bonding and binding strengths in the case of adrenaline molecular recognition by SIVSF (Serine, Isoleucine, Valine, Cysteine, and Phenylalanine). The study of polymer aggregation using TDUS shows that aggregation is preferred with a less-polar solvent, which is also consistent with the experimental observation of a tape-casting process. Therefore, TDUS can be generally useful in FES explorations from simple chemical reactions to complex processes of molecular recognition and polymer aggregation.

How accurate are approximate quantum chemical methods at modelling solute-solvent interactions in solvated clusters?

In this paper, the performance of a wide range of DFT methods is assessed for the calculation of interaction energies of thermal clusters of a solute in water. Three different charge states (neutral, proton transfer transition state and zwitterion) of glycine were solvated by 1 to 40 water molecules as sampled from molecular dynamics simulations. While some ab initio composite methods that employ insufficiently large basis sets incurred significant errors even for a cluster containing only 5 water molecules relative to the W1X-2 benchmark, the DLPNO-CCSD(T)/CBS and DSD-PBEP86 (triple zeta basis set) levels of theory predicted very accurate interaction energies. These levels of theory were used to benchmark the performance of 16 density functionals from different rungs of Jacob's Ladder. Of the Rung 4 functionals examined, the ωB97M-V and ωB97X-V functionals stood out for predicting absolute interaction energies in 40-water clusters with mean absolute deviations (MAD) ∼4 kJ mol-1. The B3LYP-D3(BJ) functional performed exceptionally well with a MAD ∼1.7 kJ mol-1 and is the overall best performing method. Calculations of relative interaction energies allow for cancellation of systematic errors, including basis set truncation and superposition errors, and the ωB97M-V and B3LYP-D3(BJ) double zeta basis set calculations yielded relative interaction energies that are within ∼3 kJ mol-1 of the benchmark. The ONIOM approximation provides another strategy for accelerating the calculation of accurate absolute interaction energies provided that the calculations have converged with respect to the size of the "high-level-layer".

QM/MM Investigation of the Role of a Second Coordination Shell Arginine in [NiFe]-Hydrogenases

[NiFe]-hydrogenases are highly efficient catalysts for the heterolytic splitting of molecular hydrogen (H₂). The heterobimetallic cysteine-coordinated active site of these enzymes is covered by a highly conserved arginine residue, whose role in the reaction is not fully resolved yet. The structural and catalytic role of this arginine is investigated here using QM/MM calculations with various exchange-correlation functionals. All of them give a very consistent picture of the thermodynamics of H₂ oxidation. The concept of the presence of a neutral arginine and its direct involvement as a Frustrated Lewis Pair (FLP) in the reaction is critically evaluated. The arginine, however, would exist in its standard protonation state and perform a critical role in positioning and slightly polarizing the substrate H₂. It is not directly involved in the heterolytic processing of H₂ but guides its approach and reduces its flexibility during binding. Upon substitution of the positively charged arginine by a charge-conserving lysine residue, the H₂ binding position remains unaffected. However, critical hydrogen bonding interactions with nearby aspartate residues are lost. In addition, the H₂ polarization is unfavorable and the reduced side-chain volume may negatively affect the kinetics of the catalytic process.

An efficient protocol for computing the pKa of Zn-bound water

We present an efficient and accurate method for computing absolute pK_w values in Zn²⁺ complexes.

Steric Effects of pH Switchable, Substituted (2-pyridinium)urea Organocatalysts: a Solution and Solid Phase Study

The study of hydrogen bonding organocatalysis is rapidly expanding. Much research has been directed at making catalysts more active and selective, with less attention on fundamental design strategies. This study systematically increases steric hindrance at the active site of pH switchable urea organocatalysts. Incorporating strong intramolecular hydrogen bonds from protonated pyridines to oxygen stabilizes the active conformation of these ureas thus reducing the entropic penalty that results from substrate binding. The effect of increasing steric hindrance was studied by single crystal X-ray diffraction and by kinetics experiments of a benchmark reaction.

In search of the best DFT functional for dealing with organic anionic species

“And the winner is…” This work assesses the ability of different Density Functional Theory (DFT) functionals for a proper treatment of organic anionic species.

Polyelectrolyte pKa from experiment and molecular dynamics simulation

The pK_a of a polyelectrolyte has been determined experimentally by potentiometric titration and computed using Molecular Dynamics (MD) constant pH (CpH) methodology, which allows the pK_a of each titratable site along the polymer backbone.

Origin of Acid-Base Catalytic Effects on Formaldehyde Hydration

The mechanisms of hydronium- and hydroxide-catalyzed formaldehyde hydrations were investigated by quantum mechanical/molecular mechanical molecular dynamics in combination with flexible coordinates. A stepwise bimolecular and a concerted termolecular mechanism were found with a hydronium catalyst. The latter is more favorable and better consistent with experiment. Structurally, a dipole-bound species initially arranges the nucleophile in a favorable configuration for both routes, significantly enhancing the reactive collisions. On the one hand, the hydronium catalyst also plays a role of a reactant in the bimolecular path. On the other hand, only a stepwise mechanism was found with a hydroxide catalyst. Overall, hydroxide is a stronger catalyst than a hydronium when it is in contact distance with formaldehyde.

Role of Deprotonation Free Energies in pKa Prediction and Molecule Ranking

A computationally efficient classical molecular simulation technique is derived for ranking the pKa values of a set of chemically similar congeneric molecules in an implicit solvent model of water. This uses the deprotonation free energy of the titratable group in the gas and aqueous phases obtained by thermodynamic integration (TI). For a series of alcohols and acids a strong linear correlation is demonstrated between the experimental pKa and the deprotonation free energy difference in the gas and liquid phases. These calculations also show that classical TI is more efficient than slow-growth TI in calculating deprotonation free energies for the series of molecules considered herein.

Advantages of Organic Halogen Bonding for Halide Recognition

The study of hydrogen bonding organocatalysis is rapidly expanding. Much research has been directed at making catalysts more active and selective, with less attention on fundamental design strategies. This study systematically increases steric hindrance at the active site of pH switchable urea organocatalysts. Incorporating strong intramolecular hydrogen bonds from protonated pyridines to oxygen stabilizes the active conformation of these ureas thus reducing the entropic penalty that results from substrate binding. The effect of increasing steric hindrance was studied by single crystal X-ray diffraction and by kinetics experiments of a benchmark reaction.

Ab Initio Molecular Dynamics Simulations of Amino Acids in Aqueous Solutions: Estimating pKa Values from Metadynamics Sampling

Changes in the protonation and deprotonation of amino acid residues in proteins play a key role in many biological processes and pathways. Here, we report calculations of the free-energy profile for the protonation-deprotonation reaction of the 20 canonical α amino acids in aqueous solutions using ab initio Car-Parrinello molecular dynamics simulations coupled with metadynamics sampling. We show here that the calculated change in free energy of the dissociation reaction provides estimates of the multiple pKa values of the amino acids that are in good agreement with experiment. We use the bond-length-dependent number of the protons coordinated to the hydroxyl oxygen of the carboxylic and the amine groups as the collective variables to explore the free-energy profiles of the Bronsted acid-base chemistry of amino acids in aqueous solutions. We ensure that the amino acid undergoing dissociation is solvated by at least three hydrations shells with all water molecules included in the simulations. The method works equally well for amino acids with neutral, acidic and basic side chains and provides estimates of the multiple pKa values with a mean relative error, with respect to experimental results, of 0.2 pKa units.

Estimating successive pKa values of polyprotic acids from ab initio molecular dynamics using metadynamics: the dissociation of phthalic acid and its isomers

pKa values of polyprotic acids estimated from ab initio molecular dynamics computation of the change in free energy upon dissociation.

Like-charge ion pairs of hydronium and hydroxide in aqueous solution?

Hydroxide can form stronger like-ion pairs than hydronium in aqueous solution mostly due to its versatile coordination ability with solvents.

Calculations of pKa's and redox potentials of nucleobases with explicit waters and polarizable continuum solvation

The SMD implicit solvation model augmented with one and four explicit water molecules was used to calculate pKa's and redox potentials of N-methyl-substituted nucleic acid bases guanine, adenine, cytosine, thymine, and uracil. Calculations were carried out with the B3LYP/6-31+G(d,p) level of theory. The same numbers of water molecules were hydrogen bonded to the neutral, protonated, and deprotonated nucleobases in their unoxidized and oxidized forms. The improvement in pKa1 involving neutrals and cations was modest. By contrast, the improvement in pKa2 involving neutrals and anions was quite significant, reducing the mean absolute error from 4.6 pKa units with no waters, to 2.6 with one water and 1.7 with four waters. For the oxidation of nucleobases, adding explicit waters did little to improve E(X(•),H(+)/XH), possibly because both species in the redox couple are neutral molecules at pH 7.

Multiscale Reactive Molecular Dynamics for Absolute pKa Predictions and Amino Acid Deprotonation

Accurately calculating a weak acid's pK_a from simulations remains a challenging task. We report a multiscale theoretical approach to calculate the free energy profile for acid ionization, resulting in accurate absolute pK_a values in addition to insights into the underlying mechanism. Importantly, our approach minimizes empiricism by mapping electronic structure data (QM/MM forces) into a reactive molecular dynamics model capable of extensive sampling. Consequently, the bulk property of interest (the absolute pK_a) is the natural consequence of the model, not a parameter used to fit it. This approach is applied to create reactive models of aspartic and glutamic acids. We show that these models predict the correct pK_a values and provide ample statistics to probe the molecular mechanism of dissociation. This analysis shows changes in the solvation structure and Zundel-dominated transitions between the protonated acid, contact ion pair, and bulk solvated excess proton.

Predicting pKa in Implicit Solvents: Current Status and Future Directions

Computational prediction of condensed phase acidity is a topic of much interest in the field today. We introduce the methods available for predicting gas phase acidity and pKas in aqueous and non-aqueous solvents including high-level electronic structure methods, empirical linear free energy relationships (LFERs), implicit solvent methods, explicit solvent statistical free energy methods, and hybrid implicit–explicit approaches. The focus of this paper is on implicit solvent methods, and we review recent developments including new electronic structure methods, cluster-continuum schemes for calculating ionic solvation free energies, as well as address issues relating to the choice of proton solvation free energy to use with implicit solvation models, and whether thermodynamic cycles are necessary for the computation of pKas. A comparison of the scope and accuracy of implicit solvent methods with ab initio molecular dynamics free energy methods is also presented. The present status of the theory and future directions are outlined.