Protex—A Python utility for proton exchange in molecular dynamics simulations

Abstract

Protex is an open-source program that enables proton exchanges of solvent molecules during molecular dynamics simulations. While conventional molecular dynamics simulations do not allow for bond breaking or formation, protex offers an easy-to-use interface to augment these simulations and define multiple proton sites for (de-)protonation using a single topology approach with two different λ-states. Protex was successfully applied to a protic ionic liquid system, where each molecule is prone to (de-)protonation. Transport properties were calculated and compared to experimental values and simulations without proton exchange.

1 Introduction

Molecular dynamics (MD) simulations have become indispensable in modern computational science. Over the last decades, major improvements have been made regarding the size and speed so that nowadays, biologically relevant systems [i.e., membrane proteins (Goossens and De Winter 2018)] and many others can be studied in acceptable timescales (Hospital et al., 2015). Polarizable MD simulations further improved the accuracy of the underlying force fields, especially for dynamic properties (Schröder and Steinhauser, 2010; Schröder et al., 2011; Bedrov et al., 2019).

One drawback of classical force fields is the fixed topology, which means bonds are not designed to build or break. There are different approaches how to deal with that: Reactive force fields (REAX-FF) (Russo Jr and Van Duin, 2011; Zhang et al., 2014; Weismiller et al., 2015) have been developed, which use bond orders to describe the formation or breaking of bonds. In condensed-phase system, proton transfer has been modeled by applying a Markov model on top of molecular dynamics simulations (Dreßler et al., 2020a; Dreßler et al., 2020b) Alternatively, alchemical mutations with single or dual topology approaches can be applied if topology changes are required (Boresch and Karplus, 1999a; Boresch and Karplus, 1999b; Shirts, 2012; Mey et al., 2020). Alchemical approaches typically use an alchemical coupling parameter λ to control the transition of one molecule into another one (including possible bond break/formation); in our case, the transition from the protonated to the deprotonated species or vice versa. In constant pH simulations, the (de-)protonation reaction may be described as an instantaneous protonation state change (Mongan et al., 2004) or using alchemical intermediates (Lee et al., 2004; Khandogin and Brooks, 2005; Mongan and Case, 2005; Radak et al., 2017; Dobrev et al., 2020).

However, almost all these approaches are usually applied to a solute with few (de-)protonation sites. Often these sites are coupled to a “proton bath” (Börjesson and Hünenberger, 2001; Donnini et al., 2016; Radak et al., 2017) or an implicit solvent (Mongan et al., 2004; Mongan and Case, 2005) to ensure charge neutrality of the simulation box. However, this coupling limits the number of (de-)protonation sites, which is fine for constant pH simulation of an aqueous protein solution but maybe not be appropriate anymore if all solvent molecules are subject to the proton transfers. This is particularly true for protic ionic liquids (PILs), as proton transfers must be adequately captured even though hundreds of (de-)protonation sites exist. PILs consist of a Brønsted base (B) and acid (HA) and, therefore, can exchange a proton, which is a reaction currently not featured by modern force fields. In general, this reaction reads as

For example, accounting for proton exchange effects will be crucial for an adequate description of the conductivity in PILs. Additionally, examining the moving proton particularly can gain insight into the mechanism. However, classical constant pH simulations cannot model proton hopping from one molecule to another. Multistate empirical valence bond models for water (Schmitt and Voth, 1998; Day et al., 2002) focus on the moving proton and its delocalization between different water molecules. Still, they cannot cope with a large number of different (de-)protonation sites in protic ionic liquids.

