Path sampling for atmospheric reactions: formic acid catalysed conversion of SO₃ + H₂O to H₂SO₄

Introduction

SO₃ dispersed in the atmosphere reacts with water, generating H₂SO₄ which is one of the main contributors to acidic rain, as well as contributing to nucleation processes in the upper atmosphere. The SO₃ · H₂O complex can form readily, however early investigations of the hydrolysis of SO₃ to form H₂SO₄ found the uncatalysed rate to be too low to explain the concentration of atmospheric H₂SO₄ (Wang et al., 1988; Hofmann & Von Ragué Schleyer, 1994; Loerting & Liedl, 2000). Further studies revealed that the formation of H₂SO₄ is catalysed by collisions with other atmospheric molecules (Akhmatskaya et al., 1997; Hazra & Sinha, 2011; Long et al., 2012; Torrent-Sucarrat, Francisco & Anglada, 2012; Bandyopadhyay, Kumar & Biswas, 2017; Sarkar, Oram & Bandyopadhyay, 2019), where the common denominator is the facilitation of a double proton transfer via the catalysing molecule (Kumar, Sinha & Francisco, 2016). The relative importance of collisions with different molecules is determined by both the atmospheric concentration of the colliding molecule and the species-dependent reaction rate. For example, while water vapour has a high atmospheric concentration, the reaction barrier remains quite high when an additional water molecule acts as the catalyst (Akhmatskaya et al., 1997; Long et al., 2012). On the other hand, NH₃ is an effective catalyst, lowering the reaction barrier significantly, but its atmospheric concentration is comparatively low (Bandyopadhyay, Kumar & Biswas, 2017; Sarkar, Oram & Bandyopadhyay, 2019).

Recent work has highlighted the importance of formic acid, HCOOH, in catalysing the SO₃ hydration reaction (Hazra & Sinha, 2011; Long et al., 2012). Static quantum chemical calculations of the initial and final states, combined with reaction kinetics and transition state theory (TST), show that the reaction becomes nearly barrier-less, and hence that HCOOH is an effective catalyst. Combined with the recent findings that formic acid concentrations in the atmosphere are higher than previously assumed (Millet et al., 2015; Chaliyakunnel et al., 2016), and the presence of formic acid in volatile organic compounds (VOCs) emitted in boreal forests (Rantala et al., 2015), it seems likely that the formic acid catalysed hydration of SO₃ is a significant contributor to sulphuric acid concentrations in the atmosphere.

Investigations of gas phase SO₃ hydration reactions, however, to our knowledge have heretofore only been performed by static calculations, with kinetics and dynamics only included via TST (Long et al., 2012) or Rice, Ramsperger, Kassel, Markus (RRKM) theory (Hazra & Sinha, 2011). In the present work, we aim to directly investigate the dynamics of the formic acid catalysed conversion of SO₃ · H₂O to H₂SO₄ via path sampling (Van Erp, 2007; Van Erp, Moroni & Bolhuis, 2003; Moroni, Bolhuis & Van Erp, 2004). The method has recently reached a sufficient level of maturity, both theoretically (Riccardi, Dahlen & Van Erp, 2017; Cabriolu et al., 2017) and computationally (Lervik, Riccardi & Van Erp, 2017; Riccardi et al., 2019a). Notable recent work using path sampling has focused on chemical or other transitions between states in aqueous systems, for example water auto-ionization (Moqadam et al., 2018), interactions between biomolecules (Riccardi et al., 2019c) and silicate condensation (Moqadam et al., 2017). The main focus of the present work is to demonstrate the application of path sampling to investigating a new class of gas-phase reactions.

