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.

Introduction

Acid–base chemistry is a crucial aspect of all chemical fields. Given this obvious importance, it is remarkable that so little is understood about the microscopic mechanisms underlying acid–base reactions. Progress in this task must involve investigations into two related but somewhat separate questions. First, we need to understand the initial bond-breaking mechanism(s) which lead to the creation of solvated protons and hydroxyl anions. Second, once the free H+ and/or OH– ions are produced, the state of these solvated ions in aqueous solution must be better understood.

In this article, we report on our work using ab initio molecular dynamics (AIMD) simulations to investigate deprotonation of formic acid and the resulting solvated proton. Because of the bond-breaking process, empirical models are not well-suited to study the hydrogen dissociation reaction. Although empirical reactive force fields such as ReaxFF[1] show some promise, in our experience, they are not yet accurate enough to study this problem in a quantitative way. Previous work using AIMD simulations has made some progress in both modeling the deprotonation and the resulting state of the hydrated proton. Using metadynamics to drive the trajectory over the reaction barrier, some studies have been published which can simulate the initial deprotonation in various aqueous acids.[2−10] In particular, some interesting insights into the difference between the behavior of an acid molecule located on the air–water interface versus bulk solution have been obtained.[5,6,10]

As for the nature of the hydrated proton, AIMD simulations and other theoretical approaches have been used to good effect.[11−16] The modern understanding of proton hydration is that the “free” proton exists in a dynamic equilibrium between Eigen (H3O+·(H2O)3) and Zundel (H5O2+) cations. Very recent work using path sampling approaches has reignited debate regarding the mechanisms involved in water autoionization.[11,16]

We study formic acid in pure water, as well as in an aqueous solution of lithium bromide. Our choice of lithium bromide was motivated by our previous work[10,17,18] to model interesting experimental results involving molecular collisions with liquid microjets,[19] where highly concentrated lithium bromide solutions are used to mitigate against water evaporation from the surface. To our knowledge, the influence of an alkali halide salt on the deprotonation of acid has never been studied in a molecular dynamics simulation. There are some experimental data on, for example, carbonic acid in artifical seawater,[20−22] and other organic acids including formic acid in different salt solutions,[23,24] showing the dependence on acidic pH as a function of salt concentration. However, a molecular understanding of this dependence is lacking. An AIMD study may lead to insights into, for example, whether the presence of salt ions may aid or suppress the formation of hydrogen bond networks believed to be important in acid deprotonation mechanisms and may also be relevant to a better understanding of the crucial issue of ocean acidification, where the influence of salt may need to be accounted for to precisely model carbonic acid in seawater.[4,5,25]

Computational Details

Initial configurations for the AIMD trajectories were generated with LAMMPS[26] using the TIP4P-Ew empirical water model,[27] the associated Joung–Cheatham models for the Li+ and Br– ions,[28] and a model for the formic acid molecule due to Jedlovszky and Turi.[29,30] Because we used the SHAKE algorithm in LAMMPS to maintain the rigid geometry of water, we were unable to maintain the rigidity of the entire formic acid molecule consistent with the Jedlovszky–Turi model. Instead, a strong harmonic bond and angle potentials as well as two improper dihedral potentials from the OPLS-AA model were used for the intramolecular interactions in formic acid, except for the OH bond which was kept rigid with SHAKE along with the water geometry.

After initial equilibration with the empirical models over several nanoseconds of simulation, AIMD simulations were started. These Born–Oppenheimer molecular dynamics simulations were run using the Quickstep module of CP2K[31] with a 0.5 fs timestep. The BLYP exchange correlation functional[32,33] was used, supplemented with Grimme’s D2 dispersion correction[34,35] which has been shown to be important in simulations of liquid aqueous systems. For the basis sets, we used the DZVP-MOLOPT-SR-GTH basis functions along with the BLYP-GTH pseudopotentials.[36] Overall, the properties of bulk water are modeled well by this combination of basis set and pseudopotentials and dispersion correction,[37,38] with the exception of an overestimated equilibrium density.[10,18,39] Based on our previous experience with aqueous lithium bromide[18] and with aqueous formic acid,[10] we are confident that this choice of simulation parameters should be sufficiently accurate for our systems and purposes.

