Simulating Metal-Imidazole Complexes

One commonly observed binding motif in metalloproteins involves the interaction between a metal ion and histidine’s imidazole side chains. Although previous imidazole-M(II) parameters established the flexibility and reliability of the 12–6–4 Lennard-Jones (LJ)-type nonbonded model by simply tuning the ligating atom’s polarizability, they have not been applied to multiple-imidazole complexes. To fill this gap, we systematically simulate multiple-imidazole complexes (ranging from one to six) for five metal ions (Co(II), Cu(II), Mn(II), Ni(II), and Zn(II)) which commonly appear in metalloproteins. Using extensive (40 ns per PMF window) sampling to assemble free energy association profiles (using OPC water and standard HID imidazole charge models from AMBER) and comparing the equilibrium distances to DFT calculations, a new set of parameters was developed to focus on energetic and geometric features of multiple-imidazole complexes. The obtained free energy profiles agree with the experimental binding free energy and DFT calculated distances. To validate our model, we show that we can close the thermodynamic cycle for metal-imidazole complexes with up to six imidazole molecules in the first solvation shell. Given the success in closing the thermodynamic cycles, we then used the same extended sampling method for six other metal ions (Ag(I), Ca(II), Cd(II), Cu(I), Fe(II), and Mg(II)) to obtain new parameters. Since these new parameters can reproduce the one-imidazole geometry and energy accurately, we hypothesize that they will reasonably predict the binding free energy of higher-level coordination numbers. Hence, we did not extend the analysis of these ions up to six imidazole complexes. Overall, the results shed light on metal–protein interactions by emphasizing the importance of ligand–ligand interaction and metal-π-stacking within metalloproteins.


INTRODUCTION
−5 The absence of specific metal ions can cause fatal deficiencies, such as carcinogenesis, severe malnutrition, and eventually death.Over 25% of proteins contain metal ions that can function in either a structural or catalytic role 6,7 and are targets for the design of new pharmaceutical agents. 8Computational chemistry has become an effective tool for examining metal ion-coordinating systems in various biological systems, like proteins, nucleic acids, carbohydrates, and lipids, to understand biological structure and function. 9−13 Among different computational techniques, force field models and methods afford the computational speed needed to study complex biomacromolecules.However, reproducing metal ions' structural features and thermodynamic properties in water or protein systems is still challenging for these methods under the prerequisite of maintaining a low computational cost 14 and being physically meaningful.Apart from the accurate machine learning models, 15−18 some commonly used physically meaningful force field models 19 include but are not limited to bonded models, nonbonded models, 20 Drude oscillators models, 21,22 cationic dummy atom (CDA) and CDA pol models, 23,24 and the ReaxFF model. 25n the bonded model, covalent bonds exist between the metal ion and the coordinated residues, in which the bond angle, dihedral, van der Waals, and electrostatic interactions are defined by classical terms. 26Although the bonded model can successfully replicate experimentally determined structures, it cannot simulate the change of coordination number or ligating residues in order to model catalytic metal centers or metal ion transport. 27he nonbonded model, in comparison, is another widely used model for metal ions, where the metal ion is represented by a sphere using both van der Waals (vdW) and Coulombic terms to interact with its surroundings. 19The vdW interactions are defined by the 12−6 Lennard-Jones (LJ) 28 or Born-Mayer potential. 29−35 The deficiency of the 12−6 model is primarily because it does not include ion-induced dipole interactions, which is an essential feature in highly polarized systems.The 12−6−4 LJ model was developed to overcome this challenge by adding the C 4 term to account for the ion-induced dipole interaction.The ion-induced dipole interaction is proportional to r −4 , where r is the distance between the two particles. 36−27,30−43 Since the 12−6−4 model is easy to apply, computationally efficient, and accurately describes the interactions between the metal ion and its coordinating ligands, it is an excellent model for simulating metal ions in molecular dynamics (MD) simulations.In earlier work, 44,45 the 12−6−4 Lennard-Jones (LJ)-type nonbonded model has successfully simulated metal ion systems using the Potential of Mean Force (PMF) method using a modified polarizability of the metal-chelating nitrogen. 44This work parametrized the HID (δ-nitrogen protonated) imidazole molecules against the experimental values for 11 metals (Ag(II), Ca(II), Cd(II), Co(II), Cu(II), Cu(I), Fe(II), Mg(II), Mn(II), Ni(II), and Zn(II)) in conjunction with the commonly used OPC water model. 46To explore the capability of these new parameters, we selected five ions (Co(II), Cu(II), Mn(II), Ni(II), and Zn(II)) with experimental values available 47 for a detailed analysis of multiple-imidazole interactions from the list of 11 ions.By computing these five ions' absolute binding free energies to imidazole(s) (ranging from one to six) and comparing the results with the experiment, this parametrization routine of using PMFs with extended simulation (40 ns per window) was shown to be reliable at producing transferable parameters.These new parameters can then serve as good priors to predicting the metal-multiple-imidazole interactions that currently lack experimental data.
To explore the new parameter's ability to accurately model the addition of multiple imidazoles to a metal center, we used thermodynamic cycles to demonstrate that we can model multiple-imidazole complexes up to six ligands.Finally, due to the lack of metal−ligand distance information, DFT calculations were conducted to obtain estimates for the metal-imidazole distances in multiple-imidazole complexes.The derived 12−6− 4 models are compatible with the AMBER class of force fields and with the recommended OPC water model to allow a range of applications.We hope that the derived parameters will be useful in studying metalloprotein structure and function and for transition metal ion transport.

