Advantages of Multidimensional Biasing in Accelerated Dynamics: Application to the Calculation of the Acid pKa for Acetic Acid

The use of accelerated sampling methods such as metadynamics has shown a significant advantage in calculations that involve infrequent events, which would otherwise require sampling a prohibitive number of configurations to determine the difference in free energies between two or more chemically distinct states such as in the calculation of acid dissociation constants Ka. In this case, the most common method is to bias the system via a single collective variable (CV) representing the coordination number of the proton donor group, which yields results in reasonable agreement with experiments. Here we study the deprotonation of acetic acid using the reactive force field ReaxFF and observe a significant dependence of Ka on the simulation box size when biasing only the coordination number CV, which is due to incomplete sampling of the deprotonated state for small simulation systems and inefficient sampling for larger ones. Incorporating a second CV representing the distance between the H3O+ cation and the acetate anion results in substantially more efficient sampling, both accelerating the dynamics and virtually eliminating the computational box size dependence.

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.

Predicting pKa Values of Quinols and Related Aromatic Compounds with Multiple OH Groups

Quinonoid compounds play central roles as redox-active agents in photosynthesis and respiration and are also promising replacements for inorganic materials currently used in batteries. To design new quinonoid compounds and predict their state of protonation and redox behavior under various conditions, their pK_a values must be known. Methods that can predict the pK_a values of simple phenols cannot reliably handle complex analogues in which multiple OH groups are present and may form intramolecular hydrogen bonds. We have therefore developed a straightforward method based on a linear relationship between experimental pK_a values and calculated differences in energy between quinols and their deprotonated forms. Simple adjustments allow reliable predictions of pK_a values when intramolecular hydrogen bonds are present. Our approach has been validated by showing that predicted and experimental values for over 100 quinols and related compounds differ by an average of only 0.3 units. This accuracy makes it possible to select proper pK_a values when experimental data vary, predict the acidity of quinols and related compounds before they are made, and determine the sites and orders of deprotonation in complex structures with multiple OH groups.

Chemistrees: Data-Driven Identification of Reaction Pathways via Machine Learning

We propose to analyze molecular dynamics (MD) output via a supervised machine learning (ML) algorithm, the decision tree. The approach aims to identify the predominant geometric features which correlate with trajectories that transition between two arbitrarily defined states. The data-driven algorithm aims to identify these features without the bias of human “chemical intuition”. We demonstrate the method by analyzing the proton exchange reactions in formic acid solvated in small water clusters. The simulations were performed with ab initio MD combined with a method to efficiently sample the rare event, path sampling. Our ML analysis identified relevant geometric variables involved in the proton transfer reaction and how they may change as the number of solvating water molecules changes.

Force-Field-Based Computational Study of the Thermodynamics of a Large Set of Aqueous Alkanolamine Solvents for Post-Combustion CO2 Capture

The ability to predict the thermodynamic properties of amine species in CO₂-loaded aqueous solutions, including their deprotonation (pK_a) and carbamate to bicarbonate reversion (pK_c) equilibrium constants and their corresponding standard reaction enthalpies, is of critical importance for the design of improved carbon capture solvents. In this study, we used isocoulombic forms of both reactions to determine these quantities for a large set of aqueous alkanolamine solvent systems. Our hybrid approach involves using classical molecular dynamics simulations with the general amber force field (GAFF) and semi-empirical AM1-BCC charges (GAFF/AM1-BCC) in the solution phase, combined with high-level composite quantum chemical ideal-gas calculations. We first determined a new force field (FF) for the hydronium ion (H₃O⁺) by matching to the single experimental pK_a data point for the well-known monoethanolamine system at 298.15 K. We then used this FF to predict the pK_a values for 76 other amines at 298.15 K and for all 77 amines at elevated temperatures. Additionally, we indirectly relate the H₃O⁺ hydration free energy to that of H⁺ and provide expressions for intrinsic hydration free energy and enthalpy of the proton. Using the derived H₃O⁺ FF, we predicted the pK_a values of a diverse set of alkanolamines with an overall average absolute deviation of less than 0.72 pK_a units. Furthermore, the derived H₃O⁺ FF is able to predict the protonation enthalpy of these amines when used with the GAFF. We also predicted the carbamate reversion constants of the primary and secondary amine species in the data set and their corresponding standard heats of reaction, which we compared with the scarcely available experimental data, which are often subject to significant uncertainty. Finally, we also described the influence of electronic and steric effects of different molecular fragments/groups on the stabilities of the carbamates.

Ab Initio Molecular Dynamics Simulations of the Influence of Lithium Bromide Salt on the Deprotonation of Formic Acid in Aqueous Solution

