Protex: Proton Transfer Simulation in OpenMM

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.