Acidities of MgO surface sites: implications for the formation mechanism of Mg(OH)₂

Abstract

The hydroxylation of periclase (MgO) to brucite (Mg(OH)₂) is thought to be an important intermediate step when using MgO to capture CO₂ from the atmosphere. However, the mechanism of hydroxylation of MgO to form Mg(OH)₂ is poorly understood. In this work, we used atomic-scale density functional tight binding simulations coupled with the metadynamics rare event method to analyze the surface chemistry of MgO and the acid dissociation equilibrium constants (pK_a) of its surface sites. The method and parameters were validated by calculating the pK_a for hydroxylation of the first shell water bound to aqueous Mg²⁺ ion. The pK_a value derived using a probabilistic method was 12.3, which is in fair agreement with the accepted value of 11.4, with the difference between them equal to a ∼5 kJ mol⁻¹ error in the calculations. We then extended these pK_a calculations to probe the hydroxylation reactions of the surface sites of the MgO(100)–water interface, arriving at pK_as of 5.4 to deprotonate terminal water molecules bound to the surface magnesium sites (η-OH₂ or 〉MgOH₂), and 13.9 to deprotonate hydroxylated bridging oxygen sites (μ₅-oxo or 〉O). Hydroxide (OH⁻) adsorption on the surface was also probed and found to be less thermodynamically favorable than deprotonation of the terminal water molecule. The plausibility of the computed pK_as was verified using an activity-based speciation model and compared to pH measurements of water equilibrated with MgO nanoparticles and single crystals. The model predicted a solution pH of 7.1 when surface sites buffered and the pH of 12.0 when MgO dissolution dominated. These are close to the experimental initial solution pHs of 7–7.5 and the long term pHs of ∼10.5. The similarity suggests that the calculated pK_a values from the DFTB+/metadynamics simulations are plausible and that these methods can be a useful tool to probe reaction mechanisms involving covalent bonds.

---

1 Introduction

A potential strategy for large-scale direct air capture (DAC) of CO₂ is to use mineral-sourced alkaline earth metal oxides in a looped process.¹ In this process, the oxide sorbent reacts with CO₂ in air and forms carbonate phases which are then calcined to separate the CO₂ and regenerate the sorbent. The process is then repeated, creating a loop. The captured CO₂ can be sequestered² and used for producing synthetic aggregates for concrete,^3–5 for enhanced oil recovery^6,7 and in food processing.⁸ DAC mineral looping technologies based on conversion of periclase (MgO) or quicklime (CaO) may be a cost effective method to separate CO₂ from the atmosphere.^9–11 While most of the industrial attention to date has focused on CaO, magnesite (MgCO₃) has a lower calcination temperature than calcite (CaCO₃). This lower regeneration temperature property could make MgO more economical for mineral looping DAC of CO₂ (ref. 12) if carbonation rates were similar, but the carbonation rate of MgO is thought to slow with time due to an armoring effect of the carbonate phases.¹³ Thus, there is a compelling reason to work to understand the mechanisms of carbonation of MgO.

In the techno-economic analysis of McQueen et al.,¹² it is assumed that within the mechanism of carbonation there is first a facile conversion of MgO to brucite (Mg(OH)₂). It is then assumed that the nucleation of magnesium carbonate phases is fast, and thus the dissolution of Mg(OH)₂ controls the overall rate of carbonation. Consistent with this interpretation, analysis of the carbonation of MgO has shown d-spacings consistent with Mg(OH)₂, prior to forming hydrous magnesium carbonate phases.¹³ It is widely accepted that the conversion of MgO to Mg(OH)₂ occurs by a dissolution and precipitation mechanism consisting of the following reaction steps:^14–17

Hydroxylation: MgO_(s) + H₂O_(l) → MgOH⁺_(surf) + OH⁻_aq (1)

Adsorption: MgOH⁺_(surf) + OH⁻_aq → MgOH⁺_(surf)·OH⁻ (2)

Dissolution: MgOH⁺_(surf)·OH⁻ → Mg²⁺_aq + 2OH⁻_aq (3)

Precipitation: Mg²⁺_aq + 2OH⁻_aq → Mg(OH)_2(s) (4)

where the suffixes s, l, surf and aq represent solid, liquid, surface and aqueous solution respectively.

In the above mechanism, interfacial water molecules dissociate and hydroxylate the surface (eqn (1)). The OH⁻_aq ions from the dissociated water molecules then adsorb on the positively charged hydroxylated surface (MgOH⁺_(surf), eqn (2)) and release Mg²⁺_aq ions from the surface by dissolution (eqn (3)). The Mg²⁺_aq ion concentration and solution pH increase (≥10) until a critical supersaturation is reached, where Mg(OH)₂ solid precipitates (eqn (4)). Mg(OH)₂ formation after the dissolution step may vary slightly depending on whether homogeneous or heterogeneous nucleation occurs.^14,15

Neither the reaction mechanism nor which surface site it is that actually hydroxylates is known. Of the MgO surfaces, the (100) has a lower surface energy than either the (110) or (111) surfaces,¹⁸ suggesting it will be the most common surface. On it, there are three likely possibilities for surface sites that could hydroxylate: the magnesium surface sites with bound water (η-OH₂ or 〉MgOH₂〉), magnesium sites without anything bound at all (〉Mg) that have been observed computationally adjacent to hydroxylated bridging oxygens (μ₅-hydroxo),^14–17 or bridging oxygen sites (μ₅-oxo or 〉O). (Note: for clarity, formal charges are used here, but in actuality the Mg surface sites will accept some electron density from neighboring oxygens that reduces the divalent positive charge, and vice versa for the μ₅-oxo sites, which will reduce their divalent negative charge.)

Surface hydroxylation: 〉MgOH²⁺₂+OH⁻ ⇌ 〉MgOH⁺ + H₂O (5)

〉Mg²⁺ + OH⁻ ⇌ 〉MgOH⁺ (6)

〉O²⁻ + H₂O ⇌ 〉OH⁻ + OH⁻ (7)

It is not clear which of these moieties are critical for dissolution to occur; however, several studies have assumed that the non-hydrated magnesium sites are the more likely candidate.^14–17 The goal of this study is to determine the acid dissociation equilibrium constants (pK_a) of these sites as presented in eqn (5)–(7) to identify which of these sites is most likely to undergo hydroxylation.

