Exploiting Cation Structure and Water Content in Modulating the Acidity of Ammonium Hydrogen Sulfate Protic Ionic Liquids

In this paper, we investigated the effect of cation structure and water content on proton dissociation in alkylammonium [HSO4]− protic ionic liquids (ILs) doped with 20 wt % water and correlated this with experimental Hammett acidities. For pure systems, increased cation substitution resulted in a reduction in the number of direct anion–anion neighbors leading to larger numbers of small aggregates, which is further enhanced with addition of water. We also observed spontaneous proton dissociation from [HSO4]− to water only for primary amine-based protic ILs, preceded by the formation of an anion trimer motif. Investigation using DFT calculations revealed spontaneous proton dissociation from [HSO4]− to water can occur for each of the protic ILs investigated; however, this is dependent on the size of the anion aggregates. These findings are important in the fields of catalysis and lignocellulosic biomass, where solvent acidity is a crucial parameter in biomass fractionation and lignin chemistry.


Hammett Acidity Measurements
IL Hammett acidity was measured using UV-Vis, combining the Beer-Lambert law (Equation 1) and a modified form of the Henderson-Hasselbalch equation (Equation 2) as in the work of Grasvik et al. 1
Where H0 is Hammett acidity,  "# ! is the basicity constant of the Hammett base,  ! is the absorptivity of the fully unprotonated Hammett base, and  is the absorptivity of the partially protonated Hammett base.
The Hammett base used in this work was 4-nitroaniline.This base has been previously reported to be suitable for the range of Hammett acidities investigated. 2The pKBH+ value of 4-nitroaniline is 1.00, obtained from literature. 1,3 he absorptivity of fully unprotonated 4-nitroaniline ( ! ) was determined experimentally to be 16,600 M -1 cm -1 , which is consistent with values previously reported by Grasvik et al. 1 This was done by measuring the UV-Vis absorbance at 380 nm of 5 different concentrations of the Hammett base in anhydrous dichloromethane, and measuring the gradient of the A vs. cl graph.
All other absorptivity values () were determined by measuring the UV-Vis absorbance of the unprotonated peak of 4-nitroaniline (380 nm) for 5 different concentrations of the Hammett base in the relevant ionic liquid and measuring the gradient of the A vs. cl graph.All Hammett acidity measurements were repeated twice, i.e. preparing new stock solutions of 4-nitroaniline in DCM and ionic liquid, which were then diluted to the desired concentrations.All UV-Vis measurements were performed using a Perkin Elmer Lambda 650, with solutions pipetted into sealed UV-clear quartz cuvettes with a path length of 0.5 cm.

Density Measurements
Aqueous IL solution densities were measured using a Mettler-Toledo DM40 density meter, between 0 o C and 50 o C. Solutions were sonicated for 20 minutes before measurements to remove any bubbles and dissolved air.The cell condition was verified before each use by measuring the density of purified water at 20 o C (0.99821 g/cm 3 ), and ensuring it was within the tolerance (0.00008 g/cm 3 ).