We present protex - an open-source Python-based tool for proton exchange in MD simulations. It works seamlessly with the OpenMM toolkit (Eastman et al., 2017) and can perform customized transfer reactions without restricting the number of (de-)protonation sites. Contrary to common Monte Carlo approaches (Baptista et al., 2002; Mongan et al., 2004) of accepting/denying the proton transfers, we perform the one-step proton hopping with a quantum-mechanically derived probability once a distance criterion between the hopping proton and the acceptor is met. However, our probabilities are determined by a Markov chain model [Jacobi et al. (2022)]. In contrast to Dreßler et al, (2020a) and Dreßler et al. (2020b) our Markov chain is applied prior to MD simulation to compute reasonable starting probabilities for the various reactions. Subsequently, the probabilities can be set manually to test several models for proton diffusion and to optimize the agreement with the experiment. The program package protex is freely available on GitHub (https://github.com/cbc-univie/protex).

2 Materials and methods

2.1 λ-states of protic ionic liquids in a single topology approach

The protic ionic liquid 1-methylimidazolium acetate [Im₁H][OAc] is in equilibrium with its neutral species 1-methylimidazole Im₁ and acetic acid HOAc as shown in Figure 1.

FIGURE 1

The program protex uses a single topology approach with two discrete λ-states to allow for the proton exchange. For imidazoles and acetate, we model the neutral species Im₁ and HOAc and the cation Im₁H⁺ and anion OAc⁻, respectively. In principle, it is also possible to model the protonated acetic acid (Ingenmey et al., 2018; Jacobi et al., 2022), which might be necessary for the Grotthus conductivity mechanism, but we restrict ourselves to the simple protonation scheme by Jacobi et al. (2022) for the sake of simplicity. The deprotonation of the Im₁H⁺ or HOAc is modeled by turning the hydrogen (HP) into a dummy atom (DM) which is part of the acetate OAc⁻ and imidazole Im₁ molecule.

In contrast to common alchemical mutations for proton transfer, the presented approach is not limited to partial charge mutations (Mey et al., 2020). As the atom types are changed to fit the DGenFF force field nomenclature (Chatterjee et al., 2019; Lin and MacKerell, 2019; Kumar et al., 2020) of the charged/neutral species, all bonded and non-bonded parameters are modified. Tables 1, 2 outline these changes in the atom types, Lennard-Jones parameters, partial charges q_iβ, and polarizabilities α_iβ. In imidazolium, both ring nitrogens share the same atom type NC. The neutral imidazole Im₁ has lone pairs at the unsubstituted ring nitrogen (NB). As a consequence of protonation, the charge of these lone pairs is set to 0 e, turning off all their interactions.

All these changes ensure that the molecules behave according to their charge state. This is particularly crucial for ionic liquids as the Coulombic interactions are neither short-ranged nor restricted to ion pairs and lead to cage-like structures (Schröder, 2011; Szabadi et al., 2022).

2.2 Polarizable force fields

During the protex update of the λ-states, significant changes in the atomic charges occur (see Table 2), which turn molecular ions into neutral molecules and vice versa. Since such drastic changes in electrostatic interactions between the molecules destabilize MD simulations, polarizable forces were applied to smoothen the transition of the Coulomb energy. These polarizable forces are anyway essential to close the gap between computational and experimental dynamics as non-polarizable force fields are usually one order of magnitude to viscous (Bedrov et al., 2019).

Among the different approaches to introduce polarizability to an MD simulation, we used the Drude oscillator model (Lamoureux and Roux, 2003; Bedrov et al., 2019). Each polarizable atom iβ is assigned an additional pair of Drude particles evoking an induced dipole . One Drude particle is located at the site of the atom itself with a charge of −q^δ. The second Drude particle carries the opposite charge + q^δ, which is generally negative as q^δ < 0e. The second Drude particle is attached by a harmonic spring to the first. The corresponding force constant of the spring is given by

And usually set to a constant value for all atoms in a simulation resulting in increasing Drude charges q^δ with increasing polarizabilities α_iβ. The induced dipole depends on the charge q^δ and the distance between the Drude particles: . It points from the Drude particle located at the polarizable atom to the mobile Drude particle, which can also be seen in Figure 2. The displacement should be much smaller than bond distances and is usually less than 0.1 Å. The total self-polarization energy reads .

FIGURE 2

As the Coulomb interaction of the Drude particles already contains the dispersion between molecules, the corresponding Lennard-Jones interactions have to be reduced to counteract double counting. Bypassing a complete reparametrization of the Lennard-Jones parameters, atomic -parameters can be scaled systematically as a function of the polarizability (Becker et al., 2016; Bedrov et al., 2019; Szabadi and Schröder, 2021; Joerg and Schröder, 2022):

Using the largest atomic polarizability α_max and the difference Δα between α_max and the polarizability α_iβ of the current polarizable atom. The scaling factor s determines the influence of the polarizability on the Lennard-Jones (Becker et al., 2016; Bedrov et al., 2019; Szabadi and Schröder, 2021; Joerg and Schröder, 2022).

2.3 The program package protex

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.

FIGURE 3

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), 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). However, as the corresponding calculations result in significant computational effort, we start with a fixed distance criterion. Dreßler et al, (2020a) and Dreßler et al. (2020b) 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)where and 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 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

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.

2.4 Computational setup

Details on the parametrization process of the molecular species involved can be found in Joerg and Schröder (2022). In short, the force field for Im₁H⁺OAc⁻ was based on the CHARMM General Force Field (CGenFF) (Kumar et al., 2020). 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 = 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), 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) 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). 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) 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) and Gong and Padua, (2021), 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).

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) but are updated applying Eq. 4. The time interval between the transfer checks was set to 10 ps