Optimization of the
where e depicts the proton charge while Q i and Q j represent the partial charge of atoms i and j, where i and j are the indices of the ion and one of the ligands, respectively.The Coulomb pair potential was utilized to represent the electrostatic interaction between atoms i and j.In contrast, the classic 12−6 LJ potential plus an extra r 4 term represented the van der Waals interactions.The C 4 terms between water and ions were taken from our previous studies. 20,36,40The C 4 terms between histidine and metal ions were optimized based on the following equation: where α 0 is an atom type-dependent polarizability, which makes the C 4 different between ligands while the metal binds to water molecules and imidazole (which mimics the side chain of histidine) simultaneously.θ represents the angle between the ion-ligand line and the induced dipole direction, which, in our case, is always zero.PMF calculations were used to optimize the pairwise parameters to reproduce the experimental free energies of each metal bound to imidazole.The HID charge used on imidazole molecules is described in Figure 1, which is the same as what ff19SB uses in AMBER20AmberTools21. 48The major features in the HID-charged histidine are the nitrogen charge values and the protonation locations.The connecting carbon (CC in Figure 1A) is not connected to histidine's α-carbon but to hydrogen in the present work. 45The related RESP 49 charge fitted mol2 file and frcmod file of this HID-charged imidazole molecule can be found at https://github.com/lizhen62017/MIRIAM/tree/main/TI/ff19sb_CC.frcmod and https:// github.com/lizhen62017/MIRIAM/tree/main/TI/cooresp.mol2.

Parametrization Processes.
A PMF simulation protocol was prepared according to the MD procedure described in Section 2.3.The parametrization was conducted iteratively, as Figure 1B indicates.For each iteration, a PMF with an umbrella sampling of either 4 or 40 ns for each window was performed to calculate the binding free energies.When the PMF-calculated interaction energy is within ±0.1 kcal/mol of the experimental binding free energy, the polarizability, associated C 4 value, and simulated binding free energy will all be recorded.If the umbrella sampling does not give a binding free energy within the range of ±0.1 kcal/mol of the experimental value, a more precise estimate of the parameter can be obtained using linear extrapolation, assuming that the parameter and interaction energy follow a linear relationship.Then, a new round of PMF simulations is initiated with a newly assigned polarizability value to continue the iteration until an accurate binding free energy is obtained.Similarly, for the PLUMED runs, the distance restriction was instead performed by PLUMED during the umbrella sampling step of the PMF. 50he same 4 or 40 ns simulation time was used to construct a solid comparison between AMBER results.