The path sampling approach here chosen is the replica exchange transition interface sampling (RETIS) (Van Erp, 2007; Cabriolu et al., 2017). It permits the investigation of the path space, and thus the determination of a reaction rate, orders of magnitude faster than conventional simulations. This computational enhancement does not come at the cost of resolution since the method provides the same results, when converged, as infinitely long standard simulations. The RETIS methodology will be combined with ab initio molecular dynamics (AI-MD) simulations to predict the reaction rate of the uni-molecular SO₃ hydration reaction with formic acid, considering different density functional theory (DFT) functionals and basis sets, and comparing the results with data available in the literature produced by more approximate methods (e.g., TST).

Materials and Methods

Description of the reaction and DFT methods

We perform DFT-based Born–Oppenheimer molecular dynamics (BO-MD) simulations using the QUICKSTEP module in CP2K (VandeVondele et al., 2005). Before beginning these simulations we considered several different combinations of DFT functionals and basis sets, first investigating the relative energies of the product and reactant states using NWChem (Valiev et al., 2010).

Snapshots of optimised geometries in the reactant (A1 and A2) and product (B1 and B2) states are displayed in Fig. 1. Other computational studies of SO₃ hydration have considered the complete reaction pathway, including the formation of the pre-reactive complex via atmospheric collisions. In this study, we focus only on the uni-molecular reaction involving the double proton transfer from the water molecule via the formic acid to form sulphuric acid, that is $${\mathrm{H}}_2\mathrm{O}\cdot \mathrm{S}{\mathrm{O}}_3\cdot \mathrm{H}\mathrm{C}\mathrm{O}\mathrm{O}\mathrm{H}\stackrel{k_{AB}}{\to }{\mathrm{H}}_2\mathrm{S}{\mathrm{O}}_4\cdot \mathrm{H}\mathrm{C}\mathrm{O}\mathrm{O}\mathrm{H}$$

This reaction involves no collisions, and can be treated as an isomerization reaction, with first-order kinetics described by a single rate coefficient k_AB with units of inverse time. In previous work, this uni-molecular rate coefficient has been computed using TST (Long et al., 2012) and RRKM theory (Hazra & Sinha, 2011). The uni-molecular reaction barrier height was found to be low (0.3 kJ mol ⁻¹ at the MP2/6-311++G(3df, 3pd) level of theory (Hazra & Sinha, 2011)) compared to the barrier for the water catalysed reaction.

Two different local minima in both reactant and product states are present, differing only in the orientation of one OH bond initially belonging to the water molecule. The minimum in the reactant state (A1) differs from that reported in other investigations of this reaction (A2), but is only ∼1 kJ mol ⁻¹ lower in energy. The minimum energy in the product state (B2) is generally several kJ mol ⁻¹ lower in energy than the other local minimum (B1). We remark that we did not consider the final re-organization of the post-reaction H₂SO₄ · HCOOH complex into a lower energy conformation which is the true global minimum (Hazra & Sinha, 2011). The uni-molecular rate coefficient we measure is unaffected by this final rearrangement since it is a spontaneous process. We display some of these results along with comparisons with the literature in Table 1.

Table 1:

Energetics of optimised reactant (A1 and A2) and product (B1 and B2) states using several different combinations of DFT functionals and basis sets.

We also compare with high-level ab initio results already available in the literature (Hazra & Sinha, 2011; Long et al., 2012). We note that the 6-311++G(3df, 3pd) Gaussian basis set is similar to the TZV2P plane-wave basis set. The zero of energy is defined to be the lowest energy reactant state. All energies are expressed in units of kJ mol ⁻¹. Only the result of (Long et al., 2012) includes zero point vibrational energy (ZPE) corrections (in brackets).

Theory/basis set	EA1	EA2	EB1	EB2
BLYP/DZVP	0.0	1.1	8.1	3.5
BLYP/TZVP	0.0	1.2	0.3	−2.9
B3LYP/TZVP	0.0	1.1	−17.4	−22.3
PBE/TZVP	0.0	1.3	−9.2	−13.0
PBE/6-311++G(3df, 3pd)	0.0	1.0	−13.9	−18.4
PBE/aug-cc-pVTZ	0.0	1.2	−12.5	−16.5
PBE0/TZVP	0.0	1.3	−27.3	−32.4
MP2/6-311++G(3df, 3pd) (Hazra & Sinha, 2011)	–	0.0	–	−32.4
CCSD(T)//MP2/aug-cc-pV(T + d)Z (Long et al., 2012)	–	0.0	–	−42.5 (−41.0)