2.5 Analysis

For the analysis of the trajectories, the MDAnalysis package (Michaud-Agrawal et al., 2011; Gowers et al., 2016) was applied and augmented by self-written Python code. For example, the combination of MDAnalysis and the voro++ library (Rycroft, 2009) allows for the computation of the coordination number N_kl and the volume V_k(shell = 1) of the first solvation shell around molecules of species k (Zeindlhofer et al., 2018; Szabadi et al., 2022). Based on this information, a shell-based potential of mean forcecan be computed from the concentration c_l (shell = 1) = N_kl/V_k(shell = 1) of species l in the first shell around species k and the bulk concentration c_l = N_l/V. Negative PMF-values indicate preferential coordination of the species l around species k, whereas positive values result from a depletion of species l around k.

The diffusion coefficient can be calculated using the Einstein relation (Allen and Tildesley, 1986). For species k, it reads

With Δ r(t) = |r(t) − r (0)|. To obtain diffusion coefficients for each species in the system, the possible proton transfers, which consequently change the residue names, had to be accounted for. Therefore, the time series for each residue was cut when a transfer occurred. Only time series with at least 25 ns length were analyzed for the final analysis. Although this reduces the statistics of the mean-squared displacement, it ensures that the mobility of the charged and neutral compounds is not mixed. The slope for Eq. 6 was taken between 2 and 6 ns Additionally, diffusion coefficients for combined Im₁H⁺/Im₁ as well as OAc⁻/HOAc were calculated.

The analysis of the conductivity σ(0) needed some extra attention. Commonly, σ(0) is obtained from the mean-squared displacement

of the collective translational dipole moment using the molecular charges q_i and the respective center-of-masses from the unfolded trajectory. However, the occurrence of proton transfer reactions is decided on the minimum distance using the periodic boundary conditions during the production of the folded trajectory. In Figure 5, a proton transfer between 1-methylimidazole (turquoise dot) and 1-methylimidazolium (red dot) is sketched. The distance between the nitrogen of the Im₁ and the hydrogen of Im₁H⁺ is below 1.55 Å considering the periodic boundary conditions. If a proton transfer occurs, the imidazole charge is set to +1e, and the imidazolium becomes neutral. However, after unfolding the trajectory, the distance between the two exchange partners is much larger, and huge jumps in the collective dipole moment occur.

FIGURE 5

The simplest way to bypass this problem is to undo this huge jump in of the unfolded trajectory and subsequently add the contributions emerging from the occurred proton transfers. This way, one does not have to bookmark all toroidal jumps. The contribution of each proton transfer is

Using the center-of-masses and of the molecules i and j. The indices b and a denote before (b) and after (a) the proton transfer. Based on our simple reaction scheme in Table 4 eight different types of these contributions exist, which are tabulated in Table 5.

3 Results and discussion

In MD simulations, proton transfer events are usually harmful non-equilibrium situations as ions may become neutral or vice versa. Consequently, one expects significant jumps in the Coulomb energy of the system. This is undoubtedly true for non-polarizable MD simulations, but fortunately, the induced dipoles in polarizable trajectories counteract these jumps and smoothen the non-bonded (NB) interactions as shown in Figure 6 for s = 0.25 (blue) and s = 0.40 (orange). Strictly speaking, the non-bonded energy also comprises the Lennard-Jones interactions, but these do not change very much during the proton transfer as hydrogens usually have no significant contributions.

FIGURE 6

The time evolution of this non-bonded energy is a constant profile for s = 0.4 and gets more negative for s = 0.25. Interestingly, the Drude self-polarization energy, on the other hand, rises about the same amount in that case. Weakening of the Lennard-Jones spheres allows for closer distances of the induced dipoles of two polarizable atoms. Consequently, the interaction of these induced dipoles with other induced dipoles and with the permanent charges results in lower energy. Since this also leads to larger distances between the mobile Drude particle and the polarizable atom, the corresponding self-polarization term U^δδ increases. These effects in the non-bonded energy and self-polarization cancel out in the total energy, which changes roughly −7 kcal/mol/ns of the complete simulation box in case of s = 0.25 which might still be acceptable, although it is almost twice the drift per Drude oscillator compared to water [Lamoureux and Roux (2003)]. However, a scaling factor s of 0.40 is preferable as no drift is observed in Figure 6. Zooming into the trajectory (see inset in Figure 6), one notices that the jumps due to the multiple proton exchanges are less compared to the fluctuations of the non-bonded energy between two proton exchange events. This clearly demonstrates the induced dipoles’ functionality for stabilizing proton transfer MD simulations.