Classical Molecular Dynamics
Classical MD simulations were carried out using the Amber16 program. 4Structural and dynamical analysis of each trajectory was accomplished using the open-source programs TRAVIS, 5, 6 AGGREGATES 7 and CPPTRAJ. 8The protic ionic liquids (PILs) - and [DMBA][HSO4] -were modelled using the General Amber force field (GAFF) 9 and the four-site TIP4P-EW 10 water model was used in simulations that included water molecules.The GAFF potential has been recently applied to similar PILs systems, 11 and has been shown to be able to predict thermodynamic and transport properties of several ILs, both pure and when combined with water (i.e.doped or aqueous IL). 12 In addition, GAFF is recommended where limited experimental data is available. 12Furthermore, Tenney et al. have previously employed GAFF for the investigation of ILs and have noted that in general the initial GAFF parameterised assigned require very minor refinements. 13onsequently, we employ the GAFF parameters as assigned.Next, partial atomic charges for each of the PIL species were assigned using the RESP 14 (restrained electrostatic potential) method at the B3LYP/6-311++g(d,p) level of theory.This approach has been previously shown to yield reliable charges for organic molecules 15 and ILs. 13 As in the literature the charges of the PIL species were scaled by 0.8 to model the average charge screening due to polarisation and improve solvent dynamics. 12,13 h simulation box comprised of 400 cation-anion ion pairs.For simulations containing 20 wt% water, the number of water molecules added for each system is provided in Table S1.Initial simulation configurations were generated using the PACKMOL program. 16Each simulation has been carried out using periodic boundary conditions and a 1 fs time step.The particle mesh Ewald (PME) method was employed for long-range electrostatics interactions and a 12 Å cut-off was included for non-bonded interactions. 17The temperature and pressure were controlled using the Langevin thermostat (γ = 1.0 ps -1 ) and Berendsen barostat (p = 1 atm, τ = 1000 fs). 18,19 ll bonds containing hydrogen atoms were constrained using SHAKE. 20 equilibrate each system, we followed a procedure adapted from our previous work on IL mixtures. 21ollowing an initial minimisation (steepest descent -1000 steps and conjugated gradient -4000 steps), each simulation was subjected to the following steps.1) A heating phase, where each system was gradually heated to 600 K over 5 ns in the NVT ensemble to remove possible energy hotspots.2) A short (1 ns) equilibration was then carried out at 600 K in the NPT ensemble.3) This was followed by a cooling phase, where each system was gradually cooled from 600 K to 300 K over 5ns (NVT).4) A further 10 ns NPT simulation was carried out to obtain a consistent density.The average density for the last 5 ns of the each of the respective NPT simulations was then used to carry out 100 ns production simulations at 300 K in the NVT ensemble and the last 50 ns of each production simulation was used for analysis.Additionally, each system was characterised by computing the density at different temperatures ranging from 273 K to 353 K in steps of 10 K (except for 300 K).Computed densities of the pure PILs and PILs doped with 20 wt% H2O, together with experimentally measured densities for the doped PILs obtained during this study are reported in Figure S and Table S2.Moreover, for the pure PILs there is a lack of experimentally measured densities.Consequently, we have only included one data point for pure [BA][HSO4].Each of our simulated results were found to be within 2 % of those measured experimentally.Given the accuracy of these results and the great body of work carried out on ILs using the GAFF force field, we are confident with the selected model and the obtained results.
Table S1.The number of water molecules added to each classical molecular dynamics simulation to obtain a final concentration of 20 wt% H2O.Each simulation contains 400 PIL ion pairs.

Ab initio Molecular Dynamics
Ab initio molecular dynamics (AIMD) simulations were carried out using the QUICKSTEP 22 module within the CP2K 23 code.The BLYP 24,25 functional was employed with the molecularly optimized doubleζ basis set (MOLOPT-DZVP-SR-GTH) 26 for each of the atoms and core electrons were represented by the corresponding BLYP Goedecker-Teter-Hutter 27-29 pseudopotentials.Dispersion interactions were accounted for using the DFT-D3 empirical dispersion correction. 30,31 he resulting BLYP-D3 functional has been previously shown to produce good results for ionic liquid and hydrogen bonding systems. 32itial simulation boxes were generated using PACKMOL 16 and contained 16 ion pairs of the corresponding protic IL and the commensurate number of water molecules.The systems were preequilibrated using following the method described above for classical MD simulations.For the AIMD equilibration phase, the plane wave cutoff was set at 200 Ry, together with multigrids number 5 (NGRID 5 and REL CUTOFF 30).Smoothing of the electron density was carried out using NN10 SMOOTH and its derivative (NN10). 22The SCF convergence evaluation criteria was set to the default (10 -5 ) and the DIIS minimiser 23 was used to reach a faster orbital transformation via direct inversion in the iterative subspace.The maximum number of SCF iterations is 100 while a maximum of 10 iterations was performed for outer SCF loops.A further equilibration over 5 ps was performed using the keyword REGION MASSIVE, which means that for every single atom is thermostated individually for faster equilibration with a time constant of the thermostat chain of 50 fs.The temperature for this phase was set at 450 K to increase the sampling.Periodic boundary conditions were applied to avoid boundary effects.
A final production run was then carried out for each system in the NVT ensemble.For these runs, the density CUTOFF criterion was set to 350 Ry, together with NGRID 4 and REL CUTOFF 40.Moreover, the SCF convergence was set to 1.0E-6.A Nose-Hoover chain thermostat 33,34 with a time constant 100 fs was used to obtain a target temperature of 300 K.Each simulation was subsequently run for 100 ps (using a 0.5 fs time step) and the last 90 ps were used for analysis.