The temperature T was held constant at 300 K by using a Nose–Hoover chain thermostat of length 3 and time constant 50 fs applied on the massive degrees of freedom.[40−42] All of the simulations were performed in a fixed cubic simulation geometry with three-dimensional periodic boundary conditions. When salt was added, the simulation box was modified to accommodate the higher equilibrium density of the salt water. The simulation box lengths L are given along with the rest of the simulation parameters in Table . We did not attempt to run simulations in the NPT ensemble, which could determine the “correct” system density at zero pressure. As mentioned above, with our chosen level of theory and basis set, we expect the equilibrium density of water to be significantly overestimated.[10,18,39]NPT trajectories must also be long enough to assure statistical convergence, and this would be challenging for our expensive AIMD simulations. Instead, we chose to fix the system size to match the experimental density of these systems.

Table 1

Summary of Simulation Parameters, Including Collective Variables and Metadynamics Parametersa

nH2O	nLiBr	L/Å	tsim/ps	# traj.	n	m	r0/Å
63	0	12.4	200	4	6	18	1.4
57	3	12.3	300	4	6	18	1.3
43	10	12.3	300	5	6	18	1.3

Columns 1 and 2 are the numbers of water molecules and lithium bromide ion pairs, respectively. Column 3 is the box length in Å = 10–10 m. Columns 4 and 5 show the total length of each individual trajectory analyzed, and the number of trajectories run for each system. Columns 6–8 show the parameters used in the collective variable function (eq ), that is, the exponents n and m and the cutoff radius r0

An initial equilibration run over 20–40 ps was completed without metadynamics to ensure that the energetics and structure of the AIMD simulations had converged. After that, we moved on to determining the free energy landscape of formic acid deprotonation using the built-in capability of CP2K to run well-tempered metadynamics.[43−45] In this method, the interatomic potential energy U0() is modified by the addition of a bias potential V(,t) which depends on one or more collective variables () = (s1, s2, ..., s). The bias potential consists of the addition of repulsive Gaussians placed along the collective variable(s) coordinates as the simulation proceeds. The precise form of the bias potential iswhere τG represents the time interval between deposited Gaussians, and ′ are respectively the instantaneous and previous values of the collective variables when Gaussians were deposited, σ is the width of the deposited Gaussians in each collective variable i, and W is the height of the Gaussians. Unlike in standard metadynamics, in well-tempered metadynamics, W is not constant but is history-dependent as well,

In this expression, ΔT is a bias temperature which is typically a few times larger than the system temperature T. As the simulation proceeds, the height lowers toward a constant value W0. The reduction in the height of the Gaussians with increasing time guards against overestimation of the free energy during long metadynamics simulations. In our work, we use a value of τG = 125 fs and single values of σ = 2.6255 kJ mol–1, and of W0 = 0.02, and we apply a bias temperature ΔT = 1800 K.

An important consideration in using metadynamics to compute free energy barriers is ensuring that the collective variables () = (s1, s2, ..., s) are well chosen.[46] They must allow the simulation to be biased effectively in order to explore all of the relevant phase space. Many functional forms for the collective variables are possible, subject to the limitation that the variables must be continuous and differentiable. In much of the previous work on formic acid or other organic acids,[7−10] a single collective variable s1 was used which approximates the coordination number between the protonated carboxylate oxygen in formic acid and its bound hydrogen atom,where rOH is the distance between the oxygen and hydrogen atoms, the parameter r0 is chosen to distinguish between the protonated and deprotonated states, and the exponents n and m determine the sharpness of the change of the function from 1 for rOH ≪ r0 down to 0 for rOH ≫ r0. Typical values for these parameters used in previous work are n = 6, m = 12, and r0 = 1.6 Å.

A function of the form of eq does work well for determining whether the initial acidic proton has broken the bond to the carboxylate oxygen. However, there are issues with this simple collective variable definition. In an ab initio simulation, in principle, there can be no distinction made between any particular atoms of the same element in the system. In a long simulation, the possibility that the initial acidic proton is replaced by a new hydrogen atom bound to either of the two carboxylate oxygens must be accounted for in the definition of the collective variables. Otherwise, a small value of the collective variable may not in fact represent a deprotonated state.

In this study, we have modified the simple expression for s1 above, in an attempt to better describe the conformational space of protonation and deprotonation in formic acid. Instead of only the initial formic acid OH bond, we consider all possible OH atom pairs, writing s2 now aswhere H and O indicate each reactive hydrogen and each carboxylate oxygen atom, respectively. When all hydrogens are included, a new complication is that they all will contribute a small amount to s2, even though the interatomic distance may be well above r0. In particular, the contribution of water molecules hydrogen bonded to the formic acid can lead to values of s2 well in excess of 1. We modify some of the equation parameters, using a larger value of m = 18 and smaller value of r0 to try to avoid these extra contributions to s2.