DOI: 10.7717/peerj-pchem.7/table-1

Hybrid functionals such as B3LYP or PBE0 which include some of the Hartree–Fock exchange correlation are necessary to obtain DFT based results which are in good agreement with fully ab initio methods. However, using hybrid functionals in BO-MD simulations would be prohibitively expensive, especially for an exploratory study such as this one. In the end, in our BO-MD simulations, we compare two Generalized Gradient Approximation (GGA) DFT methods, with three different basis sets (Goedecker, Teter & Hutter, 1996). First, we use the BLYP functional (Becke, 1988; Lee, Yang & Parr, 1988) with the DZVP-MOLOPT-GTH basis functions, which we have previously used to investigate formic acid in bulk water (Murdachaew et al., 2016; Hänninen et al., 2018; Daub & Halonen, 2019). Although the BLYP results lack accuracy in terms of the comparative energies of the initial and final states, it was a relatively inexpensive methodology to use for initial investigation. Thereafter, the PBE functional (Perdew, Burke & Ernzerhof, 1996) and both the TZVP-MOLOPT-GTH and TZV2P-MOLOPT-GTH basis sets have been used. There are some differences between the basis sets implemented in NWChem and the plane-wave basis sets employed in CP2K, but respective double-zeta and triple-zeta basis sets should be similar. In particular, we note that the 6-311++G(3df, 3pd) Gaussian basis set produced comparable results to the TZV2P-MOLOPT-GTH basis set. All simulations, both with NWChem and CP2K, used the Grimme D2 dispersion correction (Grimme, 2004, 2006). All BO-MD simulations were performed in the microcanonical (NVE) ensemble in a cubic periodic box of fixed geometry with a box length of L = 8 Å (1 Å = 10⁻¹⁰ m).

Results

As shown in Eq. 2, the computed reaction rate coefficient, k_AB, is derived from the initial flux f_A multiplied by the total crossing probability P_cross. The value of f_A × P_cross at the end of n_cycles of the RETIS simulations is reported in Fig. 2. The value of k_AB obtained via path sampling using the RETIS methodology is directly comparable with the uni-molecular rate coefficient describing the double proton transfer reaction forming H₂SO₄. We show the average result obtained at the end of n_cycles of all RETIS simulations using a given DFT method, basis set and temperature. Our most accurate estimate at T = 300 K has been obtained via PBE calculations with the TZV2P basis set, k_AB = 3.5 × 10⁻⁴ fs⁻¹, with a statistical error of ∼50%. For comparison, a computation of the rate coefficient at T = 300 K based on TST and reaction kinetics calculations (see column 6 of Table S4 in the Supporting Information for Long et al., 2012) gave an estimate of k_AB = 1.9 × 10⁻³ fs⁻¹, while calculations based on the RRKM kinetic theory (see Table 7 in the Supporting Information for Hazra & Sinha, 2011) gave an estimate of k_AB = 4.6 ± 0.1 × 10⁻⁴ fs⁻¹. Overall, our predicted rate coefficient is in reasonable agreement with previous estimates in the literature, despite the fact that we underestimate the energy difference between the reactant and product states (see Table 1). We might reasonably expect that a RETIS calculation using a hybrid functional such as PBE0 with a more accurate energy difference would result in a somewhat larger rate coefficient. Unfortunately, at the present time the computational resources required for such calculations leave the use of hybrid functionals in path sampling calculations beyond our reach.