Molecular Dynamics Simulations Using PMF Simulations.
The CUDA version of PMEMD from the AMBER20 48 package performed all the molecular dynamics (MD) simulations.The system was prepared using the tLEaP program by dissolving one ion with the designated numbers of imidazole in a 24 Å diameter octahedral box, which consists of about 5000 water molecules and corresponds to an ion concentration of approximately 10 mM.The minimization of the system was performed in three steps: (a) minimization of water molecules, with the imidazole group and ions restrained; (b) minimization of side chain hydrogens (which are imidazole hydrogens) on the nonwater molecules; (c) minimization of the whole system.Each step consists of 10,000 cycles of minimization using the steepest descent method followed by 10,000 cycles using the conjugate gradient method.In the next step, the system was heated to 300 K gradually during a 1 ns NVT simulation.A 550,000-step run was done at 300 K for system equilibration, employing the NPT ensemble.Then, a production run for 500,000 steps at 300 K under NPT conditions was performed to prepare and fully equilibrate the system.Later, a Jarzinski NMR restriction was applied to carry out a steered MD simulation between the metal ion and εnitrogen of one of the imidazole molecule(s) to pull them away.
After the steered MD is completed, about 120 snapshots will be taken from the steered MD trajectory to start 2 ns of NPT equilibrium and a 40 ns production run with umbrella sampling as implemented in both AMBER and PLUMED.Finally, the umbrella sampling results are analyzed using the WHAM (Weighted Histogram Analysis Method) program 51 to eventually get the PMF of one imidazole leaving the metal-imidazole complex.The process was automated in a workflow called MIRIAM (Metal-Imidazole Rational Interaction Analysis Method) and can be acquired via this link: https://github.com/lizhen62017/MIRIAM/PMF.
The Langevin thermostat with a collision frequency of 1 ps −1 was applied to control the temperature, and the Berendsen barostat, 52 with a pressure relaxation time of 5 ps, was employed for the pressure control.The time step was 2 fs, and the nonbonded cutoff was 10 Å.The SHAKE 53 algorithm was used to constrain bonds involving hydrogen atoms, and the time step was set to 2 fs.Cluster analysis was utilized to obtain the most representative structure from the MD simulations.The UCSF Chimera 34 and VMD programs were used to visualize and prepare the figures.

Molecular Dynamics Simulations Using Thermodynamic Integration (TI).
To validate our PMF results, we utilized TI calculations to form a thermodynamic cycle that is supposed to close (more details below).The TI simulation starts with building a system with the two desired metal ions overlapped in the initial structure.Then the two ions were masked as "timask1" and "timask2" respectively.After this, the system will undergo the same minimization, heating, and equilibration as the MD simulation for PMF does.Then, a 12window TI was applied to the system with up to six imidazole ligands and the metalloprotein 54 to slowly mutate the "timask1" ion to "timask2" and compare the relative energy change.The process-related files can also be found at https://github.com/lizhen62017/MIRIAM/TI.