The simulation period of our polarizable MD simulations is 50 ns. As we stop the production every 10 ps to check for proton exchanges, (Jacobi et al., 2022), a molecule may face 5,000 exchanges at maximum. However, the average number of proton transfers for each molecule is much lower (see in the top panel of Figure 7) and equals roughly 10 to 15 exchanges on average.

FIGURE 7

Please note that the histograms are quite broad, revealing a heterogeneous system. Stronger Lennard-Jones-s-scaling leads to more transfers in general for all 500 Im₁/Im₁H⁺ and 500 HOAc/OAc⁻ [s = 0.25 (blue): 5,725 transfers; s = 0.4 (orange): 3,585 transfers] since the overall movement is increased.

The lower panel of Figure 7 depicts the relevance of the reactions in Table 4. Interestingly, the proton transfers are dominated by proton exchanges between imidazole and imidazolium, although the reaction probability is significantly lower than that of the reactions Im₁H⁺ + OAc⁻ or OAc⁻ + HOAc. This is due to the crucial distance between the hydrogen donor and acceptor, which was set to 1.55 Å in our simulations. As imidazole and imidazolium seem to come closer to each other and have the correct mutual orientation, these reactions happen more often than those with higher probability.

Our simulation still reproduces the equilibrium value of 30% ionic:70% neutral molecules (Jacobi et al., 2022; Joerg and Schröder, 2022) for both s-values as shown in Figure 8.

FIGURE 8

However, we had to apply Eq. 4 to prohibit drifting away from the equilibrium partitioning of the molecules as the total number of proton transfer events, i.e. 5,000, is much lower than in the Markov chain analysis reported by Jacobi et al. (2022). Furthermore, due to the distance criterion r_max and the mutual orientation of the reacting species in the liquid phase, particular reactions are favored regardless of the value of the reaction probability p_ref.

So far, we have shown that our polarizable MD simulations, including proton transfer, are stable for at least 50 ns with the correct ratio of charged and neutral molecules. However, the more interesting question is: What is the difference between a polarizable simulation with fixed molecular charges q_i and our new simulations, including proton transfers?

Table 6 shows the box length L, density ρ and conductivity σ(0) of the systems for scaling factors of s = 0.25 and s = 0.4. The box length and, thus, density are very similar for the different replicas, as well as compared to the simulations without proton transfer in our previous study (Joerg and Schröder, 2022). This is expected since the same workflow was used, and no proton transfers were allowed during the NpT runs, opposite to the NVT production run, which was used for analyzing transport properties. Hence, the conductivity is expected to differ. A notable increase was found for both systems allowing proton transfers as displayed in Table 6. Since conductivity is a collective property, the statistics are challenging explaining the slight deviations for the different replicas. Interestingly, the standard deviations in the case of s = 0.25 are significantly larger. Also, the conductivity σ(0) for s = 0.25 is above the experimental value, whereas σ(0) for s = 0.4 is within the range of the measured values.

Figure 9 depicts the diffusion coefficients of the four involved species for both s-scaling factors. The horizontal solid lines are the average diffusion coefficient of the imidazole-based and carboxylate-based molecules taking into account the different mole fractions and proton transfers. The dashed lines correspond to the polarizable simulations without proton transfers (Joerg and Schröder, 2022). Although the diffusion coefficients increased compared to Joerg and Schröder (2022) (black arrows), they are still smaller than the corresponding experimental value (gray star) for both species (Thawarkar et al., 2019). As expected, the diffusion coefficients of the neutral molecules are higher than their charged counterpart because of fewer Coulombic interactions. This was also true for the polarizable simulations without proton transfers (black x in Figure 9) (Joerg and Schröder, 2022). Except for imidazole, the diffusion coefficients of the species are more or less unaffected by the implementation of the proton transfers. Overall, the molecular translational motion characterized by the diffusion coefficients is not responsible for the increase in the conductivity, which has to be due to collective effects.

FIGURE 9

Cage effects can be characterized by the shell-resolved potential of mean force PMF. Since we are interested in the different behavior of polarizable simulations with and without proton transfer, we plotted the differences ΔPMF of the mutual shell-resolved potential of mean forces PMF_kl (shell = 1) for the species k, l ∈ {Im₁H⁺, OAc⁻, Im₁, HOAc} as a heat map in Figure 10.