Our estimated rate coefficients k_AB reported in Fig. 2 include a total of at least 300 ps of simulation data for each RETIS run. In Fig. 3 we show the final value of k_AB (i.e., the value of k_AB at the maximum value of s(r)) as a function of the RETIS cycle number for the PBE calculation using the TZV2P basis set. This plot shows that our results are well-converged, and in fact that reliable results could have been obtained from a simulation using only half the data, or less. More careful optimisation of the interface positions λ_i would also help the results to converge with less total simulation time.

Comparison of results using the same basis set and level of theory (BLYP, DZVP) at two different temperatures (300 K and 250 K) showed that the RETIS simulations predict a smaller rate coefficient at lower temperature. While this may make intuitive sense, it disagrees with the TST results, which predict an increase in the uni-molecular rate coefficient at lower temperature (Long et al. (2012), see column 6 of Table S4 of the Supporting Information). Subsequently, we performed a series of simulations using the most reliable basis set and level of theory available to us (i.e., PBE, TZV2P), at a series of temperatures ranging from T = 250 to 350 K. This approach should allow us to directly compare the temperature dependence observed in the TST calculation with the results of the RETIS method. We summarise these results in Fig. 4.

The rate coefficient seems to be lower at T = 250 K, consistent with the BLYP data shown in Fig. 2, however over the full range of temperature we have sampled we cannot measure any statistically meaningful correlation between k_AB and T. This lack of temperature dependence is in contrast with the monotonic decrease in the uni-molecular reaction rate predicted by TST. It is unlikely that this temperature dependence (or lack thereof) would qualitatively change if the basis set or level of theory were changed.

In the path sampling approach, the temperature dependency of the rate coefficient arises from entropic contributions which are directly computed, and complicating factors such as vibrational anharmonicity are implicitly included. In TST or RRKM, on the other hand, entropic contributions are considered, but in a way that depends on certain assumptions about, for example, the harmonic nature of the vibrational potentials used to construct the vibrational partition function. The rather different temperature dependence of the rate coefficient on we have observed points toward how direct methods like RETIS might lead to different insights into reaction pathways. At the same time, in this case the reaction barrier is low and the uni-molecular reaction rate is very fast compared with the initial collision rate, and therefore the uni-molecular rate coefficient does not play a major role in determining the overall reaction rate (Long et al., 2012).

To investigate the reaction pathway, frame density plots can be used to visualise, in path space, the values of different descriptors. The visualisation and analysis tool of PyRETIS, PyVisA (O. Aarøen, E. Riccardi, 2019, unpublished data), has been used to plot the transition in different path spaces. Correlations, even if partial, can then be detected allowing direct analysis of the transition pathways available from the RETIS simulations. Details about the entire reaction, and not just the static reactant and product states, can be extracted for the whole transition region.

We have examined the three interatomic distances which contribute to our order parameter s(r), as well as the order parameter itself. Those interatomic distances are defined in accordance with the atom labels shown in Fig. 1. In Fig. 5 we display four different correlation plots considering only the accepted paths belonging to the outermost ensemble, i.e. with −4.6 Å < s(r) < −3.9 Å. In the top two figures, the variation in the distance between S and O_W, $r_{{\mathrm{S}\mathrm{O}}_{\mathrm{W}}}$, and the distance $r_{{\mathrm{O}}_1{\mathrm{H}}_1}$ are shown. These are both plotted versus the order parameter s(r). All three of the individual components of s(r) are strongly correlated with the order parameter, as shown by the high values of the correlation variable r².