To define the MgO(100)–water interfacial environment in which these reactions occur, a monolayer of adsorbed, surface-bound water molecules has been observed in both experimental and computational studies.^19–23 In the experiments and density functional theory (DFT) calculations of Hollerer et al.²⁴ and Włodarczyk et al.²⁵ on the MgO(100) surface at sub-ambient conditions with the highest water coverage (p(3 × 2) − 1.33/uc and c(4 × 2) − 1.25/uc) also show the presence of both surface hydroxyls and OH⁻_aq ions with their axes oriented parallel to the surface normal. These OH⁻_aq ions are located slightly above the surface-bound water monolayer. DFT simulations by Sassi and Rosso²³ showed that the p(3 × 2) − 1.33/uc of the MgO(100)–water system is the energetically most favorable structure at ambient conditions, and that 25% of the water molecules bound to surface magnesium sites are dissociated. Their simulations also confirmed the retention of the water monolayer structure and protrusion of OH⁻_aq ions slightly above the water layer at ambient conditions.²³ A similar structure of ordered, interfacial water molecules and OH⁻_aq ions located slightly further from the interface than the terminal water layer was also observed in a longer simulation (∼35 ps) using ab initio molecular dynamics (AIMD).²²

Although all of the available experimental and computational studies have confirmed the likelihood of MgO surface hydroxylation (eqn (1)), none has explored the relative favorabilities of the different surface sites, nor the OH⁻_aq ion adsorption step in the mechanism (eqn (2)). Recent progress in a quantum chemical technique, semi-empirical Density Functional Tight Binding (DFTB), coupled with the metadynamics rare event method facilitates the calculation of proton transfer energetics and OH⁻_aq adsorption reactions, which means that a better understanding of these reaction mechanisms is now possible. We have used these methods to calculate the energetics of proton transfer (ΔF, thus pK_a) reactions. The known pK_a of water was used for calibration and that of the aqueous Mg²⁺ ion was used to validate the simulation techniques, and the system parameters (Mg(H₂O)²⁺_6aq → Mg(H₂O)₅OH⁺_aq + H⁺_aq). To confirm that the measured pK_a values were plausible, we compared pHs predicted by speciation software using activity-concentration relationships and surface charge to the pH observed in aqueous solutions in contact with MgO nanoparticles and, second, larger suspensions of a crushed single crystal.

2 Methodology

2.1 Simulation

All quantum chemical molecular dynamics simulations were performed with energies and forces computed directly based on the DFTB method, as implemented in the DFTB+ software package.²⁶ DFTB is an approximate density functional theory (DFT) method that is between two and three order of magnitude faster than conventional DFT calculations.²⁷ It thereby allows routine simulations on longer time and length scales, crucial for the simulations presented in this work that require ensemble average solvation structures and a biased simulation. Interatomic Coulomb interactions were treated with the third-order self-consistent charge DFTB method (i.e., DFTB3) in conjunction with the 3OBw Slater-Koster parameter set.^28,29 At each molecular dynamics step, iteratively partial charges get updated until a pre-defined accuracy is achieved. Atomic partial charges are calculated using Mulliken charge populations at each time step. DFTB has limitations in providing correct vibrational frequencies, whereas the improved version, DFTB3, accurately describes structures and binding energies.^30–34 Prior work on pK_a estimation using DFTB includes inorganic aqueous systems and biological systems. The Gamma point approximation was used in all periodic DFTB calculations. To bias the simulation to explore deprotonation and adsorption reactions, we employed DFTB-based metadynamics using multiple walkers as implemented in the PLUMED2 (ref. 35 and 36) library interoperating with the DFTB+ software.²⁶ The initial system configurations used in all DFTB simulations were pre-equilibrated using classical molecular dynamics simulations, following the approach mentioned in Adapa et al.³⁷

Using DFTB with metadynamics, we first simulated an aqueous system containing 1 Mg²⁺ ion, 2Cl⁻ ions and 750 water molecules in a cubic box with a volume (V) of 7495 ų to verify that we could reproduce the known bulk structural properties (g(r) and Cn) and pK_a (note: eqn (7) is written as a protonation). For surface calculations, the unit cell parameters of MgO reported by Hazen³⁸ were then used to create a flat, four layer MgO(100) surface, which had the dimensions of 12.64 × 12.64 × 8.43 ų (144 slab atoms, 4 monolayers thick) and was in contact with a 20 Å thick layer of pre-equilibrated water (110 water molecules) in the surface normal (z−) direction. This system was periodic in three dimensions. At each interface 18 magnesium (Mg) and 18 oxygen (O_s) atoms of the surface were exposed to the water molecules. The system was equilibrated in DFTB molecular dynamics simulations for 100 ps at 298 K. Once the system was equilibrated, we commenced DFTB-based metadynamics to probe the deprotonation and adsorption reactions. Time integration was performed using a 0.5 fs time step in all DFTB-based molecular dynamics simulations.

Protonation and deprotonation states in our simulations were defined based on the O–H distance (r_OH), referred to as a collective variable in metadynamics simulations. The protonated state is the 1st minimum in the free energy profiles, which coincides with the 1st peak position in the O–H pair-correlation function (g(r), Fig. 3a). When the O site in the respective system is deprotonated, it was found that another proton (H⁺) from coordinated water molecules would spontaneously reprotonate the O site. We restricted this reprotonation by constraining all protons, except the one that was being biased, to remain outside of the first co-ordination sphere of the O site. The biased H⁺ ion was allowed to interact with other water molecules outside the first hydration shell of the deprotonated O site and form hydronium ions (H₃O⁺).

2.1.1 Absolute or direct pK_a estimate. For the reaction AH ⇌ A⁻ + H⁺, the acid dissociation equilibrium constant (pK_a) and the Helmholtz free energy difference between protonated and unprotonated states (ΔF) is related as:³⁹

pK_a = βΔF/ln(10); pK_a = −log₁₀ K_a