FIGURE 10

Red boxes indicate that the solvation became less favorable after including proton transfers, whereas blue boxes reveal an increased attraction compared to the simulations without proton transfer. For the conductivity, the ΔPMFs, including the charged species in the top left regions of the heat map, are the most interesting. The dark red boxes for Im₁H⁺/OAc⁻ and Im₁H⁺/Im₁H⁺ correspond to weaker coordination of these species and fewer cation/anion or cation/cation pairs increase the conductivity σ(0). The overall charge of cation/anion pairs is zero; consequently, this pair does not contribute to σ(0). If two cations stick together for a long time, their overall mobility is reduced, hence the electric current. Allowing for proton transfer in the polarizable simulations has multiple effects: First, Im₁H⁺/OAc⁻ may react and become two neutral molecules. This reaction does not increase the conductivity. Second, in the case of a cation/cation pair, one of the imidazoliums may exchange its proton with acetate. Now, the second imidazolium has a new imidazole and acetic acid neighbor and may be more mobile than in the cation/cation aggregate before. This would increase the conductivity. Quite generally, in Figure 10, neutral molecules seem to accumulate in the immediate neighborhood of charged molecules (blue boxes in the top right region of the heat map) as a consequence of the multiple proton transfer events. This fact shows the weakening of ion cages by proton transfer reactions.

The proton transfers promote the diffusion of imidazoles Im₁. If an imidazolium inside an ion cage transfers its acidic proton to one of the acetates, the emerging Im₁ still encounters many other acetates in the former ion cage. This situation is energetically unfavorable, and the imidazole will try to escape immediately, thereby increasing the diffusion coefficient. The corresponding ΔPMF is −0.16 kJ·mol⁻¹ (see Figure 10).

4 Conclusion

The lightweight open-source Python package protex was successfully implemented for polarizable MD simulations of the protic ionic liquid 1-methylimidazolium acetate. In contrast to constant pH simulation techniques handling proteins’ (de-)protonation, the current work deals with proton transfer within the solvent. Protex augments an OpenMM simulation object and is, therefore, straightforwardly usable with any polarizable OpenMM simulation and not restricted to protic ionic liquids. The transfer can either be instantaneously or through intermediate λ-states, with user-defined distances and probability criteria.

Allowing for proton transfers overcomes one of the critical limitations in classical MD simulations, i.e., the formation and breaking of bonds. However, proton transfers are essential for the meaningful simulation of protic ionic liquids or other proton-exchanging solvents. In the case of the protic ionic liquid 1-methylimidazolium acetate, a slight increase in the diffusion coefficients of all species is accompanied by a significant increase in the overall conductivity σ(0) of the system, which is now in excellent agreement with the experimental values.

References

• 1

  Allen M. P. Tildesley D. J. (1986). Computer simulation of liquids. OxfordClarndon press.

  • Google Scholar

• 2

  Baptista A. M. Teixeira V. H. Soares C. M. (2002). Constant pH molecular dynamics using stochastic titration. J. Chem. Phys.117, 4184–4200. 10.1063/1.1497164

  • CrossRef
  • Google Scholar

• 3

  Becker T. M. Dubbeldam D. Lin L.-C. Vlugt T. J. (2016). Investigating polarization effects of CO2 adsorption in MgMOF-74. J. Comput. Sci.15, 86–94. International Computational Science and Engineering Conference 2015 (ICSEC15). 10.1016/j.jocs.2015.08.010

  • CrossRef
  • Google Scholar

• 4

  Bedrov D. Piquemal J.-P. Borodin, MacKerell Jr. A. D. Roux B. Schröder C. (2019). Molecular dynamics simulations of ionic liquids and electrolytes using polarizable force fields. Chem. Rev.119, 7940–7995. 10.1021/acs.chemrev.8b00763

  • CrossRef
  • Google Scholar

• 5

  Boresch S. Karplus M. (1999b). The role of bonded terms in free energy simulations. 2. Calculation of their influence on free energy differences of solvation. J. Phys. Chem. A103, 119–136. 10.1021/jp981629f

  • CrossRef
  • Google Scholar

• 6

  Boresch S. Karplus M. (1999a). The role of bonded terms in free energy simulations: 1. Theoretical analysis. J. Phys. Chem. A103, 103–118. 10.1021/jp981628n

  • CrossRef
  • Google Scholar

• 7

  Börjesson U. Hünenberger P. H. (2001). Explicit-solvent molecular dynamics simulation at constant pH: Methodology and application to small amines. J. Chem. Phys.114, 9706–9719. 10.1063/1.1370959

  • CrossRef
  • Google Scholar

