Ionic Liquids

All samples were imaged in a JCM-5700 Scanning Electron Microscope (JEOL USA, Peabody, MA, USA) contained inside a mobile biological containment enclosure (Dycor Technologies Ltd, Edmonton, AB, Canada)8 (link). Gold coated specimens were imaged under high vacuum at 6 kV, with an 8 mm working distance and a 30 μm objective lens aperture. Images were collected using the secondary electron detector, the acquisition time per image was 160 sec and each image was 2560 × 1920 pixels. Images of ionic liquid stained samples were obtained using the above noted settings with the exception that the acceleration voltage was adjusted to 4 kV. SEM images were recorded at magnifications ranging from 3,000 x to 20,000x.

Thermogravimetric measurements were carried out using thermogravimetry/differential scanning calorimetry TGA/DSC1 analyzer (Mettler Toledo, Greifensee, Switzerland) in the temperature range of 25–600 °C, with a heating rate of 10 °C/min in an argon atmosphere (flow rate 50 ml/min.). Prior to the measurements, Thermogravimetry (TG) analyzer was calibrated using indium and zinc as standards. Additional analysis was performed using Setsys TG-DTA 16/18 analyser (SETARAM Instrumentation, Caluire-et-Cuire, France) coupled to a Balzers (Pfeiffer) mass spectrometer for evolved gas analysis.
DSC1 analyzer (Mettler Toledo, Greifensee, Switzerland), calibrated with indium and n-octane as standards, was employed to study thermal transitions of pure and DmiBr-modified fillers and the temperature of ionic liquid release/desorption from the surface of filler. The measurements were performed in the temperature range of 25–500 °C, with a heating rate 5 °C /min.
Rubber compounds of ethylene-propylene-diene elastomer (EPDM, Vistalon 8600, Exxon Mobile, Irving, TX, USA) containing 20 phr of DmiBr-modified fillers were prepared using a laboratory two-roll mill. Then, the prepared EPDM compounds were cured at 150 °C using an electrically heated hydraulic press for the optimal vulcanization time, which was determined with rotorless D-RPA 3000 rheometer (MonTech, Buchen, Germany).
SEM images of analyzed filler surface and fractures of EPDM vulcanizates were taken using an LEO1450 SEM microscope (Carl Zeiss AG, Oberkochen, Germany). Prior to the measurement, vulcanizates were broken down using liquid nitrogen; their fractures were coated with carbon and next examined. Based on the SEM images, the morphology and size of filler particles were studied, as well as their dispersion in the elastomer matrix. Energy-dispersive X-ray spectroscopy (EDS) was used to confirm the presence of DmiBr on the surface of modified fillers. Samples of pure fillers were coated with carbon to improve the quality of SEM/EDS results.