where β = 1/(k_BT), k_B is the Boltzmann constant, and T is the absolute temperature. In this approach, the energy difference between protonated and unprotonated states is defined using a cutoff radius (r_c). At separation distances between the proton and oxygen that are smaller than r_c, the moiety is considered protonated. At distances greater than r_c, it is assumed that the covalent bond between the H⁺ and O has been broken. The donor (OH⁻) and acceptor (H₃O⁺) pair are at finite separation in the simulation, and equivalent pK_a units are added to make infinite separation in our calculations.⁴⁰

2.1.2 Probabilistic pK_a estimate. Another approach to calculating pK_a is by normalizing the probability of finding a proton within a cutoff distance (r_c) by a probability of maximum possible distance (r_max).⁴¹ The equilibrium constant (K_a(r_c)) for weak acids is
where N is the number of site (=1), V is the volume of the simulation system, and c_o is the standard concentration (=1 M). The cut-off distance r_c is computed as that which gives the correct pK_w = 14 for bulk-like water using K_w(r_c) = ((1 − α(r_c))²·(N_w/c₀V)). For water, this is r_c = 1.13 Å, as determined for the 3OB water model using DFTB+.⁴²

2.2 Materials

The MgO nanometer-sized powder used in our experiments was obtained from Beantown Chemicals (BTC) (40–60 nm APS, 99 wt%). These were pre-treated by heating under vacuum at 400 °C for 24 h to eliminate the adsorbed water, surface hydroxyls and carbonate phases formed by reaction with ambient air. The surface area of the pre-treated powder was analyzed using a Flex3 BET instrument (Micromeritics, Norcross, GA) and was 32.20 ± 0.24 m² g⁻¹. The MgO single crystal used here was synthesized using a carbon arc-fusion technique at ORNL.⁴³ Immediately prior to the experiment, a sample of this material was cleaved using a razor blade roughly along the (100) direction to provide fresh surfaces. The final sample was cubic, roughly 2–3 mm in length.

2.3 Hydroxylation experiments

Hydroxylation experiments were conducted on two different MgO samples (MgO powder and the freshly-cleaved single crystal). The pH measurements in both experiments were conducted at 27 °C using a Denver instruments pH/ISE Conductivity 250 model. For powder experiments, 0.1 gram of powder was suspended in 25 mL of deionized water. We measured solution pH at each second while stirring using a magnetic stir bar. The crystals (∼0.2 g) were placed in 10 mL of deionized water and reacted in a time series for 5 to 64 days.

3 Results and discussion

To understand the MgO(100)–water interfacial structure, we calculated one-dimensional (1d) atomic density profiles normal to the surface (Fig. 1a), two-dimensional (2d) atomic density distribution of exposed surface atoms and the directly associated atoms of water molecules (Fig. 1b and c). The pair co-relations and respective co-ordination numbers are also calculated to probe the hydration of exposed surface atoms (Fig. 2b), and in the aqueous system as well (Fig. 2a).

Fig. 1 Structure of the MgO(100)–water interface. (a) Atomic density profiles of slab magnesium (Mg–green) and oxygen (O_s–red) atoms, and the water molecule oxygen (O_w–blue) and hydrogen (H–gray) atoms. Inset shows a snapshot of the interface after equilibration. At the interface, the hump at ∼0.5 Å in the O_s peak (red line), which projects towards the solution, is due to hydroxylated O_s atoms by the H⁺ at ∼1.45 Å. (b) Two-dimensional (2d) atomic density distribution isosurfaces (small region) showing H⁺ ions hydroxylating the exposed surface O_s atoms. As highlighted in the dotted rectangle, distinctive pairs of neighboring surface hydroxyls are observed at both interfaces in unbiased simulations. (c) 2d atomic density distribution of the water molecules in the first layer (small region). These water molecules hydrate the exposed surface Mg–atoms, as shown in the dotted circles. A few of the exposed Mg–atoms are in a non-hydrated state as shown in the dotted triangle. The structure of neighbouring terminal water molecules made them inaccessible for hydrating the few exposed Mg–atoms on the surface. For clarity, when highlighting surface hydroxyls in (b) terminal water molecules were not displayed; similarly in (c), surface hydroxyls are not shown.

Fig. 2 Pair correlations in aqueous and MgO–water systems. (a) Pair correlation (g(r)) and coordination number (Cn) of the aqueous magnesium ion (Mg) with the oxygen (O_w) and hydrogen (Hw) of water molecules, O_w with Hw of water molecules in the aqueous magnesium system. (b) g(r) and Cn of surface atom Mgs with other Mgs, and Mgs with O_w of terminal water molecules at the MgO–water interface.

3.1 Interfacial structure

The 1d atomic density profiles of the slab atoms (Mg and O_s) and the water molecules (O_w and H) are presented in Fig. 1a. The density peaks of the slab atoms are broadened due to atomic vibrations around their mean positions, but the slab atoms remained in their respective layers. These layers are separated by ∼2.2 Å, which is close to the MgO(100) d-spacing of 2.1 Å of Hazen.³⁸ The discrepancy is likely due to the finite thickness of the slab. Ordered water molecules at the MgO(100)–water interface form distinct high density O_w and H peaks. The atoms in the 1st O_w (2.4 Å) and H peaks (1.45 Å), and the 2nd H peak (3.0 Å) are directly associated with exposed surface atoms (Mg and O_s), as shown schematically in Fig. 1a. Some of the surface-bound water molecules dissociated into H⁺ and OH⁻ ions during unbiased 30 ps simulations. The H⁺ ions hydroxylated part of the exposed surface by adsorbing onto bridging oxygen (O_s) atoms (the distance between them is ∼0.98 Å). These are the μ₅-hydroxo moieties mentioned earlier, appearing as the 1st H peak at ∼1.45 Å from the surface (Fig. 1a). The hydroxylated O_s were displaced from their bulk-like positions, moving away from the bulk and towards the water (0.47 Å). The displaced oxygen atoms appeared as an asymmetry in the exposed layer O_s peak in the surface atomic density profile. They formed a pair of hydroxyls on the surface as shown in the two-dimensional density distribution isosurfaces (Fig. 1b). In total, ∼25% (9 out of 36 O_s) of the exposed surface (O_s) atoms became hydroxylated over the course of the unbiased simulation. To test if the magnitude of surface hydroxylation is an artifact of the simulation cell size, we equilibrated a larger system of 16.84 × 16.84 × 37.58 ų (256 slab atoms and 295 water molecules). The ∼25% fraction of hydroxylated O_s is consistent (16 out of 64 O_s) between both systems, suggesting that system size effects are minimal. This fraction is also in agreement with a prior DFT-based study using a smaller cell.²³ The 1st O_w peak at ∼2.4 Å and the 2nd H peak at ∼3.0 Å (Fig. 1a) correspond to water molecules hydrating Mg surface sites at the interface (shown in Fig. 1c). We analyzed the probable distance of these water molecules from exposed surface Mg atoms (Mgs) as well as the distance between Mgs atoms on the surface by calculating the pair correlations. We presented these pair correlations Mgs–Mgs and Mgs–O_w and respective coordination numbers in Fig. 2b. These water molecules are referred to as terminal waters (η-OH₂) and their distance is 2.16 Å from the surface. The location of the η-OH₂ moieties (Mgs–O_w in Fig. 2b and S1†) is consistent with the reported value of 2.15 Å by Ding and Selloni.²² The high density peaks further from the interface (2nd of O_w and 3rd of H) are due to structuring of the water molecules by hydrogen bond formation with η-OH₂ molecules. The OH⁻ ions from the dissociated water molecules were also present in this layer, consistent with prior DFT observations.²³

