Polarizable Molecular Dynamics Simulations of Ionic Liquid Proton Transfer

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