Predictive energetic tuning of C-Nucleophiles for the electrochemical capture of carbon dioxide.

This work maps the thermodynamics of electrochemically generated C-nucleophiles for reactive capture of CO2. We identify a linear relationship between the pKa, the reduction potential of a protonated nucleophile (E red ), and the nucleophile's free energy of CO2 binding ( Δ G b i n d ). Through synergistic experiments and computations, this study establishes a three-parameter correlation described by the equation Δ G b i n d = - 0.78 p K a + 4.28 E r e d + 20.95 for a series of twelve imidazol(in)ium/N-heterocyclic carbene pairs with an R 2 of 0.92. The correlation allows us to predict the Δ G b i n d of C-nucleophiles to CO2 using reduction potentials or pKas of imidazol(in)ium cations. The carbenes in this study were found to exhibit a wide range CO2 binding strengths, from strongly CO2 binding to nonspontaneous. This observation suggests that the Δ G b i n d of imidazol(in)ium-based carbenes is tunable to a desired strength by appropriate structural changes. This work sets the stage for systematic energetic tuning of electrochemically enabled reactive separations.

Introduction

Lowering atmospheric concentrations of CO2 to levels that will mitigate the effects of climate change will require substantial implementation of carbon-neutral and carbon-negative technologies, including carbon capture (Institute, 2020; National Academies of Sciences, Engineering, and Medicine, 2019). The development of optimal carbon capture materials requires a careful balance of the energetics of CO2 binding to the proposed capture material. However, CO2 binding free energies () are not easily experimentally accessible.

Among proposed materials for CO2 capture and conversion, imidazolium salts have gained considerable attention due to the attractive liquid properties of their alkylated variants, including their favorable solvation of CO2 gas as well as their high conductivity, which renders them uniquely suited for electrochemical applications (Kelemen et al., 2011; Sowmiah et al., 2009; Zhao et al., 2016). In addition to these favorable properties, imidazolium salts readily yield N-heterocyclic carbenes (NHCs) through deprotonation (Scheme 1B). Interestingly, electrochemical access to the free NHCs has been previously demonstrated as shown in Scheme 1A. Upon reduction, an imidazolium cation releases half of aH2 equivalent and becomes a free NHC (Gorodetsky et al., 2004). The free NHC then exhibits reactivity with CO2, as shown in Scheme 1C. Historically, the reaction of Scheme 1C was originally developed as a protection method for reactive NHCs (Voutchkova et al., 2007), and it has found applications in the synthesis of organometallics (Hahn and Jahnke, 2008) and catalysis (Luca and Fenwick, 2015; Luca et al., 2015; Riduan et al., 2009; Yang and Wang, 2014).

Scheme 1

Reactions of imidazolium cations

Generation of NHCs (Im) from their imidazolium precursors (A and B) and Im binding to CO2 (C) along with the three properties associated with these reactions: (A) E, (B) pK, and (C) ΔG. R1 and R2 substituents are described in Figure 1.

Figure 1

Imidazol(in)ium-based NHCs used in the correlation studies

Structures of each compound under study. See also Figure S1 and Table S1 for applicable experimental values associated with each NHC. Imidazolium counterparts (protonated NHCs) are abbreviated as "NHC"H in this work.

Using this reactivity, we now demonstrate a three-way correlation between the experimentally accessible reduction potential (E) of the imidazolium (Scheme 1A), the pK of the imidazolium (Scheme 1B), and the CO2 binding energy of the corresponding NHC (Scheme 1C). While the correlation between the pK of an imidazolium (Scheme 1B) and ΔG of the corresponding free NHCs (Scheme 1C) has been studied using density functional theory (DFT) on a set of 90 imidazoliums (Wang et al., 2018), it has not been related to the E of the imidazolium (Scheme 1A), which is a property that is routinely measured by standard electrochemical experiments (Gorodetsky et al., 2004). While the number of the imidazolium salts in this study is smaller (Figure 1), we map a landscape of CO2 binding free energies spanning 15 kcal, 20 pKa units, and 1.7 V of potential.

We investigate the relationship between E, pK, and ΔG associated with reactions A–C in Scheme 1. The parameters are related due to their mutual dependence upon the electron density at the apical imidazolium carbon. Basicity and nucleophilicity are intrinsically related, with basicity being a subset of nucleophilicity. All nucleophiles are Lewis bases, but the term nucleophilicity is utilized to describe the affinity of a molecule to not only protons but also other electron-deficient molecules such as carbon dioxide. In essence, our correlation of pKa with the free energy of CO2 binding ranks the nucleophilicity of an NHC against two types of electrophiles: protons and carbon dioxide.