The deprotonation of formic acid is investigated using metadynamics in tandem with Born-Oppenheimer molecular dynamics simulations. We compare our findings for formic acid in pure water with previous studies before examining formic acid in aqueous solutions of lithium bromide. We carefully consider different definitions for the collective variable(s) used to drive the metadynamics, emphasizing that the variables used must include all of the possible reactive atoms in the system, in this case carboxylate oxygens and water hydrogens. This ensures that all the various possible proton exchange events can be accommodated and the collective variable(s) can distinguish the protonated and deprotonated states, even over rather long ab initio simulation runs (ca. 200-300 ps). Our findings show that the formic acid deprotonation barrier and the free energy of the deprotonated state are higher in concentrated lithium bromide, in agreement with the available experimental data for acids in salt solution. We show that the presence of Br^- in proximity to the formic acid hydroxyl group effectively inhibits deprotonation. Our study extends previous work on acid deprotonation in pure water and at air-water interfaces to more complex multicomponent systems of importance in atmospheric and marine chemistry.

The multiple dissociation constants of glutathione disulfide: interpreting experimental pH-titration curves with ab initio MD simulations

The hexapeptide glutathione disulfide (GSSG) has six ionizable groups with six associated dissociation constants. The experimentally measured pH-titration curve, however, does not exhibit the six corresponding equivalence points and bears little resemblance to standard textbook examples of acid-base pH-titration curves. The curve highlights the difficulties in determining multiple pKa values of polyprotic acids - typically proteins and peptides - from experiment. The six pKa values of GSSG can, however, be estimated using Car-Parrinello molecular dynamics (CPMD) simulations in conjunction with metadynamics sampling of the underlying free energy landscape of the dissociation reactions. Ab initio MD simulations were performed on a GSSG molecule solvated by 200 water molecules. Using the estimated pKa values the theoretical titration curve was calculated and found to be in good agreement with experiment. The results clearly highlight how dissociation constants estimated from ab initio MD simulations can facilitate the interpretation of the pH-titration curves of complex chemical and biological systems.

A computational scheme of pKa values based on the three-dimensional reference interaction site model self-consistent field theory coupled with the linear fitting correction scheme

A scheme for quantitatively computing the acid dissociation constant, pKa, of hydrated molecules is proposed. It is based on the three-dimensional reference interaction site model self-consistent field (3D-RISM-SCF) theory coupled with the linear fitting correction (LFC) scheme. In LFC/3D-RISM-SCF, pKa values of target molecules are evaluated using the Gibbs energy difference between the protonated and unprotonated states calculated by 3D-RISM-SCF and the parameters fitted by the LFC scheme to the experimental values of training set systems. The pKa values computed by LFC/3D-RISM-SCF show quantitative agreement with the experimental data.

Development of a ReaxFF reactive force field for lithium ion conducting solid electrolyte Li1+xAlxTi2−x(PO4)3 (LATP)

Using a ReaxFF reactive force field, we investigated the composition-dependent ionic conductivity and the Li migration behaviors in Li_1+xAl_xTi_2−x(PO₄)₃ solid electrolyte.

Microhydration of a benzoic acid molecule and its dissociation

C₆H₅COOH requires at least eight H₂O molecules for its dissociation. A linear correlation is observed between pK_a and n_H2O; an unknown pK_a can thus be predicted.

Linear Response Function of Bond-Order

We present the linear response function of bond-orders (LRF-BO) based on a real space integration scheme for molecular systems. As in the case of the LRF of density, the LRF-BO is defined as the response of the bond order of the molecule for the virtual perturbation. Our calculations show that the LRF-BO enables us not only to detect inductive and resonating effects of conjugating systems, but also to predict pK_a values on substitution groups via linear relationships between the Hammett constants and the LRF-BO values for meta- and para-substituted benzoic acids. More importantly, the LRF-BO values for the O-H bonds strongly depend on the sites to which the virtual perturbation is applied, implying that the LRF-BO values include essential information about reaction mechanism of the acid-dissociation of substituted benzoic acids.

Ab Initio MD Simulations of the Brønsted Acidity of Glutathione in Aqueous Solutions: Predicting pKa Shifts of the Cysteine Residue

The tripeptide glutathione (GSH) is one of the most abundant peptides and the major repository for nonprotein sulfur in both animal and plant cells. It plays a critical role in intracellular oxidative stress management by the reversible formation of glutathione disulfide with the thiol-disulfide pair acting as a redox buffer. The state of charge of the ionizable groups of GSH can influence the redox couple, and hence the pKa value of the cysteine residue of GSH is critical to its functioning. Here we report ab initio Car-Parrinello molecular dynamics simulations of glutathione solvated by 200 water molecules, all of which are considered in the simulation. We show that the free-energy landscape for the protonation-deprotonation reaction of the cysteine residue of GSH computed using metadynamics sampling provides accurate estimates of the pKa and correctly predicts the shift in the dissociation constant values as compared with the isolated cysteine amino acid.

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.

Cluster-continuum quasichemical theory calculation of the lithium ion solvation in water, acetonitrile and dimethyl sulfoxide: an absolute single-ion solvation free energy scale

A solvation free energy scale excluding the net electrostatic potential inside the solute cavity is presented.