• 8

  Brooks B. R. Brooks C. L. III Mackerell A. D. Jr. Nilsson L. Petrella R. J. Roux B. et al (2009). Charmm: The biomolecular simulation program. J. Comput. Chem.30, 1545–1614. 10.1002/jcc.21287

  • CrossRef
  • Google Scholar

• 9

  Chatterjee P. Heid E. Schröder C. MacKerell A. (2019). Polarizable general force field for drug-like molecules Drude general force field (DGenFF). Biophys. J.116, 142a. 10.1016/j.bpj.2018.11.788

  • CrossRef
  • Google Scholar

• 10

  Chen J. Chen L. Lu Y. Xu Y. (2014). Physicochemical properties of aqueous solution of 1-methylimidazolium acetate ionic liquid at several temperatures. J. Mol. Liq.197, 374–380. 10.1016/j.molliq.2014.05.027

  • CrossRef
  • Google Scholar

• 11

  Day T. J. F. Soudackov A. V. Čuma M. Schmitt U. W. Voth G. A. (2002). A second generation multistate empirical valence bond model for proton transport in aqueous systems. J. Chem. Phys.117, 5839–5849. 10.1063/1.1497157

  • CrossRef
  • Google Scholar

• 12

  Dobrev P. Vemulapalli S. P. B. Nath N. Griesinger C. Grubmüller H. (2020). Probing the accuracy of explicit solvent constant pH molecular dynamics simulations for peptides. J. Chem. Theory Comput.16, 2561–2569. 10.1021/acs.jctc.9b01232

  • CrossRef
  • Google Scholar

• 13

  Donnini S. Ullmann R. T. Groenhof G. Grubmüller H. (2016). Charge-neutral constant pH molecular dynamics simulations using a parsimonious proton buffer. J. Chem. Theory Comput.12, 1040–1051. 10.1021/acs.jctc.5b01160

  • CrossRef
  • Google Scholar

• 14

  Dreßler C. Kabbe G. Brehm M. Sebastiani D. (2020a). Dynamical matrix propagator scheme for large-scale proton dynamics simulations. J. Chem. Phys.152, 114114. 10.1063/1.5140635

  • CrossRef
  • Google Scholar

• 15

  Dreßler C. Kabbe G. Brehm M. Sebastiani D. (2020b). Exploring non-equilibrium molecular dynamics of mobile protons in the solid acid csh2po4 at the micrometer and microsecond scale. J. Chem. Phys.152, 164110. 10.1063/5.0002167

  • CrossRef
  • Google Scholar

• 16

  Eastman P. Swails J. Chodera J. D. McGibbon R. T. Zhao Y. Beauchamp K. A. et al (2017). OpenMM 7: Rapid development of high performance algorithms for molecular dynamics. PLoS Comput. Biol.13, e1005659. 10.1371/journal.pcbi.1005659

  • CrossRef
  • Google Scholar

• 17

  Gong Z. Padua A. A. H. (2021). Effect of side chain modifications in imidazolium ionic liquids on the properties of the electrical double layer at a molybdenum disulfide electrode. J. Chem. Phys.154, 084504. 10.1063/5.0040172

  • CrossRef
  • Google Scholar

• 18

  Goossens K. De Winter H. (2018). Molecular dynamics simulations of membrane proteins: An overview. J. Chem. Inf. Model.58, 2193–2202. 10.1021/acs.jcim.8b00639

  • CrossRef
  • Google Scholar

• 19

  Gowers R. J. Linke M. Barnoud J. Reddy T. J. E. Melo M. N. Seyler S. L. et al (2016). “MDAnalysis: A Python package for the rapid analysis of molecular dynamics simulations,” in Proceedings of the 15th Python in Science Conference–105. Sebastian Benthall and Scott Rostrup.98

  • Google Scholar

• 20

  Hospital A. Goñi J. R. Orozco M. Gelpí J. L. (2015). Molecular dynamics simulations: Advances and applications. Adv. Appl. Bioinform. Chem.8, 37–47. 10.2147/aabc.s70333

  • CrossRef
  • Google Scholar

• 21

  Hou H.-Y. Huang Y.-R. Wang S.-Z. Bai B.-F. (2011). Preparation and physicochemical properties of imidazolium acetates and the conductivities of their aqueous and ethanol solutions. Acta Phys. Chim. Sin.27, 2512–2520. 10.3866/pku.whxb20111120

  • CrossRef
  • Google Scholar