DFT Calculation on the Metal-imidazole Distances and Generation of the C 4 Parameters.
Because PMF studies can provide energetic and geometric information, we wanted to evaluate the new parameter's performance in predicting metalimidazole distances and not just binding free energies.However, since the experimental metal-imidazole distances are unknown, we used DFT calculations to estimate these quantities.Induction energies were also calculated using the Symmetryadapted perturbation theory (SAPT) approach. 55−57 Both DFT  2, where extended sampling reveals the limit of one-imidazole parameters (red), which was previously not discovered by regular sampling (blue).This issue was resolved by introducing multipleimidazole parameters (green).Color matches the color scheme for the column headers in Table 2.
and SAPT are QM methods, with the latter being able to estimate the C 4 parameters based on induction energies so that the ion-induced dipole interaction can be calculated.
For the 1−4 imidazole systems, the initial structures were created using the graphical interface of GaussView 6.For these systems, we kept the relative positions of the imidazole molecules similar to that of the lowest energy structure found in the PMF study.DFT calculations were carried out using Becke's 58 three-parameter functional and the correlation function of Lee et al. (B3LYP) 59 and range-separated hybrid functional wB97XD. 60The geometry optimization of all these complexes was done using a triple-ζ basis set 6-311G++(2d,2p) with diffuse functions.Tight optimization criteria were employed with ultrafine grid integration methods.Vibrational frequency analysis was carried out to find the true minima of the complexes.
Gas phase geometries may not be able to reproduce metal ionligand distances in aqueous solution.So, we used the universal solvation model SMD 61 to model an aqueous solution.−73 Therefore, additional explicit water molecules were added to the first solvation shell to saturate the coordination of the metal ions of interest.In our calculations, we considered that these metals form high-spin, six coordinated complexes in the aqueous solution.To support our hypothesis, we calculated the energies from the optimized geometries of the metal ions that can form both high-spin and low-spin complexes (Mn(II), Ni(II), Co(II)).The high-spin systems were energetically more favorable, so only the high-spin results are presented.
The SAPT method draws a bridge between QM operators and the 12−6−4 parameters.Specifically, SAPT decomposes the molecular interaction into four different terms: electrostatic, induction, repulsion, and dispersion.Among these terms, the induction term can be interpreted as the charge-induced dipole interaction and shows a distance dependence of 1/r 4 in the longrange, which naturally fits the C 4 term in the 12−6−4 model. 57n the present study, we calibrated the C 4 parameters for metalimidazole interactions in different complexes based on the SAPT calculated induction energy either at the equilibrium distance or along a scan of the ion-imidazole distance from 1.7 to 3.0 Å.Specifically, the induction energy can be used to derive the C 4 parameters in eq 2 based on the equation C 4 = −E ind r 4 .

Overview of Target Binding Free Energies, Ion
Properties, and Parametrization Results.The 12−6−4 parameter set was used to reproduce the experimental metalimidazole binding free energies compared to the experimental results.The target binding free energies computed from experimental logK values are shown in Table 1. 47Among all the experimental values, Mn(II), Ni(II), and Zn(II) were obtained from thermometric titration, 74−77 Co(II) was obtained from potentiometric titration with perchloric acid, 78 while Cu(II) was obtained from UV−vis spectra during ligand exchange of the one-ion-three-imidazole system. 79The experimental free energies of each ion show a high dependence on both the ionic radius and the electronic configuration.Overall, the larger effective ionic radius ions have more negative binding free energies with imidazole, with a few exceptions in the transition metal series.For these exceptions, when the d-orbitals of the ions are full or half full, they tend to have higher, meaning less negative binding free energies with imidazole, indicating a connection between d-orbital symmetry, classic Crystal-Field Stabilization Energies (CFSE), and imidazole interaction energies. 47fter investigating the correlation between ion properties and their binding free energies with imidazole molecules, the parametrization process followed the protocol summarized in the method section.As noted before, the default parameters significantly underbind to imidazole, which led us to create a set of imidazole-M(II) parameters (see Table 2). 44As we began to build our multiple-imidazole parameter set, we discovered that the initial protocol of 4 ns per window missed a pi-stacked intermediate (see below), which was observed only after 40 ns per window.Missing the pi-stack intermediate led to a slight overestimation of the binding affinity for one imidazole binding to an M(II) ion.However, if multiple ligands are coordinating with the metal ions, this overestimation will accumulate and lead to more significant deviations (see the one-imidazole-compatible parameter column in Table 2). 44Therefore, lower polarizability values for nitrogen are needed.After a reparameterization process using extended sampling, the new parameters were found to be better converged and give better predictions for multiple-imidazole systems.Although the previous parameters 44 continue to be effective while being applied to single-imidazole systems, we recommend the newer parameter set to be used when modeling systems involve multiple imidazoles.The results of six replicate runs for each parameter set (three AMBER and three PLUMED) are summarized in Table 2.
With the exception of Ag(I), the majority of parametrized C 4 values have been found to decrease slightly from the ones used in the previous one-imidazole study. 44After decreasing, these newly parametrized C 4 values give binding free energy values closer to the experimental ones using the extended sampling protocol in both AMBER and PLUMED.By analyzing the PMF profile as illustrated in Figure 3, we found that longer simulations allow the formation of a tilted "cation-pi-stacking" conformation when the imidazole molecule is located near the second solvation shell of the metal ion. 81From this figure, we observe the generation of cation-pi-like interactions, while the previous runs using 4 ns did not observe this state.With the inclusion of this cation-pi-stacking state, the imidazole will have a higher energy barrier to move from the first solvation shell to the second due to the hydrogen bonds between water-imidazole.Interestingly, the imidazole-imidazole hydrogen bond does not