Not all exposed Mg atoms have an η-OH₂ molecule coordinated to them (dotted triangle in Fig. 1c). This is reflected in a reduction of the average Mg coordination of 0.7 rather than 1.0 that one might expect (Cn_Mgs–Ow in Fig. 2b). 30% of the exposed Mg in non-hydrated state is likely due to steric constraints on the surface. Specifically, the distance between nearest coordinated Mg atoms is 3.0 Å (g(r)_Mgs–Mgs in Fig. 2b) and this distance is too small to accommodate terminal waters on each site with full rotational freedom (∼3.36 Å for the 1st shell, see g(r)_O–O in Fig. S5†). In addition, the presence of hydroxylated bridging oxygens of surface sites in the vicinity (Fig. 1b and c), which are closer to the surface layer than η-OH₂, limits the orientation probability and thus forces the η-OH₂ molecules into an ordered structure. Both the steric constraints and limited orientation probability caused the few non-hydrated states of Mg–atoms at the interface.

What is most significant about the interfacial atomic density analysis above is the observed spontaneous dissociation of adsorbed water molecules to form MgOH⁺_(surf) and OH⁻_aq at T = 298.15 K. This is the first step in the hydroxylation reaction described in the introduction (eqn (1)).^14–17 To determine the thermodynamics of the subsequent surface hydroxylation step, we performed free-energy calculations and determined the pK_a of the site specific reactions mentioned in eqn (5)–(7).

3.2 Free energy and pK_a

To validate the method to calculate the energetics of the proton transfer reactions on surface sites, we computed the pK_a for the aqueous Mg²⁺ ion and compared the result with experimentally reported value of pK_a = 11.4.⁴⁴ The specific reaction is: Mg(H₂O)²⁺_6aq → Mg(H₂O)₅OH⁺_aq + H⁺.

For the aqueous magnesium ion, the pair correlations of Mg–O_w, Mg–O_w and O_w–H_w in the aqueous system are presented in Fig. 2a. The first maximum in g(r) for Mg–O_w is observed at 2.09 Å (Fig. 2a), in agreement with the value of 2.09 Å derived from X-ray diffraction,⁴⁵ and close to those observed in previous simulations at 2.104 Å,⁴⁶ 2.13 Å (ref. 47) and 2.11 Å (ref. 48) from DFT, and 1.96 Å from classical MD.^37,49

During the deprotonation, a distance constraint was maintained between Mg²⁺ and O_w of the hydrating water molecules at r_Mg–Ow ≤ 2.5 Å to prevent exchange of water molecules during the deprotonation. This was necessary because the deprotonation reaction is endothermic and reprotonation would occur spontaneously. The other H_w atom of the sampling water molecule was also restrained by using an upper constraint at r_O–H = 1.0 Å to avoid hopping and occupying the deprotonating site. During these biased simulations the first hydration shell water molecules had rotational degrees of freedom. In the biased simulations, hopping of H atoms between water molecules was allowed, except with the sampling O site of the water molecule. The magnitudes of all the distance constraints were selected based on the aqueous g(r) profiles (Fig. 2a).

The converged free-energy profile for deprotonation of the aqueous Mg²⁺ ion is presented in Fig. 3a (full profiles in Fig. S2†). The time evolution of the free energies is shown in the ESI (Fig. S3†). The free energy profile has a global minimum at the unbiased O–H distance, r = 0.97 Å. This distance corresponds to the covalent bond length between O and H atoms of the water molecules. A steep increase in the free energy of the aqueous Mg²⁺ ion is observed at r ≤ 2.5 Å (Fig. 3a and S2†). At separations beyond r > 2.5 Å, the free energy profile flattened. As shown in the g(r) of O_w–H_w, H is free from O coordination at r_OH = 1.16 Å (Fig. 2a). This means the covalent bond between O and H is broken, but a strong electrostatic interaction still exists between them at this distance. The H⁺ ion is associated with the O of water molecules in a further hydration shell to form a hydronium (H₃O⁺) ion when r_OH ≥ 1.5 Å (Fig. 2a). Beyond the second hydration shell (r > 2.5 Å), the probability, and the energetic contribution for forming H₃O+, is similar to that of the initial dissociation, which is reflected in the flattening of the free energy profiles at large distances.

Fig. 3 Deprotonation free energy and probabilistic pK_a. (a) Free energy profiles for deprotonation of the aqueous Mg²⁺ ion, the hydroxylated bridging oxygen and the terminal water of exposed Mg atoms. Inset shows the deprotonation free energy profiles up to 5 Å. (b) Estimated probabilistic pK_a as a function of the distance at which deprotonation is considered to have occurred. The intercept of the vertical dotted line with the hydroxylated bridging oxygen and terminal water pK_a profiles allows the estimation of the pK_a values as 13.5 and 5.5 respectively.