• 22

  Ingenmey J. Gehrke S. Kirchner B. (2018). How to harvest Grotthuss diffusion in protic ionic liquid electrolyte systems. ChemSusChem11, 1900–1910. 10.1002/cssc.201800436

  • CrossRef
  • Google Scholar

• 23

  Jacobi R. Joerg F. Steinhauser O. Schröder C. (2022). Emulating proton transfer reactions in the pseudo-protic ionic liquid 1-methylimidazolium acetate. Phys. Chem. Chem. Phys.24, 9277–9285. 10.1039/d2cp00643j

  • CrossRef
  • Google Scholar

• 24

  Joerg F. Schröder C. (2022). Polarizable molecular dynamics simulations on the conductivity of pure 1-methylimidazolium acetate systems. Phys. Chem. Chem. Phys.24, 15245–15254. 10.1039/d2cp01501c

  • CrossRef
  • Google Scholar

• 25

  Khandogin J. Brooks C. L. III (2005). Constant pH molecular dynamics with proton tautomerism. Biophys. J.89, 141–157. 10.1529/biophysj.105.061341

  • CrossRef
  • Google Scholar

• 26

  Kumar A. Yoluk O. Collectivityrell A. D. Jr. (2020). FFParam Standalone package for CHARMM additive and Drude polarizable force field parametrization of small molecules. J. Comput. Chem.41, 958–970. 10.1002/jcc.26138

  • CrossRef
  • Google Scholar

• 27

  Lamoureux G. Roux B. (2003). Modeling induced polarization with classical Drude oscillators: Theory and molecular dynamics simulation algorithm. J. Chem. Phys.119, 3025–3039. 10.1063/1.1589749

  • CrossRef
  • Google Scholar

• 28

  Lee M. S. Salsbury F. R. Jr Brooks C. L. III (2004). Constant-pH molecular dynamics using continuous titration coordinates. Proteins Struct. Funct. Genet.56, 738–752. 10.1002/prot.20128

  • CrossRef
  • Google Scholar

• 29

  Lill M. A. Helms V. (2001). Compact parameter set for fast estimation of proton transfer rates. J. Chem. Phys.114, 1125–1132. 10.1063/1.1332993

  • CrossRef
  • Google Scholar

• 30

  Lin F.-Y. MacKerell A. D. (2019). “Force fields for small molecules,” in Biomolecular simulations: Methods and protocols. Editors BonomiM.CamilloniC. (New York, NY: Springer New York), 21–54.

  • Google Scholar

• 31

  MacFarlane D. R. Pringle J. M. Johansson K. M. Forsyth S. A. Forsyth M. (2006). Lewis base ionic liquids. Chem. Commun., 1905–1917. 10.1039/b516961p

  • CrossRef
  • Google Scholar

• 32

  Martínez L. Andrade R. Birgin E. G. Martínez J. M. (2009). Packmol: A package for building initial configurations for molecular dynamics simulations. J. Comput. Chem.30, 2157–2164. 10.1002/jcc.21224

  • CrossRef
  • Google Scholar

• 33

  Martyna G. J. Klein M. L. Tuckerman M. (1992). Nosé–hoover chains: The canonical ensemble via continuous dynamics. J. Chem. Phys.97, 2635–2643. 10.1063/1.463940

  • CrossRef
  • Google Scholar

• 34

  Mey A. S. J. S. Allen B. K. Bruce Macdonald H. E. Chodera J. D. Hahn D. F. Kuhn M. et al (2020). Best practices for alchemical free energy calculations [article v1.0]. Living J. Comp. Mol. Sci.2, 18378. 10.33011/livecoms.2.1.18378

  • CrossRef
  • Google Scholar

• 35

  Michaud-Agrawal N. Denning E. J. Woolf T. B. Beckstein O. (2011). MDAnalysis: A toolkit for the analysis of molecular dynamics simulations. J. Comput. Chem.32, 2319–2327. 10.1002/jcc.21787

  • CrossRef
  • Google Scholar

• 36

  Mongan J. Case D. A. (2005). Biomolecular simulations at constant pH. Curr. Opin. Struct. Biol.15, 157–163. 10.1016/j.sbi.2005.02.002

  • CrossRef
  • Google Scholar

• 37

  Mongan J. Case D. A. McCammon J. A. (2004). Constant pH molecular dynamics in generalized Born implicit solvent. J. Comput. Chem.25, 2038–2048. 10.1002/jcc.20139

  • CrossRef
  • Google Scholar