When looking at the correlation between the individual inter-atomic distances themselves (see bottom 2 plots of Fig. 5), the strongest correlation exists between $r_{{\mathrm{S}\mathrm{O}}_{\mathrm{W}}}$ and $r_{{\mathrm{O}}_1{\mathrm{H}}_1}$. Comparing to Fig. 1, this seems sensible; as the reaction proceeds the distance between the (initial) water oxygen and the sulphur atom varies at the same time as the water hydrogen atom nearest the formic acid begins to transfer over to the formic acid molecule. Correlations with the third component of s(r), that is, $r_{{\mathrm{O}}_2{\mathrm{H}}_2}$, the distance between the acidic hydrogen and the closest SO₃ oxygen, are considerably weaker. According to the bottom right plot of Fig. 5, it seems that $r_{{\mathrm{O}}_2{\mathrm{H}}_2}$ is almost constant until $r_{{\mathrm{O}}_1{\mathrm{H}}_1}$ < 1.1 Å, at which point $r_{{\mathrm{O}}_2{\mathrm{H}}_2}$ suddenly changes from ∼1.4 Å down to 1.2 Å. In other words, we suggest that the final step in the double proton transfer reaction is the change in $r_{{\mathrm{O}}_2{\mathrm{H}}_2}$, and it does not begin without the other parts of the reaction having occurred. These findings are consistent with, and further explain the limited data about the transition region available from a TST investigation, which showed that r_O1H1 ∼ 1.2 Å at the transition state, while $r_{{\mathrm{O}}_2{\mathrm{H}}_2}$ is larger, ∼1.4 Å at the transition state (Long et al., 2012).

Conclusions

We have made use of path sampling BO-MD simulations to study the hydrolysis of SO₃ when catalysed by formic acid. We discuss in detail the reaction pathway(s) that can describe the conversion from SO₃ + H₂O to H₂SO₄. The mechanism has been highlighted and the overall rate computed with different levels of theory.

Our estimate of the uni-molecular rate coefficient for the reaction compares well with previous estimates using TST or RRKM theory (Hazra & Sinha, 2011; Long et al., 2012). Both methods have strengths and weaknesses. On the one hand, the static TST and RRKM calculations can be done with higher levels of theory and larger basis sets than we can handle with BO-MD and the currently available computational resources. It is also impossible, at the current state-of-the-art, to include the contributions of zero-point energy and tunnelling in the path sampling simulations. On the other hand, the kinetic theories are limited by their use of harmonic approximations to compute the vibrational partition functions which enter into the kinetic calculations (Long et al., 2012; Duncan, Bell & Truong, 1998). Path sampling, instead, is inherently dynamical, and the impact of vibrational anharmonicity is implicitly included. Being able to visualise the entire reactive landscape with path sampling, and not just the transition state itself, can also lead to new physical insights unavailable in traditional analyses based on TST or other kinetic theories. A more concerted attempt to carefully compare methodologies for reaction rate calculations would be a worthwhile subject for further study.

We want to emphasise that our results do not consider the complete SO₃ hydrolysis reaction, for example, the formation of the pre-reactive complex by atmospheric collisions. Our promising results encourage a broader application of the RETIS approach, which will be a goal pursued in forthcoming works. As the outset of this study, our stated focus was to demonstrate the effectiveness of path sampling to study a relatively simple uni-molecular gas-phase reaction. We acknowledge that in this case the uni-molecular reaction rate is much faster than the initial collision rate, and so overall H₂SO₄ formation in the atmosphere is not controlled by the uni-molecular reaction. The overall reaction was studied in great detail by previous researchers (Hazra & Sinha, 2011; Long et al., 2012), and our current results do not challenge their general conclusions.

The path sampling approach has here been successfully applied with BO-MD to study an atmospheric reaction. Our work highlights the advantage of studying relatively small gas-phase clusters compared with bulk aqueous phases, as the application of methods for studying rare events in combination with BO-MD is significantly eased in such small systems. We showed that path sampling can be exploited to efficiently generate ensembles of reactive paths and precise estimates of the reaction rate in gas-phase reactions of interest in atmospheric chemistry. The promising approach and results here reported are intended as an initial step towards a broader study of atmospherically relevant reactions using these novel sampling methods, including their application in our ongoing investigations into formic acid deprotonation (Murdachaew et al., 2016; Daub & Halonen, 2019).

Supplemental Information