In most ion pairing and ion adsorption studies, the reactant and product states have energy minima separated by well-defined transition states.^50–53 This contrasts with the behavior in the deprotonation free energy profile observed here, where there is no distinguishable energy minimum for the deprotonated state and the transition state is not well defined. As discussed above, the H⁺ ion leaves first hydration at 1.16 Å and associates with O of higher hydration shells at r_OH ≥ 1.5 Å (Fig. 2a). We therefore defined the midpoint of the 1st and 2nd hydration shells r_c = 1.33 Å as the cutoff distance for estimating the absolute pK_a. However, due to the system size limitation in simulations, the deprotonated O site (OH⁻) and probable hydronium forming O site (H₃O⁺) are not truly infinitely separated. We roughly estimated the energy required to separate this pair (OH⁻–H₃O⁺) to infinity from 2.73 Å (1st maximum in g(r)_OO, Fig. S5†), which is the equivalent of 1.17 pK_a units (detailed calculation are given in ESI†), which is added to the absolute pK_a values.⁴⁰

Thomsen and Shiga defined the probabilistic cutoff distance as r_c = 1.13 Å to achieve a bulk water pK_a of 14.0 using the same 3OB water model and the DFTB method.⁴² In our simulations water molecules are represented by a further refined (O–H bonding) water model 3OBw which produces bulk water density. At the cutoff r_c = 1.13 Å, we observed that the probabilistic pK_a of the aqueous Mg²⁺ ion is 12.3 (Fig. 3b), which is marginally higher than the experimentally reported value of 11.4.⁴⁴ The choice of the transition state distance may also have affected this result in that pK_a = 11.4 is obtained if r_c = 1.12 Å is used instead of 1.13 Å. Since all our investigations of deprotonation involve surface sites directly bonded to Mg, we defined the cutoff distance for subsequent probabilistic approach calculations as r_c = 1.12 Å. This choice may result in fortuitous cancellation of errors if both the energies and the structure in the DFTB calculations are incorrect, but we view this as the best estimate of the transition state distance that can be obtained practically. The computed pK_a values using these defined cutoffs are presented in Table 1.

Table 1 pK_a estimates of the aqueous Mg²⁺ ion, hydroxylated bridging oxygen (μ₅-hydroxo) and terminal water (η-OH₂). In the probabilistic approach the cutoff (r_c) value is derived from that required to achieve an experimental aqueous Mg²⁺ ion pK_a of 11.4^a (ref. 44)

Deprotonating species	Probabilistic pKa (rc = 1.12 Å)	Absolute pKa (rc = 1.33 Å)
a pKa values at rc = 1.13 Å, which is based on bulk water pKa of 14.0, are presented in brackets.
Aqueous Mg^2+	11.4 (12.3)	14.9
μ5-hydroxo	13.5 (13.9)	13.8
η-OH2	5.5	5.6

---

The estimated aqueous Mg²⁺ ion pK_a value from the absolute method (14.9) is higher than the experimental estimate (11.3). This difference is primarily due to the steep rise of free energy with r_c and the difference in the r_c values selected (r_c = 1.13 Å vs. r_c = 1.33 Å)(Fig. 3b and Table 1). Regardless of these small differences, however, this free energy dependent pK_a calculation proves that the parameters selected in the simulations are useful for probing the deprotonation reaction.

We extended the metadynamics simulations to probe the MgO surface site-specific hydroxylation reactions mentioned earlier (eqn (5)–(7)). To restrict proton hopping from neighbouring water molecules, constraints similar to aqueous Mg²⁺ were also used while deprotonating the η-OH₂ of the MgO surface to form 〉Mg–OH⁺ (〉MgOH²⁺₂ + OH⁻ ⇌ 〉MgOH⁺ + H₂O) and hydroxylating the non-hydrated surface Mg–atoms (〉Mg²⁺ + OH⁻ ⇌ 〉MgOH⁺). During the deprotonation of the μ₅-hydroxo (〉OH⁻ + OH⁻ ⇌ 〉O²⁻ + H₂O), its surface O and, when deprotonating the η-OH₂, the surface Mg atom, are constrained with sub-surface atoms to prevent them from dissolving into the water due to the applied bias. The deprotonation free energy profiles and respective probabilistic pK_a of these η-OH₂ and μ₅-hydroxo sites are shown in Fig. 3a and b, and their absolute and probabilistic pK_a values are also presented in Table 1. A negligible difference between the estimated probabilistic and absolute pK_as is observed for both η-OH₂ and μ₅-hydroxo sites. The deprotonation free energy profiles have a concavity and close to flattening behaviour at 1.13 Å ≤ r_c ≤ 1.33 Å resulting in negligible pK_a differences. The pK_a of the μ₅-hydroxo site (13.5) is much higher than that of the η-OH₂ site (5.5) which means bridging oxygen hydroxyls are much more stable than terminal water molecules. These pK_a values (Table 1) reveal to us the possibility of exchanging H⁺ from the μ₅-hydroxo site and η-OH₂ site at the MgO(100)–water interface. In the absence of another pH buffers in the macroscopic solution, we would expect that the MgO(100) surface would contain hydroxylated η-OH₂ and μ₅-hydroxo sites.

To test if the calculated pK_a values for the surface sites are plausible, we experimentally measured the pH of solutions in contact with MgO as a function of time. In otherwise unbuffered solutions of a MgO nanoparticulate powder, or in contact with single crystal MgO (Fig. 4a and b, respectively), the initial pH was 7.5 and 7.1 respectively. The pH rose over time to a value of 10.5 in nanoscale powder hydration and 10.4 in single crystal hydration. To relate the calculated surface acidity constants (pK_a) to these experimentally-observed pHs, we simulated the MgO–water interface and estimated the resultant pH using the aqueous speciation model PHREEQC⁵⁴ using the included WATEQ4F thermodynamic database for aqueous solutes. PHREEQC relies on the Davies equation for the relationship between activity and concentration and the Gouy–Chapman diffuse layer for balancing the charge from a surface. The WATEQ4F database solubility products were relied upon for aqueous speciation (e.g., carbonate species, aqueous magnesium hydroxylation), but solubility products for the magnesium-containing solid phases were defined in the input files based on literature data since these were not available in the database. The calculated deprotonation constants (pK_a) using the probabilistic method were chosen, 13.5 for μ₅-hydroxo and 5.5 for μ-OH₂ deprotonation as site-specific reaction input parameters. The BET surface area of the MgO powder (32.20 ± 0.24 m² g⁻¹), solubility data^55,56 and estimated active sites on the surface based on the (100) plane were also taken as surface parameters. The input script and result from the PHREEQC are given in the ESI.†