Increasing the electron density of an imidazolium core causes it to have a more negative E (Scheme 1A); a rise in electron density at a site of deprotonation lowers the stability of the deprotonated species and thus increases the pK of the conjugate acid(Scheme 1B); and finally, a rise in electron density on a nucleophile corresponds to an increase in its nucleophilicity and thus its ability to attack an electrophile such as CO2 is enhanced (Scheme 1C). Although we describe these general relationships to assert that the E, pK, and ΔG of the twelve imidazol(in)ium cations presented in Figure 1 are correlated, analogous relationships generally apply to other potential molecular CO2 capture substrates and thus the approach we use herein may be applied to other, more complex systems for redox-enabled CO2 capture, such as quinones and N-nucleophiles (Wilcox, 2020). Unlike traditional electroswing systems which operate through cycles of binding CO2 with cathodically generated nucleophiles followed by release of CO2 in an electrochemical oxidation step, the electrochemical step of NHC generation and binding of CO2 can be coupled with a thermal release step rather than an oxidation. The thermal release of CO2 from NHC adducts is known (Van Ausdall et al., 2009; Luca and Fenwick, 2015), thus setting the stage for hybrid electrothermal CO2 capture-release methodologies.

Results

Before performing DFT calculations to investigate the relationship between free energy of CO2 binding, reduction potential, and pKa for the set of imidazolium cations, we first sought to experimentally validate several assumptions about the system. Because imidazolium cations ImH+ are aromatic, the reduction of ImH+ to add an electron to the imidazolium π-system results in de-aromatization of the imidazolium ring. Prior studies have shown that aromatic hydrides exhibit a linear correlation between the first reduction potential of the oxidized hydride and the nucleophilicity of the hydride (Alherz et al., 2018; Ilic et al., 2018). To verify this correlation for imidazolium salts, we measured E for a family of imidazolium salts (Figure S1 and Table S1) using normal pulse voltammetry (NPV). While the influence of ion pairing and association on properties of imidazolium salts has been well documented in the ionic liquid literature (Tsuzuki et al., 2013; Voroshylova et al., 2018; Wang et al., 2015), the nature of the counterion's effects on our measurements is minimized through the use of NBu4PF6, an electrolyte with a chemically stable, non-coordinating PF6− anion. Furthermore, we validated the observation that one-electron reductions of imidazolium cations produce free NHCs and verified that this applies to the analogous imidazolinium cations with non-aromatic, saturated backbones. We confirmed that the voltammetry responses of imidazolinium cations are associated with one e− reductions using the method of Donadt et al. (Donadt et al., 2018) based on diffusion ordered spectroscopy (DOSY NMR) and NPV to quantify the number of electrons associated with the assumed one-electron reduction step (Scheme 1A and Figure S2); this is consistent with the observed one e− reduction of imidazolium cations to produce free NHCs.

With these experimental validations in hand, we next computed the values of E, pK, and ΔG for each of the imidazolium derivatives using DFT based on the MN15/6-311+G(d,p) level of theory paired with the SMD solvent model as implemented in the Gaussian 16 quantum chemistry code (Frisch et al., 2016; Krishnan et al., 1980; Marenich et al., 2009; McLean and Chandler, 1980; Yu et al., 2016). The DFT-computed pK and E values we obtained agree with experiment (Chu et al., 2007; Dunn et al., 2017), with mean absolute errors (MAE) below 0.15 pKa units and 0.05 V, respectively. For additional computational details, please refer to the STAR methods, Figure S3 and S4, Tables S2–S4.