• 38

  Radak B. K. Chipot C. Suh D. Jo S. Jiang W. Phillips J. C. et al (2017). Constant-pH molecular dynamics simulations for large biomolecular systems. J. Chem. Theory Comput.13, 5933–5944. 10.1021/acs.jctc.7b00875

  • CrossRef
  • Google Scholar

• 39

  Russo M. F. Jr Van Duin A. C. (2011). Atomistic-scale simulations of chemical reactions: Bridging from quantum chemistry to engineering. Nucl. Instrum. Methods Phys. Res. B Beam Interact. Mat. At.269, 1549–1554. 10.1016/j.nimb.2010.12.053

  • CrossRef
  • Google Scholar

• 40

  Rycroft C. H. (2009). Voro++: A three-dimensional voronoi cell library in C++. Chaos19, 041111. 10.1063/1.3215722

  • CrossRef
  • Google Scholar

• 41

  Schmitt U. W. Voth G. A. (1998). Multistate empirical valence bond model for proton transport in water. J. Phys. Chem. B102, 5547–5551. 10.1021/jp9818131

  • CrossRef
  • Google Scholar

• 42

  Schröder C. (2011). Collective translational motions and cage relaxations in molecular ionic liquids. J. Chem. Phys.135, 024502. 10.1063/1.3601750

  • CrossRef
  • Google Scholar

• 43

  Schröder C. Sonnleitner T. Buchner R. Steinhauser O. (2011). The influence of polarizability on the dielectric spectrum of the ionic liquid 1-ethyl-3-methylimidazolium triflate. Phys. Chem. Chem. Phys.13, 12240–12248. 10.1039/c1cp20559e

  • CrossRef
  • Google Scholar

• 44

  Schröder C. Steinhauser O. (2010). Simulating polarizable molecular ionic liquids with Drude oscillators. J. Chem. Phys.133, 154511. 10.1063/1.3493689

  • CrossRef
  • Google Scholar

• 45

  Shirts M. R. (2012). “Best practices in free energy CalDculations for drug design,” in Computational drug discovery and design. Editor BaronR. (New York, NY: Springer New York), 425–467.

  • Google Scholar

• 46

  Son C. Y. McDaniel J. G. Cui Q. Yethiraj A. (2019). Proper thermal equilibration of simulations with Drude polarizable models: Temperature-grouped dual-Nosé–Hoover thermostat. J. Phys. Chem. Lett.10, 7523–7530. 10.1021/acs.jpclett.9b02983

  • CrossRef
  • Google Scholar

• 47

  Szabadi A. Honegger P. Schöfbeck F. Sappl M. Heid E. Steinhauser O. et al (2022). Collectivity in ionic liquids: A temperature dependent, polarizable molecular dynamics study. Phys. Chem. Chem. Phys.24, 15776–15790. 10.1039/d2cp00898j

  • CrossRef
  • Google Scholar

• 48

  Szabadi A. Schröder C. (2021). Recent developments in polarizable molecular dynamics simulations of electrolyte solutions. J. Comput. Biophys. Chem.21, 415–429. 10.1142/s2737416521420035

  • CrossRef
  • Google Scholar

• 49

  Thawarkar S. Khupse N. D. Shinde D. R. Kumar A. (2019). Understanding the behavior of mixtures of protic-aprotic and protic-protic ionic liquids: Conductivity, viscosity, diffusion coefficient and ionicity. J. Mol. Liq.276, 986–994. 10.1016/j.molliq.2018.12.024

  • CrossRef
  • Google Scholar

• 50

  Weismiller M. R. Junkermeier C. E. Russo M. F. Salazar M. R. Bedrov D. Van Duin A. C. (2015). ReaxFF molecular dynamics simulations of intermediate species in dicyanamide anion and nitric acid hypergolic combustion. Model. Simul. Mat. Sci. Eng.23, 074007. 10.1088/0965-0393/23/7/074007

  • CrossRef
  • Google Scholar

• 51

  Zeindlhofer V. Berger M. Steinhauser O. Schröder C. (2018). A shell-resolved analysis of preferential solvation of coffee ingredients in aqueous mixtures of the ionic liquid 1-ethyl-3-methylimidazolium acetate. J. Chem. Phys.148, 193819. 10.1063/1.5009802

  • CrossRef
  • Google Scholar

• 52

  Zhang B. van Duin A. C. T. Johnson J. K. (2014). Development of a ReaxFF reactive force field for tetrabutylphosphonium glycinate/CO2 mixtures. J. Phys. Chem. B118, 12008–12016. 10.1021/jp5054277

  • CrossRef
  • Google Scholar