Fig. 4 Variation of the solution pH in MgO hydration experiments. (a) Change in pH as a function of time during the hydration of the (a) BTC MgO powder and (b) MgO single crystal. Time resolution of pH measurements in the powder hydration experiments is 1 second.

We tested two scenarios using PHREEQC. In the first case we inhibited MgO dissolution but allowed surface sites to equilibrate with the solution to mimic our single crystal (SC) and MgO nanoscale powder hydration experiments. The model resulted in a solution pH of 7.1 which agrees with the observed starting pHs of 7.5 for the nanoparticles (Fig. 4a) and 7.1 for single crystals (Fig. 4b). In the second case, we allowed MgO dissolution to bring the system to equilibrium. This resulted in a calculated pH of 12.0 which is somewhat comparable with our observed final pH value of 10.5 for the nanoparticle and 10.4 for the single crystal experiments (Fig. 4). This similarity between the pH values from the fundamental aqueous model and those observed experimentally suggests that the calculated pK_a values from the DFTB simulations are indeed plausible, and that the initial pH of a suspension of MgO and water is buffered by the surface sites, whereas the final pH of the solution is driven more towards that governed by dissolution of MgO.

Beyond the surface pK_as, the mechanism in the introduction includes adsorption of hydroxide to vacant Mg surface sites (eqn (2)). To verify the energetics of formation of MgOH_(surf) +·OH⁻ in the second step of the mechanism, we performed a biased simulation to drive aqueous OH⁻ to adsorb onto an uncoordinated surface Mg atom. The free energy profile is shown in Fig. 5, with that of deprotonation of the η-OH₂ for comparison. The energy barrier for the OH⁻ ion to adsorb at the MgO surface is significant (∼52 kJ mol⁻¹) and larger than the barriers to deprotonate the terminal water of a magnesium surface site. The global minimum is at the separation distance of r = 6.5 Å which is far away from the hydration shell of Mg. This indicates that OH⁻ prefers to stay in the solution instead of adsorbing onto the MgO surface. The free-energy profile of terminal water deprotonation shows a local minimum at r = 1.55 Å (Fig. 5 and S2†), and the energy barrier between the global (protonated state) and local minimum (deprotonated state) is ∼26 kJ mol⁻¹. This is ∼26 kJ mol⁻¹ less than the energy required to bring an OH⁻ from the bulk fluid region to form an MgOH⁺_(surf) ·OH⁻ by the OH_(aq)– molecule on the surface. The probabilistic pK_a of the terminal water molecule is 5.4 (absolute pK_a = 5.5), which makes deprotonation of the terminal water molecule highly favorable at our experimentally observed pH > 10.4 when MgO dissolution buffers pH. In contrast, the proposed adsorption step in the reaction mechanism (eqn (2)), aqueous OH⁻ forming MgOH⁺_(surf) ·OH⁻,^14,16 is energetically unfavorable and therefore less likely. This finding, of favourable terminal water deprotonation, is inconsistent with a prior interpretation of water dissociation and hydroxide adsorption on surface Mg sites.⁵⁷ In addition, on the surface the μ₅-hydroxo (〉OH⁺) with pK_a of 13.5 is also stable at the experimentally observed pH > 10.4 (Table 1). This indicates these μ-hydroxo sites are present in the vicinity of deprotonating η-OH₂ sites and facilitate forming MgOH_(surf) +·OH⁺ on the surface. Thus, both our computational and experimental results suggest replacement of the adsorption step (eqn (2)) in the mechanism by:

Deprotonation: 〉MgOH²⁺₂ + OH⁻ ⇌ 〉MgOH⁺ + H₂O (8)

Fig. 5 Free energy profile of OH⁻ adsorption on the MgO surface to form MgOH⁺_(surf)⋅OH⁻. The presence of a global minimum at r = 6.4 Å shows that the adsorption of aqueous OH⁻ on a surface Mg is not energetically favorable, as a significant free energy barrier (∼52 kJ mol⁻¹) has to be overcome for OH⁻ adsorption. The terminal water (η-H₂O) deprotonation free-energy profile has the local minimum at r = 1.55 Å at an energy barrier between the local minimum and hydration region of ∼26 kJ mol⁻¹.

4 Summary and conclusion

Density functional tight binding atomic-scale simulations, coupled with the metadynamics rare event method, were used to analyze the surface chemistry of MgO and the protonation of its surface sites. The computed probabilistic pK_a = 12.3 for the aqueous Mg²⁺ ion with the reference of bulk water pK_a of 14.0 is in fair agreement with a reported experimental value of 11.4, with the difference between them likely due to a small error in the choice of transition state distance (0.01 Å). We defined the cut off r_c = 1.12 Å for probabilistic estimation on the basis of the experimental aqueous Mg²⁺ ion. The pK_a for hydroxylation of a terminal water bound to magnesium surface sites was found to be 5.5 and that of deprotonation for a hydroxyl formed from a bridging oxygen was 13.5. Support for the plausibility of these values was derived by using these pK_as in an aqueous speciation/surface model in PHREEQC and comparing the solution pH to experiments on MgO nanoparticles and a single crystal. The model was consistent with initial solution pHs of 7.1 (vs. experimental pH = 7.4 and 7.1, respectively) when surface sites buffered solution pH, suggesting these were plausible pK_as for the surface sites.