The calculated free energies of binding to CO2 range from +5 to −10 kcal/mol, indicating that the equilibrium constant for CO2 binding varies by over 10 orders of magnitude at 298 K. Consequently, this energy window suggests that the thermodynamics of NHC binding to CO2 can be tuned over a wide range, for instance, by functionalization of NHCs with appropriate substituents to yield carbenes that bind CO2 weakly or even nonspontaneously (e.g. 2CN-IMes) to strongly (e.g. Butyl-methyl imidazolium BMIM). Notably, our calculations predict that four carbenes known to be stable in the presence of oxygen (2Cl-IMes, 2Br-IMes, 2Cl-IPr, and 2Br-IPr) (Arduengo et al., 1997; Furfari et al., 2015) have slightly negative values of ΔG (−1.60 kcal/mol, −0.43 kcal/mol, −0.76 kcal/mol, and −0.50 kcal/mol, respectively). This observation has implications for direct air capture (DAC) of carbon dioxide, where stability to oxygen is a critical prerequisite for candidate materials. Weak CO2 binders have an advantage for DAC processes whose thermodynamic efficiency also depends on low-energy release of CO2 from the sorbent material (Lackner, 2013; Petersen and Luca, 2021).

Computations were then used to investigate three sets of correlations based on reactions a–c of Scheme 1 for twelve imidazol(in)ium cations. Figure 2 includes imidazolium cations that are not easily synthetically accessible. In these cases, correlations were performed using the calculated DFT values for E and pK.

Figure 2

Linear correlations of pKa, free energy of CO2 binding, and reduction potential

The three linear two-way correlations of (left) reduction potential and pK, (middle) free energy of CO2 binding and pK, and (right) reduction potential and free energy of CO2 binding for the compounds under study. Optimized molecular coordinates are available in Data S1.

For imidazol(in)ium cations with known pK or E, estimated ΔGbind for the corresponding NHC are readily obtained using the equations produced by the linear regressions shown in Figure 2. For instance, the ΔG in kcal/mol of the commercially available NHC precursor IMesH is estimated from either Equations 1 or 2:

In the case of Equation 1, an experimental E for the compound is obtained within minutes using a routine electrochemical experiment such as NPV. Using the E value we measured of −2.64 V vs. Fc/Fc+ for IMesH, we then obtained an estimated ΔG for CO2 binding from Equation 1 of −5.94 kcal/mol. Similarly, with a known pK value of IMesH in DMSO of 19.4(Wang et al., 2018), Equation 2 is used to estimate a ΔG of −4.87 kcal/mol. These values compare favorably with the DFT-calculated ΔG of −5.78 kcal/mol, yielding errors of 0.16 kcal/mol (3%) and 0.91 kcal/mol (15%), respectively, and bracketing the calculated value.

Alternatively, in the case of imidazol(in)ium cations for which both the experimental pK and E are known, Equation 3 (Figure 3),is used to produce an estimated value for ΔG with higher confidence than the individual correlations presented in Figure 2. For IMesH, this results in an estimated ΔG of −5.44 kcal/mol. This again compares favorably to the DFT-calculated value of ΔG, this time yielding an error of only 0.34 kcal/mol (6%) and improving substantially upon the value obtained from the pK-only estimation. The results of using Equation 3 for all six compounds for which experimental values of E and pK are reported in Table 1.

Figure 3

3D correlation provides greater predictive power for estimating the free energy of CO2 binding

Three-parameter linear relationship providing stronger correlation between Ered, pKa, and ΔGbind for the compounds under study

Table 1

Comparison of DFT-calculated ΔG values with those estimated using experimental reduction potentials and pK value

ImH	Exp. Ered (V vs. Fc/Fc+)	Exp. pKa	Est. ΔGbind (kcal/mol)	DFT ΔGbind (kcal/mol)	Error (%)
IMesH	−2.64	19.4	−5.4	−5.8	6
IPrH	−2.63	19.3	−5.3	−6.0	12
IMeH	−2.83	22.0	−8.3	−8.3	0.2
EtMeH	−2.87	22.1	−8.5	−8.2	4
BMIM	−2.84	22.0	−8.3	−8.4	1

Finally, as reduction potential is the most experimentally accessible value, it is desirable for various applications to predict the pK of an imidazol(in)ium cation from its measured E. To do this, the two-way correlation equation of Figure 2 right is used (Equation 4),

Table 2 below compares computed and experimental pKas.

Table 2

Comparison of estimated pK values from Equation 4 and experimental pK values

ImH	Exp. Ered (V vs. Fc/Fc+)	Est. pKa	Exp. pKa	Error (%)
IMesH	−2.64	22.2	19.4	15
IPrH	−2.63	22.2	19.3	15
IMeH	−2.83	23.8	22.0	8
EtMeH	−2.87	24.1	22.1	9
BMIM	−2.84	23.9	22.0	8

Discussion