The definition of s2 we use does not allow us to explicitly distinguish protonation of the different oxygen atoms. To allow this distinction, another approach could be to use two different collective variables, one for each carboxylate oxygen.[3] Other, more complicated schemes with as many as three collective variables have been introduced in previous work.[2] However, the simplicity of a single collective variable is a considerable advantage of our approach. Using multiple collective variables requires much longer simulations to ensure complete exploration of the free energy landscape.

Results and Discussion

We studied formic acid in pure water and in two different concentrations of aqueous lithium bromide. Long simulations and multiple trajectories are required to make reliable estimates of deprotonation barriers and to model the deprotonated state. The details of our simulations are given in Table .

In Figure , we plot all of the free energy surfaces (FESs) generated by the metadynamics simulations for all three of our systems. In all trajectories, we see that the metadynamics biasing is able to force the initial OH bond in the formic acid to stretch beyond the cutoff radius of r0 = 1.3 Å indicated by the exploration of small values of the collective variable s2 defined in eq . However, there is tremendous variability in the likelihood of the trajectory investigating these small values. In many trajectories, a small value of s2 seems like a reasonable indication that the simulation is authentically representing a deprotonated state. In others, s2 only briefly remains small before either the initial OH bond length becomes smaller again, or there is a proton exchange event where a different proton forms a new OH bond to either of the carboxylate oxygens.

Figure 1

Deprotonation free energy landscape of formic acid calculated via metadynamics simulations as a function of the collective variable s2 (see eq ), relative to the protonated state. Left side (insets): Individual trajectories (thin lines) and average over all trajectories (thick line) computed for (a) pure water, (b) low concentration LiBr solution, and (c) high concentration LiBr. Right side: (d) Comparison of average free energy landscapes in all three systems. Pure water is the solid orange line, low concentration LiBr is the dash-dotted red line, and high concentration LiBr is the dashed black line. Error bars are one standard error in the mean of all trajectories for that system (see Table ).

In Table , we show our best estimates for the barrier height ΔGbarrier relative to the protonated state, and the free energy difference ΔGprot between the protonated and deprotonated states. ΔGprot can then be converted to estimate the pKa value of the acid,

Table 2

Barrier Heights, Free Energy Difference between the Protonated and Deprotonated States, and Estimated pKa for formic acid

solvent/source	ΔGbarrier/kJ mol–1	ΔGprot/kJ mol–1	pKa
pure water	14.8 ± 0.8	8.0 ± 4.4	1.4 ± 0.8
pure water, ref (10)	14.8
pure water, ref (3)	17.2
pure water, ref (7)		22.2	3.86
pure water, expt.		21.7	3.78
low conc. LiBr	18.4 ± 0.8	13.4 ± 3.4	2.3 ± 0.6
high conc. LiBr	27.1 ± 2.0	22.8 ± 7.5	4.0 ± 1.3

Our results for ΔGbarrier in pure water are in good agreement with some of the previous results obtained via metadynamics.[3,10] However, there is a clear discrepancy between the barrier height and the free energy difference between protonated and deprotonated states ΔGprot predicted by the experimental pKa, in that our prediction of ΔGbarrier < ΔGprot. Only the work of Tummanapelli and Vasudevan correctly matches the experimental pKa,[7] but as described above there may be issues with their definition of the collective variable which may make their estimate of the free energy landscape unreliable.

Figure and Table also show the effect of lithium bromide on the free energy barriers according to our simulations. We see that as the salt concentration increases, both ΔGbarrier and ΔGprot increase. This is generally in agreement with available experimental data, which show that weak acid ionization is enhanced by very low salt concentration (<0.5–1.0 M) before being inhibited at higher concentrations.[23,24]

To gain some more insight into the details of the trajectories, it is helpful to do some more analysis. In Figures a and 4a,b, we show the results of tracking a number of different interatomic distances. A schematic representing which distances are shown in which colors is shown in Figure . Only the initial formic acid OH bond length shown in black refers to a specific pair of atoms. All of the other interatomic distances refer to different atom pairs over the course of the simulation. This allows us to track the most relevant interatomic distances, which are usually the minimum distances between two particular atom types. Figures b and 4c show two more quantities for the same trajectories. The quantity s2 is the collective variable defined in eq , while Fprot is a flag which tracks whether formic acid is protonated (Fprot = 1), or if the acidic proton should be considered to be hydrated. If Fprot = 2, 3, or 4, there is a water oxygen with three hydrogen atoms closer than 1.3 Å, one of which is the acidic proton. The flag Fprot = 2 if the acidic proton is still closely associated with the formic acid (rOH < 1.3 Å). The flag Fprot = 3 if the hydrated proton is best described as an Eigen cation (H3O+·(H2O)3), while Fprot = 4 if the hydrated proton is part of a Zundel cation (H5O2+). Finally, rare configurations with Fprot = 5 define states where the acidic proton is bound neither to the formic acid nor to a water oxygen. States with Fprot = 5 can be regarded as transient states where our simple geometric criteria defining the status of the proton are inadequate.