As discussed in the manuscript, the first reaction in the process of conversion of MgO to Mg(OH)₂ is the hydroxylation of surface bridging oxygens, which is observed even in unbiased simulations and is supported by the above pK_a values. The subsequent reaction is hydroxide adsorption and has two possible mechanisms, absorption of aqueous OH⁻ ions or terminal water deprotonation. We found the free energy barrier for the adsorption of aqueous OH⁻ ions (∼52 kJ mol⁻¹) is significantly greater than that for the terminal water deprotonation (∼26 kJ mol⁻¹). The latter also has a pK_a of ∼5.4, which makes the terminal deprotonation reaction highly favorable at the observed experimental pH (10.4). Thus, the deprotonated terminal water bound to a magnesium surface site, 〉MgOH⁺_(surf), most likely forms MgOH_(surf) +·OH⁻ when the terminal water loses a proton and not through hydroxide ion adsorption. Our pK_a calculations showed that surface bridging oxygen hydroxyls (〉OH⁺_(surf)) are also stable at the alkaline pHs typically found when Mg(OH)₂ forms, which is also where MgO dissolution buffers the solution pH. These results suggest that to precipitate Mg(OH)₂, a Mg²⁺ion must dissolve, potentially from a hydroxylated magnesium site.

Data availability

The unit cell source for building the MgO surface is cited and given in the reference section. The source of publicly available DFTB software is also cited and given in the reference section. The input file for running the aqueous-speciation PHREEQC model and output data are given in the ESI.† Information provided in the manuscript and ESI† is sufficient to reproduce the results of PHREEQC calculations. All required files for running DFTB+ simulations with metadynamics of (i) aqueous simulations, (ii) MgO(100)–water simulations, (iii) SK parameter files representing systems are available at https://code.ornl.gov/sr4/mineral_water_dftb_metadynamics.git. For better understanding of each command line and structure of code, we will upload all these file with DOI number in the PLUMED-NEST (a public repository of the PLUMED consortium).

Author contributions

Sai Adapa: conceptualization, data curation, formal analysis, investigation, methodology, software, validation, visualization, writing – original draft. Ke Yuan: conceptualization, methodology, resources, software, writing – review & editing. Barbara R. Evans: investigation, writing – review & editing. Juliane Weber: funding acquisition, investigation, writing – review & editing. Stephan Irle: software, writing – review & editing. Lawrence M. Anovitz: writing – review & editing. Andrew G. Stack: conceptualization, data curation, formal analysis, investigation, methodology, resources, software, supervision, validation, visualization, writing – review & editing.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was primarily supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences (BES), Materials Sciences and Engineering Division. S.I. provided advice and correction on DFTB and the 3OBw implementation using support from the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, Chemical Sciences, Geosciences, and Biosciences Division. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231. This research used resources of the Compute and Data Environment for Science (CADES) at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725.

References