We now establish a three-way correlation between the reduction potential, pK, and free energy of CO2 binding for a set of imidazol(in)ium cations and their corresponding NHCs and use it to demonstrate the facile calculation of CO2 binding energies from experimentally accessible electrochemical measurements. Four halogenated NHCs known to be stable to oxygen, a necessary sorbent trait for carbon capture applications in which oxygen may be present, were predicted to bind CO2 with negative values of ΔG. Thermodynamic calculations for carbon capture (Petersen and Luca, 2021) suggest that materials that target energetically efficient reversible binding of CO2 may be weak binders because strong binding results in a substantial energy penalty for the release phase. Although in this study we demonstrate our predictions on imidazol(in)ium-based NHCs, similar correlations will apply to other families of electrochemically responsive CO2 capture molecular materials, thus enabling the mapping of energetics of a new chemical space for this application. Coupled with the known thermally triggered release of CO2 from NHC-CO2 adducts, our findings set the stage for energetic tuning of reactive separations in novel hybrid redox-thermal swing chemistries.

Limitations of the study

Strong ionic association may have an impact on observed reduction potentials, resulting in values that differ substantially from those predicted by the correlations presented herein. While the subset of organic salts investigated is limited, we expect the correlation to hold across a large library of chemically similar species. However, steric effects may cause significant deviations from the established trends, as has been observed in the case of 1,3-di-t-butylimidazolium chloride (ItBuHCl).

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Oana Luca (oana.luca@colorado.edu).

Materials availability

This study did not generate new unique reagents.

Method details

Electrochemical experiments

General methods

Reagents for use in electrochemical experiments were purchased from commercial sources and used as received except where otherwise noted. 1,3-bis(2,4,6-trimethylphenyl) imidazolium chloride (IMesHCl), 1,3-Bis(2,4,6-trimethylphenyl)-4,5-dihydroimidazolium chloride (SIMesHCl), 1,3-bis(2,6-diisopropylphenyl)imidazolium chloride (IPrHCl), 1,3-bis(2,6-diisopropylphenyl)-4,5-dihydroimidazolium chloride (SIPrHCl), 1,3-bis(2,6-diisopropylphenyl)-4,5-dihydroimidazolium tetrafluoroborate (SIPrH BF4), 1,3-dimethylimidazolium chloride (IMeHCl), 1-ethyl-3-methylimidazolium hexafluorophosphate (IEtMeHCl), 1,3-di-t-butylimidazolium chloride (ItBuHCl), and 1-butyl-3-methylimidazolium bis(trifluoromethylsulfonyl)imide (BMIM TFSI) were purchased commercially. Acetonitrile (MeCN) used in electrochemical experiments was freshly distilled and sparged with argon before use. Tetrabutylammonium hexafluorophosphate (NBu4 PF6) was recrystallized from methanol and dried overnight under high vacuum before use as supporting electrolyte.

Normal Pulse Voltammetry (NPV) experiments were performed using a silver wire single junction pseudo-reference electrode, a 3 mm diameter glassy carbon working electrode (MF-2012, BASI), and a platinum wire counter electrode. Standard NPV experimental conditions were 50 ms pulse width, 200 ms step time, and −10 mV pulse height. A vial filled with MeCN was used as a pre-bubbler to prevent solvent evaporation over the course of the experiment. Each imidazolium salt was dissolved in a 0.1 M NBu4 PF6 solution in MeCN and was degassed with argon prior to analysis. All glassware and needles used in the experiment and in solution transfer were oven-dried prior to use, and solvents and supporting electrolyte were dried as described above. The ferrocene/ferrocenium couple was used as an external standard reference for the applied potential.

Experimental reduction potentials

Each imidazolium salt was dissolved as a 1 mM solution in acetonitrile with 0.1 M NBu4 PF6 and subjected to the standard NPV conditions described above. The chloride salt was used for each imidazolium of interest, with the exception of BMIM TFSI. Experimental imidazolium reduction potentials were obtained by finding the potential at which the inflection point occurs in the NPV of the compound in question. This method was chosen over standard use of cyclic voltammetry E1/2 values due to the irreversible nature of the reduction and the improved ability of NPV to provide precise current responses with minimal interference from charging currents. This value was then averaged across three trials (with the exception of SIPr, 6 trials).

Electron quantification