Journal of Chemical Theory and Computation
contribute significantly to this energy shift, as the black/gray PMFs in Figure S3 indicate.This is likely due to the weak, but not negligible, deshielding effect of ε-Carbon, which contributes to a C−H NMR shift of 7.7 ppm in imidazole, where a normal C−H NMR shift should be around 7.3 ppm. 82In summary, this energy barrier discovered by extended sampling caused the previous parameters to overestimate the binding free energy, which led us to create the new parameter sets shown in the left column of Table 2.

Using the Parameters for Thermodynamic Cycle Analysis.
One way to further test the derived parameters is through the use of thermodynamic cycles.Using this approach provides another check on the quality of the derived parameters.In particular, we use the thermodynamic cycle shown in Figure 4 (Ni(II) → Co(II) example here) where we convert Ni(II) → Co(II) in aqueous solution and the presence of 1 to 6 imidazole molecules.Finally, the sum of the free energies around the thermodynamic cycle should be as close to 0 as possible (since free energy is a state function).With the help of this new set of parameters generated by the extended sampling method, we can better close the thermodynamic cycle.As Table 3B,C show, the previous one-imidazole-compatible and default parameters do not converge (MAE of 2.30 and 2.75 kcal/mol, which is (1.40 + 2.53 + 3.31 + 2.62 + 0.11 + 3.72)/6 and (2.05 + 3.49 + 4.45 + 4.45 + 1.76 + 0.32)/6, respectively, as the yellow column indicates), with the deviation from 0 usually being well above the uncertainty estimated from the triplicate runs.In contrast, Table 3A shows that the new parameters can generate smaller net thermodynamic cycle values (MAE of 0.61 kcal/mol, which is (0.33 + 0.18 + 0.03 + 0.37 + 1.11 + 1.64)/6 as the yellow column indicates), and the deviation from 0 is almost always within the uncertainty over the triplicate runs.
Further tests were conducted on the relative free energies between all possible 10 pairs among the 5 biologically relevant ions (Co(II), Cu(II), Mn(II), Ni(II), and Zn(II)).One example of the thermodynamic cycles using TI for the ΔG 2 leg and PMF for the ΔG 1 and ΔG 3 legs is given in Figure 4.The associated PMFs related to ΔG 1 and ΔG 3 are given in Figure S3, and the remaining thermodynamic cycles are given in Figures S4∼S9.For convenience, a master table of all the data from Figures S4∼S9 is provided in Table S1.
Overall, two end points were evaluated in this study: (1) Whether the metal-imidazole(s) interaction free energy matches experimental values and (2) whether the thermodynamic cycles close.By comparing Table 3A with Table 3B, we find that the newly derived parameters greatly outperform the earlier ones using longer simulation times.Importantly, goals (1) or (2) were Table 2. Parametrized Polarizability (Å 3 ) of Nitrogen, Final C 4 Value (kcal/mol•Å 4 ), and Energy (kcal/mol) Comparison between the Extended Sampling Parameter, One-Imidazole-Compatible 12-6-4, and Default 12-6-4 Parameter Sets a,b a Bolded data for Ni(II) are associated with the PMF in Figure 3, and the rest of that column has their PMFs given in Figure S1.Bolded data on the left are associated with PMFs given in Figure S2, and the rest of that column is associated with PMFs given in Figure S3 because these five ions were used in the thermodynamic cycle convergence analysis.b Colors used in the column headers match the color scheme in Figure 2B.Journal of Chemical Theory and Computation not as well achieved using parameters derived using 4 ns umbrella sampling, as shown in Tables 2 and 3. Overall, the extended sampling parameter sets supersede the earlier parameters and should be used in all cases.If only one imidazole is involved, the earlier parameter set can be used without significant errors, but, that said, many protein systems have multiple-imidazole side chains from His residues, making it essential to use the present parameter set. 54,83.3.QM Calculation and Metal-Imidazole Distance Comparison.As a final validation of parameter reliability, DFT optimizations were conducted to obtain the metal-imidazole distances since no experimental distance values were available.The results are listed in Table 4.By comparing the left-side "PMF" column and the right-side "DFT" columns, the results show that the PMF-calculated distances strongly agree with the DFT calculated results with different combinations of density functional and basis set.A few exceptions still exist, such as Ag(I) and Cu(II), which always have their distances overestimated by PMF studies compared to QM, and Mn(II), which, on the other hand, has the distance values underestimated by PMF studies.For Ag(I), this is likely because the coordination number for Ag(I) is commonly overestimated in MD (see Figure S1).For Cu(II), significant Jahn−Teller axial distortions are common, making reproducibility with these spherical models difficult.The Mn(II) deviations are, however, likely due to the higher relative error from weaker binding energies rather than ligand field stabilization.For example, Hexakis(imidazole)cobalt(II) cation is high spin with a low-lying low spin excited state.84 Therefore, access to the Hexakis(imidazole)cobalt(II) cation low-spin excited state, which has Jahn−Teller bond length distortions, 85 leads to variations in PMF bond lengths.The Hexakis-(imidazole)manganese(II) cation has also been shown to be high-spin.,86,87 so with access to the low-spin state, Jahn−Teller distortions would also be possible, resulting in bond length variations for the Mn complex.
For Zn(II) especially, the SAPT calculation was conducted following the flowchart displayed in Figure 5, where the SAPT calculations were performed for each of the Zn-(imidazole) 6−n (H 2 O) n systems, with n = 0∼5.5, following the methods described in Section 2.5 and Figure 5. Comparing Tables 5 and 2, it can be concluded that the average C 4 determined by SAPT is of similar magnitude to the C 4 used in the simulation work (e.g., Zn(II) C 4 is 424 kcal/mol•Å; 4 the average value of the SAPT calculated C 4 at the equilibrium distance is 470.8 kcal/mol•Å 4 ; average C 4 is 411.9 kcal/mol•Å 4 for the fit over 1.7 to 3.0 Å).Overall, we can make two conclusions: First, the C 4 value obtained using the SAPT method likely includes ion-induced dipole and also charge transfer terms.Therefore, the imidazole number will affect both the energy and metal-ligand distance significantly.Second, for the PMF parametrized C 4 , modeling the system with a constant value is computationally convenient, but will be less effective when facing spin state changes or Jahn−Teller effects in Cu(II).For Co(II) and Mn(II), the ground-state spin can change as the ligand field shifts from M(H 2 O) 6 2+ to M(imidazole) 6 2+ .Despite the changes in ligand field stabilization energy, and enhanced imidazole affinities, the 12-6-4 model of Co(II) can well reproduce the experimental binding affinities with imidazole and agrees well with the DFT-derived bond lengths.This suggests the extended 12−6−4 model qualitatively accounts for ligand field stabilization in divalent transition metal ions.Additional work is needed to increase the precision of 12−6−4 parameters and nonspherical bond lengths derived from Jahn−Teller distortions for Cu(II) and Co(II).Ligand Field Molecular Table 3. Binding Free Energy Calculated by PMF (Green and Blue, ΔG 1 and ΔG 3 in Figure 4) a,b a From top to bottom: (A) using extended sampling PMF on extended sampling parameters, (B) using extended sampling PMF on one-imidazolecompatible parameters, and (C) using regular PMF on one-imidazole-compatible parameters.b Relative energy calculated by TI (cyan, ΔG 2 and ΔG 4 in Figure 4).Thermodynamic cycle value (yellow, ΔΔG in Figure 4) for the Co(II) − Ni(II) pair.