Quantum Chemical Calculations
DFT calculations were carried out using Gaussian 16 (version c01). 35Geometry optimisations were performed using the B3LYP 24,36 functional and the 6-311++G(d,p) basis set.To account for dispersion interactions, Grimme's -D3 30,31 dispersion correction with the Becke and Johnson dampening (BJdampening) function [37][38][39] was used.This combined method is referred to as B3LYP-D3BJ in the text.The combination of functional and basis set (B3LYP-D3BJ/6-311++g(d,p)) has been previously shown to provide good quality structures and energies for ionic liquid clusters.The Conductor-like Polarisable Continuum Model (CPCM) implicit solvent model was employed. 40,41 he dielectric constant of ethanol was selected to represent the protic IL-water environment, based on the literature. 42All structures have been fully optimised under no symmetry constraints and have been confirmed as minima using vibrational analysis.Optimisation convergence criteria were set to 10 -9 on the density matrix and 10 -7 on the energy matrix.The numerical grid was improved from the default to a pruned (optimised) grid of 99 radial shells and 590 angular points per shell.Vibrational frequencies and zero-point vibrational energy corrections (ZPE) were obtained within the harmonic approximation for each structure.Topological analysis of the electron density, within the quantum theory of atoms (QTAIM) 43 framework was carried out using AIMALL 44 and NBO analysis was carried out using NBO 6.0.Table S3: Energetic data for the lowest energy structures for the anhydrous PIL ion pair dimers in the gas-phase.These structures are labelled according to the orientation of the two [HSO4] -anionseither cyclic (AA -contains two anion-anion hydrogen bonds) or chain (AHA -1 x anion-anion hydrogen bond).See Figure S12 for a clearer view of these structural conformers.All energy differences are reported in kJ.mol -1 and DEZPE includes the ZPE correction.All of E, G, H, S and ZPE have been rounded to fifth decimal point.S3) for each of the anhydrous PIL ion pair dimers.

IL
Table S4: Energetic data for the lowest energy structures for the protic IL ion pair dimers containing 6 H2O molecules (~20 wt% H2O) and include the solvent environment via CPCM (ethanol).These structures are identified according to the orientation of the two [HSO4] -anions and waters -cyclic (AA -contains two anion-anion hydrogen bonds), chain (AHA -1 x anion-anion hydrogen bond) intercalated waters (IC -2 x waters between the anions) and proton transfer (PT).All energy differences are reported in kJ.mol

Figure S1 :
Figure S1: Variation of the Hammett acidity of [BA][HSO4] with excess butylamine.All IL measurements were conducted using 20 wt % water.Excess base was measured on a molar basis.

Figure S2 :
Figure S2: Variation of Hammett acidity of [DMBA][HSO4] with excess sulfuric acid.All measurements were conducted using 20 wt % water.Excess acid was measured on a molar basis.

Figure S3 :Figure S4 :
Figure S3: Variation of IL density with temperature.All IL measurements were conducted using 20 wt % water.Excess sulfuric acid was measured on a molar basis.

Figure S9 :
Figure S9: Size distributions of Anion-Anion aggregates for (a) [BA][HSO4], (b) [HA][HSO4], (c) [MBA][HSO4] and (d) [DMBA][HSO4], obtained from the classical molecular dynamics simulations of the respective anhydrous (blue) and doped with 20 wt % water (green) protic ILs.Size distribution data was computed using the AGGREGATES software, employing the corresponding S-S distance cutoff determined from the first minimum of the anion-anion radial distribution functions.

Figure S12 :
Figure S12: Representative lowest energy structures for the anhydrous Protic IL ion pair dimers.In these structures the two [HSO4] -anions are found to orientate in either a cyclic (AA -contains two anion-anion hydrogen bonds) or chain (AHA -1 x anion-anion hydrogen bond) conformation (see Figure S11).DEZPE (black) and DG (green) energies are all reported in kJ.mol -1 .

Figure S13 :
Figure S13: QTAIM and NBO data for the Anion-Anion H-Bonding interactions found for the lowest energy structure (AA_0 -TableS3) for each of the anhydrous PIL ion pair dimers. 45,46

Table S2 .
Computed and experimental densities for pure PILs and PILs doped with 20 wt % H2O determined in this work.

Table S5 :
-1and DEZPE includes the ZPE correction.All of E, G, H, S and ZPE have been rounded to fifth decimal point.Energetic data for the two low energy structures where proton transfer (X_PT0; X is the parent cation, i.e.BA, HA, MBA and DMBA) has occurred or where proton transfer has not occurred (X_0) for the protic IL ion pair trimers containing 9 H2O molecules (~20 wt% H2O) and include the solvent environment via CPCM (ethanol).All energy differences are reported in kJ.mol -1 and DEZPE includes the ZPE correction.All of E, G, H, S and ZPE have been rounded to fifth decimal point.