NPV was used to determine the number of electrons associated with the reduction of SIPr using the method of Donadt et al. The number of electrons associated with a given electrochemical feature can be calculated using the Cottrell Equation, , in which il is limiting current (the current at the plateau of a feature), n is the number of electrons associated with the feature, F is Faraday’|'s constant, A is the electrode surface area, Co is the bulk concentration of the analyte, D is the diffusion coefficient of the analyte, and t is the pulse width. The diffusion coefficient of SIPrH was first obtained by conducting Diffusion-Ordered SpectroscopY (DOSY) NMR on a sample of SIPrH BF4 in deuterated MeCN. The diffusion coefficient was found to be 1.11∗10−5 cm2/s. NPV (Figure S2) was then collected using a sample of 3.53 mg of the imidazolium salt IPrH BF4 dissolved in a known amount of MeCN. The total sample mass was obtained before and just after the experiment to account for solvent loss during degassing and was used to calculate a final bulk concentration of 1.68 mM. Using the geometric electrode surface area of the 3 mm glassy carbon working electrode and the respective time parameter for each scan, the value of n can then be obtained with the Cottrell Equation.

Computational data

Quantum chemical density functional theory calculations have been conducted to compute the optimized geometries and energies of the carbene precursor, its protonated and reduced intermediates, and their adducts with CO2. Frequency calculations were performed to ensure that the optimized geometries have no imaginary frequencies and obtain thermal and entropic corrections to the electronic energy to calculate Gibbs free energies. The DFT calculations utilize the Minnesota 15 (MN15) functional, Pople basis sets with polarization and diffuse functions 6-311+G(d,p), and acetonitrile solvent as described by the SMD solvent model (Krishnan et al., 1980; Marenich et al., 2009; McLean and Chandler, 1980; Yu et al., 2016) implemented within the Gaussian 16 software (Frisch et al., 2016). The energies obtained are used to calculate 3 main observable properties: pKa values, reduction potentials (of the protonated carbene intermediate), and the free energies of binding CO2 to the resultant carbenes.

Similar computations were performed using the wB97xD/6-311+G(d,p)/CPCM method for comparison (Barone and Cossi, 1998; Chai and Head-Gordon, 2008; Cossi et al., 2003). The default cavity settings were additionally modified, as recommended by Mirzaei et al. (Mirzaei et al., 2019), such that the atomic radii for cavity estimates are modelled according to the Bondi radii rather than the fitted SMD-Coulomb parameters, and the calculations were carried out again. The three overall computational methods compared herein are therefore MN15/6-311+G(d,p)/SMD, MN15/6-311+G(d,p)/SMD-Bondi, and wB97xD/6-311+G(d,p)/CPCM.

All three DFT methods investigated are rather accurate in predicting pKa values compared to experimental values obtained from Ref. 18, as shown in Table S2. While the pKa and Ered values of ItBuH were accurately predicted, this imidazolium was omitted from the correlation plots because it deviates significantly from the trends obeyed by the other species where its computed ΔGbind value deviates by more than 6 and 7 kcal/mol from the expected values provided by the Ered and pKa correlations, respectively. This significantly different behavior most likely stems from the bulky tBu groups of ItBuH that hinder CO2 binding more than planar Ph-based or smaller alkyl groups, resulting in a much less favorable ΔGbind value at comparable Ered and pKa values that are not shifted so dramatically because electron and proton transfers are not hindered as much as CO2 binding.

While the SMD-Bondi cavity treatment has not been benchmarked rigorously in literature, we find that it succeeds in providing accurate pKa and Ered values (Tables S2 and S3), probably due to the nature of our fitting techniques discussed below. However, it seems to significantly over-stabilize the binding energies of CO2 to carbenes. The SMD-Bondi cavity treatment significantly exaggerates the likelihood of binding (Table S4). Since this cavity method is the least tested across different systems and thus least trustworthy, we concluded that the default SMD treatment with the MN15 functional provides more realistic binding free energy values and opted to use this method for our correlations. It is also more in agreement with the well-established wB97xD/CPCM method, within expected DFT errors (MAE = 3.2 kcal/mol).

Acidity (pKa)

pKa values are obtained by creating a linear fit between experimentally determined pKa values and the computed free energy difference – . While it is common to estimate pKa values using an isodesmic approach or use a pre-determined free energy of solvated protons (Alherz et al., 2018; Thapa and Schlegel, 2016), we argue that a linear fitting approach avoids any systematic errors in pKa predictions. Figure S3 shows the linear correlation described. The linear correlation is used to predict pKa values for other carbenes by computing the energies of the XH+ and X intermediates. This computation is performed for the energies obtained from each of the three DFT levels of theory for each imidazolium (Table S4). This can also be used to extrapolate the expected free energy of an acetonitrile-solvated proton, , in close agreement with the experimentally determined value of −265.9 kcal/mol (Camaioni and Schwerdtfeger, 2005; Isse and Gennaro, 2010; Kelly et al., 2006; Marenich et al., 2014).