CONCLUSIONS
In the present work, we have parametrized the 12−6−4 Lennard-Jones (LJ) nonbonded model for 11 biologically relevant metal ions for HID (δ-nitrogen protonated)-imidazole for the OPC water model.Among these 11 ions, five had their parameters rigorously tested using thermodynamic cycles, which demonstrated that we could both reproduce the experimental binding free energy and close the thermodynamic cycles with up to six imidazole molecules coordinating the central metal ion.Moreover, using DFT calculations, we show that the obtained C 4 values reproduce M(II)-imidazole distances well.Overall, with the help of sufficient AMBER24 100 tutorials, the obtained metalimidazole parameter sets will provide the basis to conduct accurate MD simulations of metal-imidazole complexes at a higher level of accuracy than previously possible, both from an energy and geometric perspective.This work is also compatible with the previously established metal-acetate 101 and metalphosphate 102 parameters to systematically describe ion behavior in complex protein systems.The main current limitation is the lack of parameters for ion-thiol or ion-thiolate interactions, which is largely due to the lack of experimental thermodynamic information for this class of interactions.It is worth noting that proteins that use thiol/thiolate metal interactions to bind transition metal ions make up a relatively small portion of the PDB. 103Importantly, the majority of ion-protein interactions are still constructed by ion-imidazole interactions, ion-acetate interactions and ion-phosphate interactions, and they present validated parameter sets that open up the possibility of accurately studying these systems.