Figure 3

Tracking of some relevant time-dependent quantities from an AIMD metadynamics trajectory for formic acid in pure water, with the FES displayed as the red line in Figure a. (a) Indicated distances rOH, including the shortest distance from a formic acid oxygen to the hydrated proton (violet points). (b) Collective variable s2 used to drive the metadynamics simulation and the flag Fprot indicating the state of the hydrated proton (see text).

Figure 4

Tracking of some relevant time-dependent quantities from an AIMD metadynamics trajectory for formic acid in concentrated LiBr solution, with the FES displayed as the red line in Figure c. (a) Indicated distances rOH, including the shortest distance from a formic acid oxygen to the hydrated proton (violet points). (b) Distances rXI from the indicated formic acid atoms X to the indicated ion I. (c) Collective variable s2 used to drive the metadynamics simulation and the flag Fprot indicating the state of the hydrated proton (see text).

Figure 2

Schematic showing colors of different interatomic distances plotted in Figures and 4.

Figure analyzes one of the trajectories of formic acid in pure water, with the resultant free energy function displayed as the red line in Figure a. Notable events along this trajectory are indicated by vertical dashed lines and a short description. Around t = 20 ps, we see a clear and fast proton exchange event, where the initial acidic proton dissociates from one carboxylate oxygen and is immediately replaced by a different proton on the other carboxylate oxygen. The next interesting event occurs at around t = 105 ps, where a truly deprotonated formate ion is formed. Hydrogen bonds with distances less than 2 Å still exist between water hydrogens and carboxylate oxygens, but the acidic proton is consistently found at a distance ∼5–7 Å from the formate ion for approximately 40 ps. During this time period, rapid fluctuations between Eigen and Zundel cations are detected. At t ≃ 145 ps, another proton exchange event sees a new proton form a bond to the initially protonated oxygen and the formic acid molecule is again protonated.

Figure analyzes one of the trajectories in concentrated aqueous LiBr, with the free energy function shown as the red line in Figure c. Just as in the pure water trajectory, we see a clear deprotonation event around t = 125 ps, followed by a reprotonation event around t = 175 ps. Again this reprotonation event involves a proton exchange where a proton originally in a water molecule now protonates the initially unprotonated carboxylate oxygen. During the lifetime of the formate anion, the hydrated proton is located quite far from the formate (>5 Å). By tracking the distances from the formic acid atoms to the salt ions (Figure b), we can attempt to understand their influence on the proton exchange dynamics. We can see that at around t = 55 ps, a fluctuation occurs whereby a Br– ion shifts from being located close (ca. 2 Å) to the acidic proton to a further separation. After this event, it is clear that larger fluctuations of the formic acid OH bond length become possible, until eventually the formic acid deprotonates. In other words, when a bromine ion is in close contact with the acidic proton, its presence inhibits large amplitude motions of the acidic proton required to undergo deprotonation. We also note that a higher likelihood of finding states with Fprot = 5, compared with the pure water simulations, is often due to close associations between the free proton and a Br– ion, instead of a water oxygen or the carboxylate group. This can in particular be seen in the early part of the trajectory before the closest Br– ion moves away from the vicinity of the formate (see Figure c).

Another way to investigate the possible influence of the bromine ion on deprotonation is to include it in the metadynamics scheme by introducing a second collective variable s2,Br, which has a similar functional form to that used for s2, but, instead of summing over the hydrogen atoms, we sum over the bromine ions,

The value of r0,Br is also different, so that the collective variable can distinguish configurations with a Br– ion in the hydration shells of the formic acid oxygen atoms. By experimentation, we found a value of r0,Br = 3.9 Å to work well.

In Figure , we show the two dimensional FES resulting from a metadynamics simulation using both s2 and s2,Br as collective variables. The added dimension makes it considerably more challenging to run a trajectory sufficiently long to adequately explore the complete FES. Another complication is the difficulty of clearly distinguishing a single close contact between Br– and formic acid oxygens, leading to values of s2,Br significantly larger than 1 due to contributions to s2,Br from bromine ions which may be well outside of the first solvation shell. Nevertheless, we see that when s2,Br is large, states with s2 being small are never seen. In other words, when Br– is in close contact with protonated formic acid, its presence inhibits the ability of the proton to leave the formic acid and participate in the hydrogen bond network in such a way that Eigen or Zundel cations might be formed.