Reduction potentials (Ered)

Reduction potentials are first calculated according to the equation where is the absolute redox potential of the ferrocene/ferrocenium reference couple (4.988 V in acetonitrile) (Namazian et al., 2010; Tossell, 2011). To avoid systematic underestimation of Ered (−0.3 V vs. experimental values), we then utilize a linear fitting approach between experimentally determined Ered values and the free energy differences . This computation is performed for the energies obtained from each of the three DFT levels of theory for each imidazolium (Table S4).

CO2 binding free energy

The free energy of CO2 binding is computed according to the equation , where is the free energy of the NHC-CO2 adduct, is the free energy of the free NHC, and is the free energy of free CO2. This computation is performed for the energies obtained from each of the three DFT levels of theory for each imidazolium (Table S4).

Coordinates

Data S1: Optimized molecular coordinates. All optimized molecular coordinates (MN15/6-311+G(d,p)/SMD level of theory) for the carbenes and intermediates of interest used in Figure 2 are provided in.xyz format as a supplemental data.

Correlation

Best-fit analyses of each of the correlated parameters (Ered vs. pKa, Ered vs. ΔGbind, pKa vs. ΔGbind, and Ered vs. pKa vs. ΔGbind) were carried out in OriginPro 2019b.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Chemicals, peptides, and recombinant proteins
1,3-bis(2,4,6-trimethylphenyl)-imidazolium chloride (IMesHCl)	Strem Chemicals	Cat#07-0299; CAS-No. 141556-45-8
1,3-Bis(2,4,6-trimethylphenyl)-4,5-dihydroimidazolium chloride (SIMesHCl)	Strem Chemicals	Cat#07-4011; CAS-No. 173035-10-4
1,3-bis(2,6-diisopropylphenyl)-imidazolium chloride (IPrHCl)	Strem Chemicals	Cat#07-0590; CAS-No. 250285-32-6
1,3-bis(2,6-diisopropylphenyl)-4,5-dihydroimidazolium chloride (SIPrHCl)	TCI Chemicals	Cat# B3157; CAS-No. 258278-25-0
1,3-bis(2,6-diisopropylphenyl)-4,5-dihydroimidazolium tetrafluoroborate (SIPrH BF4)	Sigma Aldrich	Cat# 693553-1G; CAS-No. 282109-83-5
1,3-dimethylimidazolium chloride (IMeHCl)	TCI Chemicals	Cat#D3341; CAS-No. 79917-88-7
1-ethyl-3-methylimidazolium hexafluorophosphate (IEtMeHCl)	Alfa Aesar	Cat# L19762; CAS-No. 155371-19-0
1,3-di-t-butylimidazolium chloride (ItBuHCl)	Strem Chemicals	Cat#07-0368; CAS-No. 157197-54-1
1-butyl-3-methylimidazolium bis(trifluoromethylsulfonyl)imide (BMIM TFSI)	Sigma Aldrich	Cat#77896-1G-F; CAS-No. 174899-83-3
Tetrabutylammonium hexafluorophosphate	Sigma Aldrich	Cat#86874-100G; CAS-No. 3109-63-5
Acetonitrile	Pharmco	Cat#300000DIS; CAS-No. 75-05-8
Methanol	Macron Fine Chemicals	Cat#3016-16; CAS-No. 67-56-1
Software and algorithms
MN15 Functional	Yu et al., 2016	N/A
6-311+G(d,p) Basis Set	Krishnan et al., 1980	N/A
SMD Solvent Model	Marenich et al., 2009	N/A
McLean-Chandler (12s,9p) → (621111,52111) basis sets	McLean and Chandler, 1980	N/A
ωB97X-D Functional	Chai and Head-Gordon, 2008	N/A
CPCM Solvent Model	Barone and Cossi, 1998	N/A
Gaussian 16	Frisch et al., 2016	N/A
Optimized molecular coordinates	Data S1	Used in Figure 2
OriginPro 2019b	OriginLab	N/A