■ ASSOCIATED CONTENT
12−6−4 Potential.In this work, we have used a 12−6−4 nonbonded model along with the AMBER force field:

Figure 1 .
Figure 1.(A) Illustration and comparison of the charge distribution for both HID imidazole molecules.Heavy atom-type names in AMBER are also marked, where NA denotes the atom type of δ-nitrogen, and NB denotes the atom type of ε-nitrogen.Here, Zn(II) was used as an example, so the distance between ion and ε-nitrogen is 2.06 Å. (B) Illustration of the parametrization process.

Figure 2 .
Figure 2. (A) Visualization of the data fromTable1.(B) Visualization of the energy data from Table2, where extended sampling reveals the limit of one-imidazole parameters (red), which was previously not discovered by regular sampling (blue).This issue was resolved by introducing multipleimidazole parameters (green).Color matches the color scheme for the column headers in Table2.

Figure 3 .
Figure 3.Comparison among AMBER, PLUMED, extended sampling AMBER, and extended sampling PLUMED PMF results using Ni(II) and one HID imidazole in OPC water as an example.Green arrow indicates the "cation-pi-stacking" conformation captured by the extended sampling methods for both old and new parameters.PMFs for other ions are available in Figure S1.Only one PMF of the three duplicates is presented.

Figure 4 .
Figure 4. Thermodynamic cycle of mutating Co(II) to Ni(II), while the ion has one imidazole coordinating; here, ΔG 4 was direct taken from a previous work.41

3 . 4 . C 4
Dependence On Ligand Number and its Correlation with Low-High Spin States.According to the induction energies calculated based on SAPT, we have found that the ion-induced dipole interaction term, which represents C 4 in the 12−6−4 model, shows a strong dependency as the imidazole number increases.The result is displayed in Table

a
If multiple-imidazole molecules are present, the average distance was taken.High spin states were used since they are energetically favored under the current level of theory and basis sets.Gaussian 16 was used. 89b Metal-hexaimidazole distance data are presented in "equatorial, axial" format if applicable.c DFT optimization on MD-generated structures for faster optimization convergence.d Optimized distances at the B3LYP-D3/6-31G* level of theory.These equilibrium distances were used to calculate the C 4 values based on SAPT.

Figure 5 .
Figure 5. Illustration of the derivation of C 4 terms based on SAPT calculations using Psi4.88

Table 1 .
Experimental Interaction Energy of Ions With One Imidazole in Aqueous Environment a a Related data are visualized in Figure2.

Table 4 .
Distances between the Metal Ion and ε-Nitrogen of the Imidazole a