Figure 5

Contour plot of the FES (in kJ mol–1, relative to the protonated formic acid) for formic acid deprotonation with two collective variables, one for the O–H distance (s2, horizontal axis) and one for the O–Br– distance (s2,Br, vertical axis), using a long metadynamics simulation of ∼450 ps duration. The contour lines are separated by 8 kJ mol–1.

Conclusions

We have completed a computational study of formic acid deprotonation using metadynamics with AIMD simulations. We carefully consider different definitions of collective variable(s) and argue that the collective variable scheme used must at least be able to describe proton exchange events involving both carboxylate oxygens as well as all reactive protons. We have performed comparatively long AIMD simulations (200–300 ps) and multiple trajectories to ensure the trustworthiness of our results. Our estimate of the free energy barrier for deprotonation, ΔGbarrier = 14.8 ± 0.8 kJ mol–1, is close to two previous estimates using similar density functional theory (DFT) methodology.[3,10] However, there is still some work to be done to match experimental results because the experimental formic acid pKa = 3.78 requires a free energy difference ΔGprot = 21.7 kJ mol–1, which is larger than our estimate of the barrier height. The only simulation results which agree with the experimental data are those of Tummanapelli and Vasudevan.[7] However, as we have discussed above, their use of a single OH atom pair to define the reaction coordinate, and short metadynamics trajectories, makes it difficult to view their results as predictive.

All things considered, we contend that our results for ΔGbarrier are an accurate estimate for this combination of density functional and basis sets, being based on long simulations and collective variables which can accommodate all possible proton transfer events between the carboxylate group and reactive protons. Discrepancies that remain between the experimental data and the simulation results can likely be assigned to the limitations of the DFT methodology we have used. For example, the use of a triple-zeta basis set would tend to increase the barrier height significantly and bring our results into closer agreement with the experimental data.[3,10] As far as the ΔGprot and resultant estimate of the pKa are concerned, it is clear that our simulation results show large variability over multiple independent trajectories. Overall, we must conclude that our collective variable definition is not completely able to distinguish the protonated acid from the deprotonated state.

We find that the addition of lithium bromide salt has the effect of significantly increasing both the free energy barrier and the energy difference between protonated and deprotonated states. Our results are in qualitative agreement with the experimental data, which shows an increase in pKa in concentrated ionic solutions compared with acid in pure water.[23,24] We have shown some trajectory analyses and supplementary metadynamics simulations, suggesting that the presence of a large bromine ion in close contact with the acidic proton strongly inhibits the pathways for deprotonation. Although we believe our findings regarding the effect of salt ions should be rather general, additional studies of different ionic species would be worthwhile in order to proceed further in our understanding of how salt ions affect acidity.

Notwithstanding our difficulties in achieving quantitative agreement with the experimental data, one of our goals in completing this study was to carefully consider and critique previous attempts to use metadynamics to study deprotonation and point out that collective variable schemes which are fundamentally unable to include, for example, proton transfers between different carboxylate oxygens, cannot fully describe deprotonation in carboxylic acids. Similar issues are likely to be important in describing other proton transfers in other organic acids as well as in amines.[7−9,47] Our position is that we have done as well as can be done with metadynamics simulations, in particular using a simple one-dimensional collective variable (s2).

Further advancement in the understanding of deprotonation reactions, and more generally, the nature of solvated protons, may benefit from new theoretical approaches. A key limitation of the metadynamics method is the requirement that the collective variable(s) defined to fit the reaction pathway must be continuous functions. This leads to requiring rather complicated functions if one hopes to be able to describe, for example, the distance from the acid to the solvated proton.[2] As we have already seen, while simple, our single collective variable approach seems to be inadequate for describing the deprotonated state when the solvated proton is far from the acidic anion. Another approach to rare event simulation is transition path sampling, and related more advanced methods such as transition interface sampling (TIS), where Monte Carlo methods are used to bias the trajectories themselves in order to explore the reaction pathway, rather than altering the potential energy surface.[48,49] Another key advantage is that it allows the use of discontinuous collective variables. It seems likely that we can take a cue from the approach already shown to be useful to model water autoionization[11,16] and apply it to the problem of acidic deprotonation. Recent work with TIS has led to new insights into the mechanisms of deprotonation and aqueous proton transport.[16] We are currently experimenting with the PyRETIS library[50] to implement this novel approach, and we look forward to reporting on our progress soon.