1. J. Wang, L. Huang, R. Yang, Z. Zhang, J. Wu, Y. Gao, Q. Wang, D. O'Hare and Z. Zhong, Energy Environ. Sci., 2014, 7, 3478–3518.
2. M. E. Boot-Handford, J. C. Abanades, E. J. Anthony, M. J. Blunt, S. Brandani, N. Mac Dowell, J. R. Fernández, M.-C. Ferrari, R. Gross and J. P. Hallett, et al., Energy Environ. Sci., 2014, 7, 130–189.
3. M. Hanifa, R. Agarwal, U. Sharma, P. C. Thapliyal and L. P. Singh, J. CO2 Util., 2023, 67, 102292.
4. L. A. Hollingbery and T. R. Hull, Thermochim. Acta, 2012, 528, 45–52.
5. N. T. Dung, R. Hay, A. Lesimple, K. Celik and C. Unluer, Cem. Concr. Compos., 2021, 115, 103826.
6. V. A. Kuuskraa, M. L. Godec and P. Dipietro, Energy Procedia, 2013, 37, 6854–6866.
7. K. Jessen, A. R. Kovscek and F. M. Orr, Energy Convers. Manage., 2005, 46, 293–311.
8. E. J. Beckman, Ind. Eng. Chem. Res., 2003, 42, 1598–1602.
9. A. Sinha, L. A. Darunte, C. W. Jones, M. J. Realff and Y. Kawajiri, Ind. Eng. Chem. Res., 2017, 56, 750–764.
10. D. W. Keith, G. Holmes, D. S. Angelo and K. Heidel, Joule, 2018, 2, 1573–1594.
11. N. McQueen, K. V. Gomes, C. McCormick, K. Blumanthal, M. Pisciotta and J. Wilcox, Prog. Energy, 2021, 3, 032001.
12. N. McQueen, P. Kelemen, G. Dipple, P. Renforth and J. Wilcox, Nat. Commun., 2020, 11, 3299.
13. J. Weber, V. Starchenko, K. Yuan, L. M. Anovitz, A. V. Ievlev, R. R. Unocic, A. Y. Borisevich, M. G. Boebinger and A. G. Stack, Environ. Sci. Technol., 2023, 57, 14929–14937.
14. L. F. Amaral, I. R. Oliveira, R. Salomão, E. Frollini and V. C. Pandolfelli, Ceram. Int., 2010, 36, 1047–1054.
15. Z. Xing, L. Bai, Y. Ma, D. Wang and M. Li, Materials, 2018, 11, 1835.
16. J. M. Rimsza, E. G. Sorte and T. M. Alam, ACS Omega, 2019, 4, 1033–1044.
17. M. Pettauer, A. Baldermann, S. Eder and M. Dietzel, Cryst. Growth Des., 2024, 24, 3085–3092.
18. D. O. Scanlon, A. Walsh, B. J. Morgan, M. Nolan, J. Fearon and G. W. Watson, J. Phys. Chem. C, 2007, 111, 7971–7979.
19. E. Carrasco, M. A. Brown, M. Sterrer, H. J. Freund, K. Kwapien, M. Sierka and J. Sauer, J. Phys. Chem. C, 2010, 114, 18207–18214.
20. E. Carrasco, A. Aumer, J. F. Gomes, Y. Fujimori and M. Sterrer, Chem. Commun., 2013, 49, 4355–4357.
21. N. M. Adhikari, A. Tuladhar, Z. Wang, J. J. D. Yoreo and K. M. Rosso, J. Phys. Chem. C, 2021, 125, 26132–26138.
22. Z. Ding and A. Selloni, J. Chem. Phys., 2021, 154, 114708.
23. M. Sassi and K. M. Rosso, Phys. Chem. Chem. Phys., 2023, 26, 2269–2276.
24. M. Hollerer, D. Prochinig, P. Puschnig, E. Carrasco, H. J. Freund and M. Sterrer, J. Phys. Chem. C, 2019, 123, 3711–3718.
25. R. Włodarczyk, M. Sierka, K. Kwapień, J. Sauer, E. Carrasco, A. Aumer, J. F. Gomes, M. Sterrer and H. J. Freund, J. Phys. Chem. C, 2011, 115, 6764–6774.
26. B. Hourahine, B. Aradi, V. Blum, F. Bonafé, A. Buccheri, C. Camacho, C. Cevallos, M. Y. Deshaye, T. Dumitric, A. Dominguez, S. Ehlert, M. Elstner, T. V. D. Heide, J. Hermann, S. Irle, J. J. Kranz, C. Köhler, T. Kowalczyk, T. Kubař, I. S. Lee, V. Lutsker, R. J. Maurer, S. K. Min, I. Mitchell, C. Negre, T. A. Niehaus, A. M. Niklasson, A. J. Page, A. Pecchia, G. Penazzi, M. P. Persson, J. Řezáč, C. G. Sánchez, M. Sternberg, M. Stöhr, F. Stuckenberg, A. Tkatchenko, V. W. Yu and T. Frauenheim, J. Chem. Phys., 2020, 152, 124101.
27. K. H. Lee, U. Schnupf, B. G. Sumpter and S. Irle, ACS Omega, 2018, 3, 16899–16915.
28. X. Lu, M. Gaus, M. Elstner and Q. Cui, J. Phys. Chem. B, 2015, 119, 1062–1082.
29. P. Goyal, H. J. Qian, S. Irle, X. Lu, D. Roston, T. Mori, M. Elstner and Q. Cui, J. Phys. Chem. B, 2014, 118, 11007–11027.
30. T. Krüger, M. Elstner, P. Schiffels and T. Frauenheim, J. Chem. Phys., 2005, 122, 114110.
31. H. A. Witek, S. Irle, G. Zheng, W. A. de Jong and K. Morokuma, J. Chem. Phys., 2006, 125, 214706.
32. M. Gaus, Q. Cui and M. Elstner, J. Chem. Theory Comput., 2011, 7, 931–948.
33. J. Ono, M. Imai, Y. Nishimura and H. Nakai, J. Phys. Chem. B, 2020, 124, 8524–8539.
34. H. Nakai, A. W. Sakti and Y. Nishimura, J. Phys. Chem. B, 2016, 120, 217–221.
35. M. Bonomi, D. Branduardi, G. Bussi, C. Camilloni, D. Provasi, P. Raiteri, D. Donadio, F. Marinelli, F. Pietrucci, R. A. Broglia and M. Parrinello, Comput. Phys. Commun., 2009, 180, 1961–1972.
36. G. A. Tribello, M. Bonomi, D. Branduardi, C. Camilloni and G. Bussi, Comput. Phys. Commun., 2014, 185, 604–613.
37. S. Adapa and A. Malani, Sci. Rep., 2018, 8, 12198.
38. R. M. Hazen, Am. Mineral., 1976, 61, 266–271.
39. J. E. Davies, N. L. Doltsinis, A. J. Kirby, C. D. Roussev and M. Sprik, J. Am. Chem. Soc., 2002, 124, 6594–6599.
40. R. Wang, V. Carnevale, M. L. Klein and E. Borguet, J. Phys. Chem. Lett., 2020, 11, 54–59.
41. M. Schilling and S. Luber, Inorganics, 2019, 7, 7060073.
42. B. Thomsen and M. Shiga, J. Chem. Phys., 2021, 154, 084117.
43. M. M. Abraham, C. T. Butler and Y. Chen, J. Chem. Phys., 1971, 55, 3752–3756.
44. C. Baes and R. Mesmer, The Hydrolysis of Cations, Kriegeer, Malabar, FL, 1986.
45. Y. Marcus, Chem. Rev., 1988, 88, 1475–1498.
46. L. Bernasconi, E. J. Baerends and M. Sprik, J. Phys. Chem. B, 2006, 110, 11444–11453.
47. F. C. Lightstone, E. Schwegler, R. Q. Hood, F. Gygi and G. Galli, Chem. Phys. Lett., 2001, 343, 549–555.
48. K. Yuan, N. Rampal, S. Adapa, B. R. Evans, J. N. Bracco, M. G. Boebinger, A. G. Stack and J. Weber, ACS Appl. Mater. Interfaces, 2024, 16, 64233–64243.
49. S. Mamatkulov, M. Fyta and R. R. Netz, J. Chem. Phys., 2013, 138, 024505.
50. J. R. Rustad and A. G. Stack, J. Am. Chem. Soc., 2006, 128, 14778–14779.
51. A. G. Stack, P. Raiteri and J. D. Gale, J. Am. Chem. Soc., 2012, 134, 11–14.
52. M. N. Joswiak, B. Peters and M. F. Doherty, Cryst. Growth Des., 2018, 18, 6302–6306.
53. A. Silvestri, P. Raiteri and J. D. Gale, J. Chem. Theory Comput., 2022, 18, 5901–5919.
54. D. L. Parkhurst, C. A. J. Appelo and U. S. G. Survey, Description of Input and Examples for PHREEQC Version 3: a Computer Program for Speciation, Batch-Reaction, One-Dimensional Transport, and Inverse Geochemical Calculations, 2013.
55. W. Hummel, U. Berner, E. Curti, F. J. Pearson and T. Thoenen, Radiochim. Acta, 2002, 90, 805–813.
56. R. Hill, J. Canterford and F. Moyle, Mineral. Mag., 1982, 46, 453–457.
57. P. Yang, J. N. Bracco, G. Camacho Meneses, K. Yuan, J. E. Stubbs, M. D. Boamah, M. Brahlek, M. Sassi, P. J. Eng and M. G. Boebinger, et al., Environ. Sci. Technol., 2025, 59(7), 3484–3494.

---

Footnotes

† Note to the Publisher: The U.S. government retains and the publisher, by accepting the article for publication, acknowledges that the U.S. government retains a non-exclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this manuscript or allow others to do so for U.S. government purposes. The U.S. DOE will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (https://www.energy.gov/downloads/doe-public-access-plan).

‡ Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d5ta03981a

---