Protex augments an OpenMM simulation object and is not restricted to simulations of ionic liquids. The two main parts of the program are the ProtexSystem and Update classes. The former gathers the simulation object and additional information on the update process, wrapped in the ProtexTemplates class. The latter is responsible for the actual update process and handles the logic during an update. Figure 3 gives an overview of the program package protex.
The system object was created using CHARMM topology and parameter files in this work. A condition to perform proton exchange between residue protonation states is that the residues prone to a proton exchange have a one-to-one mapping between their atoms in the protonated and deprotonated form of the topology file. Please find a detailed example in the documentation at GitHub (https://github.com/cbc-univie/protex).
The ProtexTemplates class is used to gather the additional information needed for the simulation. The user may specify which transfer reactions should occur by specifying the residue names, the maximum distance, and the probability of this reaction. This way, the back-and-forth reaction of, for example, Im₁H⁺ + OAc⁻ → Im₁ + HOAc, can be defined independently of the reaction Im₁ + HOAc → Im₁H⁺ + OAc⁻. Additionally, the atom name of the donor/acceptor atom needs to be specified. This would be the hydrogen for Im₁H⁺ and the nitrogen for Im₁ or the hydrogen of the acetic acid and both oxygens of the acetate. An example for the concrete definition of these variables can be found in the SI.
The ProtexSystem class combines the two former objects. It serves as an anchor for the actual propagation of the simulation, stores all information on the individual molecules (e.g., current name, charges, parameters, … ) in a separate Residue class, and can be used for loading and saving the current state and a PSF file. Two additional reporters are available, one reporting the current charge of all molecules in the system (ChargeReporter) and one reporting the energy contributions of the individual force objects (EnergyReporter). They can be used similarly to any other OpenMM reporter.
The Update classes handle everything connected to the update process during the simulation. The abstract base class Update serves as an anchor for different concrete implementations. NaiveMCUpate was used in this study, which checks for updates based on the distance and probability criterion. If the distance between the acceptor and donor falls below the distance criterion (as defined in ProtexTemplates), the proton exchange will happen with the given probability. The StateUpdate is responsible for the actual updates. It can be called anytime during the propagation of the trajectory. The update can either happen instantaneously between protonation states or using a non-equilibrium protocol in which multiple intermediate λ-states are used to interpolate between a source and target protonation state smoothly. The user can specify if only partial charges or all non-bonded and bonded interactions should be changed between protonation states.
As found in our previous study, the equilibrium for the Im₁H⁺/OAc⁻ system is around 30% charged and 70% neutral species. Hence an optional mechanism to stay around this equilibrium was implemented. As reported by Lill and Helms (Lill and Helms, 2001 (link)), the energy barrier for (de-)protonation is a function of the local environment and is not restricted to the exchanging molecules. Strictly speaking, the position of the barrier maximum is also a function of the local environment (Lill and Helms, 2001 (link)). However, as the corresponding calculations result in significant computational effort, we start with a fixed distance criterion. Dreßler et al, (2020a) (link) and Dreßler et al. (2020b) (link) introduced a Fermi function based to model the probability as a function of the distance, which will be included in future versions of protex.
The current probability p_ref is updated at each proton exchange event (see Figure 4) $p=p_{\mathit{ref}}+c\cdot {\left(\frac{n_k^{\mathit{now}}}{n_k^{\mathit{ref}}}-1\right)}^3$ where $n_k^{\mathit{now}}$ and $n_k^{\mathit{ref}}$ are the current and reference (initial) number of molecules of species k and c is a tunable prefactor. The power of three ensures the sign stays the same and allows for increased or decreased probabilities p: A ratio $n_k^{\mathit{now}}/n_k^{\mathit{ref}}$ below unity indicates that the number of the corresponding species k is below average. Hence, a reaction of that species should occur less often, which is realized by the reduced probability of this reaction due to the negative bracket in Eq. 4. On the other hand, more molecules than the reference indicate too few reactions. Hence the positive factor increases the probability of the reaction. Protex is designed to model proton transfers in a solvent at room temperature. Quantum effects at lower temperatures may only be indirectly modeled by changing the distance criterion and probability for particular reactions.
Figure 4 shows the typical workflow of a protex simulation. Each number depicts the trajectory of one species in the system. After some specified simulation time (A), protex checks for possible proton transfers and executes them (indicated by the black arrows in Figure 4). Then the simulation is propagated until the next update event (B). Here, some of the molecules may have stayed close to each other and exchanged the proton back (see trajectory (7) and (8) in Figure 4). However, it is also possible that the proton is transferred to the next molecule [see trajectory (3)–(4)–(5)]. A significant amount of molecules never face a proton exchange [see trajectory (1), (6), (9), and (10)] which may be due to unfavorable orientations or no corresponding partner. The number of protonations equals the number of deprotonations, as the overall system is neutral at all times. Consequently, the number of up arrows is the same as that of down arrows in Figure 4. Also, the total number of protonations/deprotonations may differ between two exchange events. For example, (C) in Figure 4 has fewer exchanges than (A) or (B).
Benchmark tests on a NVIDIA RTX3090 and AMD Threadripper with a typical setup of 10 ps simulation time between the updates, showed that the protex routine takes about 25% of the total simulation time. Details can be found in the SI.

Details on the parametrization process of the molecular species involved can be found in Joerg and Schröder (2022) (link). In short, the force field for Im₁H⁺OAc⁻ was based on the CHARMM General Force Field (CGenFF) (Kumar et al., 2020 (link)). Since the ionic liquid is not fully featured in the standard force field, electrostatic and bonded parameters were optimized based on quantum-mechanical reference calculations. For the calculation of dynamics properties, polarizable MD simulations were utilized. The polarizability was implemented using the Drude model, which adds an additional harmonic spring to all non-hydrogen atoms to emulate the induced forces. Due to their low mass, hydrogen atoms cannot be made polarizable, so the respective polarizabilities are added to their corresponding parent atoms. Drude particles were assigned a mass of 0.4 μ and a force constant $k_{i\beta }^{\delta }$ = 1,000 kcal/mol/Ų, (squared Angstrom). For stability reasons, the maximum distance for the mobile Drudes was set to 0.25 Å. Lennard-Jones interactions were reduced as described in Joerg and Schröder (2022) (link), using Eq. 3. Scaling factors s of 0.25 and 0.4 were employed, each with five replicas and a simulation time of 50 ns Each system contained 1,000 molecules, resulting in 150 Im₁H⁺/OAc⁻ and 350 Im₁/HOAc each (representing the initial 30%:70% equilibrium) as shown in Table 3.
Packmol (Martínez et al., 2009 (link)) was used to pack the initial simulation boxes, which were subsequently subject to energy minimizations using CHARMM, removing possible clashes or very unfavorable configurations of molecules (Brooks et al., 2009 (link)). Then, the polarizable system was equilibrated with OpenMM for 5 ns applying a Monte-Carlo barostat at 1.0 atm to determine the final box length. The production runs in the NVT ensemble were done in OpenMM with a time step of 0.5 fs for 50 ns Temperature control of polarizable systems with the conventional Dual-Nosé-Hoover thermostat (Martyna et al., 1992 (link)) is challenging, due to heat flow from the degrees of freedom of real atoms to Drude atoms. This causes the center-of-mass temperature to be overestimated. Hence, we used a temperature-grouped Dual-Nosé-Hoover thermostat as described by Son et al. (2019) (link) and Gong and Padua, (2021) (link), which adds an additional group for center-of-mass translations, thus improving the accuracy of the simulations. The temperature was set to 300 K for the real atoms and 1 K for the Drude particles. Electrostatic interactions were treated using the Particle Mesh Ewald method: The cut-off distance was set to 11 Å and the switch distance to 10 Å. All simulations were run on the CUDA platform in single precision. Further details on the setup can be found in Joerg and Schröder, (2022) (link).
Four possible transfer reactions were defined, including the forward and backward reaction described by Eq. 1 as well as the transfer between Im₁H⁺/Im₁ and HOAc/OAc⁻. In this work, the protonation states were switched instantaneously, with no additional λ states between the initial and final state. In the first step at each transfer event (see Figure 4), distances between transferable hydrogen atoms and hydrogen acceptors (nitrogen/oxygen) of the other molecules were checked, and only those pairs with a distance lower than 1.55 Å considered for the next step. The second step involves proton transfers with a particular probability. The initial probability of Table 4 are in accordance with Jacobi et al. (2022) (link) but are updated applying Eq. 4. The time interval between the transfer checks was set to 10 ps