BILFF: All-Atom Force Field for Modeling Triazolium- and Benzoate-Based Ionic Liquids

We present an extension of our previously developed all-atom force field BILFF (Bio-polymers in Ionic Liquids Force Field) to three different ionic liquids: 1-ethyl-3-methyl-1,2,3-triazolium acetate ([EMTr][OAc]), 1-ethyl-3-methyl-1,2,3-triazolium benzoate ([EMTr][OBz]), and 1-ethyl-3-methylimidazolium benzoate ([EMIm][OBz]). These ionic liquids are of practical importance as they have the ability to dissolve significant amounts of cellulose even at room temperature. Our force field is optimized to accurately reproduce the strong hydrogen bonding in the system with nearly quantum chemical accuracy. A very good agreement between the microstructure of the quantum chemical simulations over a wide temperature range and experimental density data with the results of BILFF were observed. Non-trivial effects, such as the solvation shell structure and π–π stacking of the cations, are also accurately reproduced. Our force field enables accurate simulations of larger systems, such as solvated cellulose in different (aqueous) ionic liquids, and is the first to present the optimized parameters for mixtures of these solvents and water.


Introduction
Ionic liquids (ILs) have become an attractive alternative to conventional solvents in various industries due to their unique properties like their non-volatility, non-flammability, and their very good thermal and chemical stability [1,2].The properties of ILs can be tailored to suit specific applications by changing the cation or anion components [3,4].This flexibility has led to their use in a range of different applications, including in catalysis, electrochemistry, material science, energy, and biotechnology [5][6][7][8][9][10][11].One of the promising usage of ILs is the dissolution and processing of cellulose [12][13][14][15][16][17], which enables the production of cellulose-based derivatives in a homogeneous manner [18,19] and the processing of biomass by the separation of wood components [20,21].
Cellulose is characterized by its highly structured composition of glucose monomers, forming long chains that can range from several hundred to tens of thousands of units in length.These chains are highly ordered and arranged in parallel.Intramolecular hydrogen bonding stiffens the polymer chains, while intermolecular hydrogen bonding facilitates the arrangement of these linear polymers into sheet-like structures.These sheets are densely packed together through hydrophobic interactions and form crystalline structures [22,23].To dissolve cellulose in ILs, the intermolecular hydrogen bonds between the cellulose strands must be disrupted and replaced by interactions with the solvent.This can be achieved through the formation of new hydrogen bonds to the solvent molecules.Several studies have demonstrated that the effectiveness of ILs in dissolving cellulose depends on several factors such as the type of cation and anion, the basicity as well as the position and length of the alkyl chains of the cation [15][16][17]24,25].
The two cations, [EMTr] + and [EMIm] + , are distinguished by the replacement of one ring carbon atom in [EMTr] + by a third nitrogen atom.As a result, the third ring proton in the 1,2,3-triazolium is missing.This structural change leads to a lower probability of proton abstraction and thus to a lower formation of N-heterocyclic carbenes (NHC).NHCs are highly reactive and form dimers or undesirable byproducts.Thus, they can also cause the unwanted decomposition of cellulose [26].The potential use of triazoliumcontaining ILs as non-derivatizing cellulose solvents is nearly unexplored and was first described by Brehm et al. in 2019 [15].By understanding the complex mechanisms involved in cellulose dissolution in different ILs, it is possible to develop more efficient and sustainable processes for the production of cellulose-based materials.This requires a thorough understanding of the interactions between ILs and cellulose, which can be achieved through molecular dynamics (MD) simulations.However, the complexity of the system, together with the problem that solubility processes take place on time scales that exhaust the resources of quantum chemical MD simulations, necessitates the use of force field MD (FFMD) simulations.To accurately simulate the interactions between ILs and cellulose, it is crucial to have well-optimized force field parameters for the individual ions.Until now, there has been no optimized force field available for the ions [EMTr] + and [OBz] − , which limited the accuracy of simulations involving these ILs.In this article, we present the development of a force field for (aqueous) [EMTr][OAc], [EMIm][OBz] and their combination [EMTr][OBz], which is an extension of our previously published force field for [EMIm] + and [OAc] − called BILFF (Bio-Polymers in Ionic Liquids Force Field) [27,28].The force field development is focused on the accurate reproduction of the microstructure of the ILs, especially the hydrogen bonds, based on results of ab initio MD (AIMD) simulations.To optimize the force field parameters, we compared the results of different analysis like radial and combined distribution functions of the hydrogen bonds along with statistical analyses of bond lengths, angles, and dihedral angles with the results from AIMD simulations for both the pure and aqueous ILs.Based on this comparison, the force field parameters were adjusted and optimized for a simulation of the ILs both in the absence and presence of water at 350 K.For the aqueous system, a ratio of IL ion pairs to water of 1:3 was chosen, at which cellulose is already precipitated again [29][30][31].
The results demonstrate that the new force field effectively reproduces the structural and dynamic properties of the ILs, including ion pairing, solvation shell structure, and hydrogen bonding with respect to their geometry and lifetime, over a wide temperature range.Furthermore, we contrasted and compared the results of the three ILs The optimized force field parameters enable the simulation of larger cellulose systems in various ILs on larger time scales, which would otherwise be too computationally expensive when using AIMD simulations.This enables a deeper understanding of the underlying interactions for the solubility of cellulose in the different ILs and thus the identification of novel and improved cellulose solvents.

Optimization Procedure
In a previously published article [27], the optimized force field parameters for (aqueous) [EMIm][OAc] were presented.The aim of this article is the extension of the parameter set to the ions [EMTr] + and [OBz] − using the already optimized parameter set for were performed.Since it is known that many imidazolium-and triazolium-based ionic liquids are hygroscopic and that this might drastically affect their properties [32,33], the optimization of the force field parameters was carried out simultaneously in the presence and absence of water.For this purpose, simulations of pure [EMIm][OBz] and [EMTr][OAc] as well as the combination of both ions as [EMTr][OBz] in water, were performed (cf. Figure 1).The latter simulation is used to investigate both parameter sets of [EMIm][OBz] and [EMTr][OAc] in the presence of water.The optimization of the force field parameters was performed by comparing various analyses of force field MD simulations with reference AIMD simulations of these three systems, which are discussed in more detail in the following sections.The aim of optimizing the force field parameters for [EMTr] + and [OBz] − in the presence and absence of water is to enable force field MD simulations of these ions with nearly the accuracy of quantum chemical calculations to compute larger systems with less computational effort.
The focus was set on the correct reproduction of the hydrogen bonds between the investigated ions [EMTr] + as well as [OBz] − and their counterions [OAc] − and [EMIm] + calculated from quantum chemical simulations.[EMIm] + and [OAc] − were chosen as counterions for [OBz] − and [EMTr] + , since optimized force field parameters are already available for these two ions (see our previously published articles on BILFF [27,28]).
The hydrogen bonds in the aqueous system were investigated using an AIMD simulation of aqueous [EMTr][OBz].A molar fraction of χ IL = 0.25 was simulated in order to cover a wide range of applications as it is known that, at this water content, cellulose is already precipitated again in other ILs such as [EMIm][OAc] [29][30][31].
The force field is based on the functional form of the OPLS-AA force field [34][35][36]: i,j = i j , σ i,j = σ i σ j . ( In order to optimize the force field parameters, radial distribution functions of the reference AIMD simulations were compared with the results of the force field MD simulations.To retain the properties of the underlying force fields as much as possible, the force field parameters q and σ were iteratively adjusted using a trial-and-error approach.The force field parameters were varied until the deviations of the RDFs were minimized.The results are discussed in Section 2.1.1.In addition, σ was optimized with respect to a good reproduction of the experimental density (see Table 7).
The bonded interactions were adjusted based on a comparison of the statistical occurrence of bond lengths, bond angles and dihedral angles, and by varying the force constants and height of the potential barrier.The nomenclature of the atom types can be found in Figure 2. To correctly reproduce the reference AIMD simulations, the equilibrium bond length of the NA-NR bond in [EMTr] + , for example, was modified by +0.3%.The NA-CW bond length was reduced by about 2%, while the bond force constant was increased by about 12% compared to the value in [EMIm] + .
In [OBz] − , the greatest adjustment in the bond length and force constant occurred in the CO-O2 bond, with changes of about +2% and +56%, respectively.The equilibrium bond angle, on the other hand, was not changed for both ions.
A flowchart of the force field optimization process can be found in Figure 3.The force field parameters and nomenclature for the atom types and classes can be found in Tables 1-5.

Adjustment of the Nonbonded Force Field Parameter
Comparison of the results with reference AIMD and experimental density and variation of partial charge and Lennard-Jones parameter

Adjustment of the Bonded Force Field Parameter
Comparison of the statistical occurrence with reference AIMD and variation of the force constants and potential barrier height Good Agreement? No. Yes.

Good Agreement?
No.

Yes.
(Adapted) Literature Force Field BILFF Figure 3. Flowchart of the optimization process of the force field parameters for BILFF.As a starting point for the optimization of the force field parameters, BILFF [27,28] for [EMTr] + was used.For [OBz] − the Lennard-Jones parameters and force constants were adapted from OPLS-AA [34][35][36].The partial charges were calculated using the restrained electrostatic potential (RESP) methodology.The equilibrium bond length was calculated using a geometry optimization.2), while the atom classes are used for the bonded interactions (see Tables 3-5).

Atom Type
Atom Class   In the following, the radial distribution functions (RDFs) in the different ILs are compared with the results of the corresponding AIMD simulations.Furthermore, conclusions are drawn about the differences in the formation of hydrogen bonds between the cations [EMTr] + and [EMIm] + and their respective bonding partners, as well as between the anions [OBz] − and [OAc] − and their respective hydrogen bond donors.First, the hydrogen bonds of the cation [EMTr] + are analyzed (see Figures 4-6).The ring protons of [EMTr] + form a strong hydrogen bond to the anion [OAc] − with a particle distance of 195 pm (cf. Figure 4) and a g(r) = 5.4, as can be seen from the first maximum of the RDF.The agreement of the RDF with the AIMD result is excellent.To achieve this, the partial charge of the ring protons was increased by about 27%, from 0.150 to 0.191, compared to the force field parameters for [EMIm] + [27].The value of σ was decreased from 1.62 Å to 1.48 Å by about 12%, reducing the atom repulsion.Furthermore, it was investigated whether the two aromatic ring protons in both cations differ in terms of the intensity of hydrogen bonding.The corresponding RDFs are shown in Figure 6 (black/grey curve and red/orange curve).An almost identical behavior of the curves can be observed, which occurs due to the equivalence of the ring protons in each of the two molecules.

Torsion Angle
In the following, the hydrogen bonds of the second considered ion [OBz] − are analyzed (cf.Figures 7-9).With a g(r) of 3.7 and a particle distance of 205 pm in the first maximum of the RDF, the anion forms strong hydrogen bonds to [EMIm] + (cf. Figure 7).By increasing the partial charge of the oxygen atom of [OBz] − by about 5% from −0.524 to −0.550, and reducing the atom repulsion by modifying σ by about 5% from 2.96 Å to 2.80 Å compared to the adapted force field parameter of OPLS-AA [34][35][36] (cf.Table S1), an excellent agreement with the AIMD simulation was obtained.In contrast, the adapted literature force field [34][35][36] shows slight deviations in the position and height of the first maximum.However, benzoate forms an even stronger hydrogen bond to the protons of water than to [EMIm] + (cf.Figures 7 and 8).Besides the obvious difference in the O-H vs. C-H dipole moment of water and [EMIm] + , another reason can be steric effects.The oxygen atoms of benzoate are more easily accessible by the small water molecules than by the large imidazolium cation.
In comparison, the adapted literature force field gives rise to a significantly weaker peak intensity of the corresponding RDF (cf.Figures 7 and 8, green curve).
In addition, a comparison of the cation-anion hydrogen bonding intensities of all four ion pair combinations is analyzed in Figure 9.The first maximum in the RDF of in the anhydrous system (black and grey curves).This indicates that [OBz] − forms slightly stronger hydrogen bonds to the cation than [OAc] − , as already observed before [40].We can therefore confirm the previous reports in the literature.This effect is even more pronounced in the triazoliumcontaining ILs with a difference in the first maximum of a g(r) of 1.1.The same trend is observed when comparing the corresponding RDFs of the AIMD simulations (cf. Figure S3 in the Supplementary Materials).
In the presence of water, the respective cation• • • anion hydrogen bonds of both anion species are weakened (cf. Figure 9, red and orange curves).This is caused by the additional formation of an anion• • • water hydrogen bond, which is much stronger than the cationanion hydrogen bond (cf. Figure S4 in the Supplementary Materials).
Due to the stronger hydrogen bond, the hydration shell around [OBz] − is more pronounced than that of [OAc] − .As a result, there is a greater attenuation of the cation• • • [OBz] − interaction relative to the cation• • • [OAc] − hydrogen bond in the presence of water.
The hydration shell around the ions, as well as the coordination numbers, are discussed in Section 2.2.
In summary, the following trends are thus obtained from the comparison of the RDFs of the different ILs:

•
Triazolium forms stronger hydrogen bonds with both the anion and water than imidazolium.

•
Benzoate forms stronger hydrogen bonds with the cation than acetate in the anhydrous system.

•
Benzoate forms stronger hydrogen bonds with water than acetate, resulting in a greater attenuation of the cation-anion interaction in the aqueous system.

Combined Distance-Angle Distribution Functions
While radial distribution functions provide an important initial assessment of the quality of the simulated hydrogen bonding structure, the directional aspect should not be neglected.
Figure 10 shows the combined distance-angle distribution functions (CDFs) of the distance H EMIm • • • O OBz and the corresponding angle (CH EMIm , H EMIm , O OBz ) as illustrated in the sketch.The overall pattern of the CDFs of the AIMD simulation (top) and the FFMD simulation with BILFF (bottom) agree very well.The main hydrogen bonding peak at a particle distance of about 200 pm shows a broader angular distribution in the BILFF simulation, which may be the result of better statistical sampling compared to the AIMD simulation.
The second peak at about 400 pm originates from both the "second" oxygen of [OBz] − and the "second" hydrogen of the [EMIm] + molecule.Again, the AIMD simulation yields a slightly stronger localization than the force field.
The third peak, observed at approximately 600 pm, occurs at an angle of around 15 • and corresponds to the hydrogen bonding involving the isolated CH group, leading to an increased presence of benzoate ions on the "upper" side of the molecule.
Considering the CDFs of the 1,2,3-triazolium cation and the anions acetate and benzoate, a similar good agreement of the CDFs of the AIMD and FFMD simulations can be obtained (cf.Figures S8 and S9 in the Supplementary Materials).Additional CDFs regarding the cation• • • anion, cation• • • cation, and anion• • • water interactions in [EMIm][OAc] can be found in the previously published article [27].
However, for a correct reproduction of the microstructure of the AIMD simulation, not only the hydrogen bonds are important, but also the description of the π-π interactions between the aromatic rings of Considering the homomolecular distance-angle distribution functions of the aromatic cations, [EMTr] + and [EMIm] + (cf. Figure 11, left panel and Ref. [27]), symmetric peaks are observed at 400 pm, 0 • and 400 pm, 180 • , but no residence probability is found around 400 pm, 90 • .This pattern indicates the occurrence of π-π stacking.The absence of a significant probability for the cations being at a 90 • angle to each other implies that a T-shaped arrangement between the cations is unlikely.

Spatial Distribution Functions
Providing valuable insights into the overall structural arrangement, spatial distribution functions (SDF) provide a comprehensive picture of the relative orientations of molecules within the complex system involving the different cations, anions, and water.Figure 12 shows the results of the force field MD simulation using BILFF for all four systems.The oxygen atoms of water and the anion arrange themselves competitively at the respective ring protons of the cation.Comparing the imidazolium and triazolium patterns, the isolated aromatic ring proton of [EMIm] + acts as a strong attractor for the oxygen atoms, resulting in deflation around the aliphatic side chains, while the substituted nitrogen is avoided, resulting in a much more pronounced oxygen density around the methyl/ethyl side chains.) and the oxygen atom of the anion (green, 7 nm −3 ) as well as the ring center of the cation (red, 7 nm −3 ) and anion (cyan, 10 nm −3 ) in all four systems calculated with BILFF.
The spatial distribution function also shows distinct π-π stacking motifs above/below the aromatic cations.Interestingly, the stacking partner is determined by the anion type; for [OAc] − a homomolecular cation-cation stacking (red areas in Figure 12) is observed, while for [OBz] − the corresponding regions are preferentially occupied by the benzoate anions (cyan areas).The ring stacking of the benzoate is in agreement with results in the literature [16].The results from the force field MD simulation of aqueous [EMTr][OBz] also agree well with the results obtained in the reference AIMD simulation (cf. Figure S10 in the Supplementary Materials).

Competing Hydrogen Bonds
In the ternary mixtures of ionic liquids in water, all direct interactions between the ions are brought into competition with hydrogen bonding interactions with water.The topology of the interaction pattern is analyzed in the form of a Sankey diagram (cf. Figure 13).The coordination numbers were calculated from the integral under the curve of the corresponding RDFs to the first minimum.These numbers indicate the number of hydrogen bonds per atom, which is proportional to the thickness of the bars in the Sankey diagram.Considering the anions, benzoate is surrounded by about 4.5 other molecules of the cation or water in the corresponding ILs (cf.Figures 12 and S12 in the Supplementary Materials).The smaller acetate, on the other hand, has a slightly higher coordination number of about 5 (cf.Figures 12 and S12 in the Supplementary Materials).
Compared to the results for [EMTr][OBz] from the AIMD simulation, the coordination numbers of BILFF can be well reproduced, especially for the cation.It should be noted, that the force field parameters were not explicitly optimized for this combination, but for [EMTr][OAc] and [EMIm][OBz].

Radial Distribution Functions at Higher Temperatures
To verify the transferability of BILFF to [EMTr][OBz] force field MD simulations were performed and compared with reference AIMD simulations.Also, the application of the force field at elevated temperatures has been investigated and is discussed in this section.The radial distribution functions (RDFs) of the hydrogen bonds between the ring protons of the cation and the oxygen atoms of the anion are shown in Figure 14.
At 350 K, the first maximum of the RDF occurs at a particle distance of 205 pm and a g(r) = 2.5.As expected, with increasing temperature, the maximum decreases slightly to a g(r) = 2.1 at a temperature of 550 K.As the temperature rises, the mobility of the molecules increases, so that the hydrogen bonds become more short-lived, resulting in a smaller g(r).
Similar to the AIMD simulation, the height of the maxima in the RDFs of the force field MD simulation does not change drastically with increasing temperature.Although the force field parameters were not directly optimized for aqueous [EMTr][OBz], the agreement with the reference AIMD simulation is very good, even at high temperatures, underlining the high accuracy of the optimized force field parameters of BILFF.

Hydrogen Bond Lifetime
Up to this point, a number of structural quantities have been computed, mainly expressed as distribution functions.Exemplary dynamical properties have also been calculated, in particular the hydrogen bond lifetime.
In Table 6, the hydrogen bond lifetime calculated from the AIMD and FFMD simulation using BILFF are compared.Furthermore, a comparative analysis of the lifetimes of the different hydrogen bonds in the investigated ILs The calculation of the lifetime is differentiated in the intermittent and continuous calculation whereby, in the intermittent calculation, the breaking and reforming of the hydrogen bond is allowed.The continuous lifetime, on the other hand, represents the duration until the very first breakage of the hydrogen bond.To calculate the lifetime, geometric criteria were defined on the basis of the first maximum in the underlying distance-angle distribution functions.These criteria are listed in Table S5 in the Supplementary Materials.In Figure 10, the definition of the geometric criterion was exemplified as a black rectangle.The lifetime values presented in Table 6 are, respectively, averaged over all ring protons of the imidazolium (three ring protons) and 1,2,3-triazolium cation (two ring protons).The lifetime of the hydrogen bonds of the single protons H Me and H Eth can be found in Table S6 in the Supplementary Materials.(The lifetime of the hydrogen bonds [EMIm] + • • • [OAc] − as well as [EMIm] + • • • water, with respect to all three individual ring protons can be found in the previous article on BILFF [27]).a Calculated from the MD simulations of our already published article [27].
The following conclusions can be drawn: The longest hydrogen bond lifetime is observed in the anhydrous system, primarily due to strong Coulomb interactions, with [EMTr] + forming more durable hydrogen bonds with the anion than [EMIm] + .However, this trend is reversed in the aqueous ILs.
The introduction of water results in the reduction of the hydrogen bond lifetime between the cations and anions in all ILs, due to a competition for oxygen atoms surrounding the cation, as well as a lower viscosity.This competition between water and the anion can also be observed in the radial and spatial distribution functions (see Sections 2.1.1 and 2.1.3)as well as in the Sankey diagrams (see Section 2.2).At the same time, the anion [OBz] − forms more long-lasting hydrogen bonds to water than [OAc] − .This is also reflected in the frequency of the occurrence of the hydrogen bond of the ILs measured by the height of the first maximum of the RDF.
However, the lifetime of the hydrogen bonds can be affected not only by water, but also by an increased temperature.Due to the resulting higher mobility of the molecules, the lifetime of the hydrogen bonds decreases, as seen in aqueous [EMTr] [OBz].
In all four ILs, the first breakup of the hydrogen bonds take place after only a few picoseconds (about 1-7 ps).Therefore, the two cation ring protons H Me and H Eth show similar values of lifetime (see Table S6 in the Supplementary Materials).In BILFF, they are regarded as equivalent.
It should be noted that dynamical properties, such as lifetimes, depend exponentially on the potential energy surface and are therefore subject to larger variations (with respect to the experiment, but also with respect to variations of the theoretical methods).From this perspective, the agreement between the force field MD simulation with BILFF and the AIMD simulation is still very good.It is also noteworthy that some of the lifetimes are in excess of the total simulation time of the AIMD simulation.

Validation of the Density and Diffusion
To test the reproducibility of experimental quantities, one approach is to compare the calculated system densities of the different ILs with corresponding experimental data (see Table 7).As expected, a decrease in density is observed in the presence of water.The mobility of the molecules also plays a crucial role in accurately describe the behavior of the ILs.Below, the self-diffusion coefficients in all investigated ionic liquids calculated from the force field MD simulations are compared with the reference AIMD simulation (see Table 8).The following observations can be made: a Calculated from the MD simulations of our already published article [27].
• According to the FFMD results, the bulky cation generally diffuse faster than the anion [42].An exception is [EMIm][OAc] in water, where the anion diffuses faster.This phenomenon is already known and is in agreement with experimental data for [EMIm][OAc] [43].

•
The diffusion coefficient is higher at elevated temperatures.

•
In the presence of water, the diffusion coefficient of all ionic liquids increases significantly.

•
The diffusion coefficient is influenced by both the cations and anions in the IL.When comparing the diffusion of the cations in the different ILs containing either [EMTr] + or [EMIm] + , it is observed that the triazolium cations exhibit a slower diffusion compared to imidazolium.Furthermore, when comparing the diffusion rates of the anions, it is found that benzoate diffuses slower than acetate in both pure and aqueous systems.This is due to the different occupied volumes and thus the bulkiness of the molecules as well as their more pronounced hydrogen bond formation.

•
The lowest diffusion rate, and, at the same time, the highest density, can be found in pure The results of the force field MD simulations reflect the trends of the reference AIMD simulations for the imidazolium-containing ILs.In the case of the triazolium-containing ILs, the AIMD simulation shows a slightly faster diffusion of the anion compared to the cation at 350 K and 450 K, and thus an opposite trend to the FFMD.Despite occasional numerical fluctuations in the specific diffusion values, the overall agreement of this dynamic quantity with the first-principles reference simulation is very good.The temperature dependence is reproduced over more than one order of magnitude.It should be noted that [EMTr][OBz] has been synthesized recently, so no experimental data for the diffusion coefficient are available yet.

Computational Details
In this section, the computational details of the underlying reference AIMD simulations and force field MD simulations are presented.As a starting configuration for the AIMD simulations, the last time step of a force field MD simulation using the bonded and Lennard-Jones parameters of OPLS-AA [34][35][36] and the CL&P force field [37][38][39] and RESP partial charges (see Ref. [15]) as well as TIP4P-EW [45] (with constrained bonds and angles using the SHAKE algorithm [46,47]) for water was used.The physical simulation time was 5 ns for [EMIm][OBz], 0.5 ns for [EMTr][OAc] and 20 ns for aqueous [EMTr][OBz].The final AIMD simulations were performed in analogy to the previous articles on BILFF [27,28] with CP2k [48][49][50] using the quickstep method [51] and orbital transformation (OT) [52].Under application of the density functional theory [53,54] using the BLYP functional [55,56], dispersion correction D3(BJ) of Grimme [57,58] and the revised damping parameters of Smith et al. [59], the electron density was calculated.The DZVP-MOLOPT-SR-GTH [60] and GTH pseudopotentials [61,62] with a plane-wave energy cutoff of 350 Ry were used as basis sets.To equilibrate the AIMD simulations, the first 30 ps for [EMTr][OAc] and the first 50 ps for the systems with [OBz] − were discarded, respectively.The physical simulation time, size of the simulation box, and particle numbers can be found in Table 9.
The force field MD simulations for optimizing the force field parameter were also performed analogously to the previous articles on BILFF [27,28].After constructing the simulation box using Packmol [63], an equilibration was performed at 500 K in the NVT ensemble using a Berendsen thermostat [64] (coupling constant of 1.0 fs) in a 25 ps simulation followed by a 100 ps simulation at a temperature of 350 K in the NpT ensemble using a Nosé-Hoover thermostat [65][66][67] (coupling constant of 100 fs) and a Nosé-Hoover barostat (coupling constant of 2000 fs).Using the Langevin thermostat [68,69], the acoustic shock waves were damped after the equilibration runs and adjustment of the final box volume.The production run was finally computed in the NVT ensemble using the Nosé-Hoover thermostat [65][66][67] (coupling constant of 100 fs) and a time integrator of 0.5 fs.The PPPM long-range solver of LAMMPS [70] was applied with a Lennard-Jones cutoff radius of 800 pm to calculate the electrostatic interactions.
As a starting point for the optimization of the force field parameters BILFF [27,28] for [EMTr] + was used.The reference bond length and angle of atom type NR were determined using a geometry optimization.For [OBz] − , the Lennard-Jones parameters and the force constants were adapted from OPLS-AA [34][35][36].The partial charges were calculated from quantum chemical calculations via the restrained electrostatic potential (RESP) methodology.The equilibrium bond length was calculated using a geometry optimization.The starting force field parameters are listed in the Supplementary Materials.The ionic charge was set to ±0.82, as in the force field development of [EMIm][OAc] [27].Reducing the total ion charge is discussed as an alternative to the usage of computationally more expensive polarizable force fields in the literature and is widely applied [71][72][73][74][75].
To validate the force field, additional force field MD simulations of pure and aqueous [EMTr][OBz] at 350 K as well as force field and reference AIMD simulations of aqueous [EMTr][OBz] at elevated temperatures of 450 K and 550 K were calculated.It should be noted that, for this IL, the force field parameters have not been optimized and the transferability of the parameters for [EMTr] + in [EMTr][OAc] and [OBz] − in [EMIm] [OBz] to this IL will be tested.The same simulation parameters as described above were used.The first 30 ps of the AIMD simulations were again discarded for equilibration.The resulting physical simulation time is 88 ps in each case.An analysis of the temperature dependence of the simulations can be found in Section 2.3.

Conclusions
In this article, we present the extension of our all-atom force field BILFF (Bio-Polymers in Ionic Liquids Force Field) to the ionic liquids [EMTr][OAc], [EMIm][OBz] and [EMTr][OBz] in absence and presence of water.BILFF was previously tuned for an optimal balance of the competing interactions in the ternary mixture of [EMIm][OAc] and water in view of its application to the challenge of modeling the solvation structure of cellulose (as solute) by [EMIm][OAc]/water [27,28] [EMTr][OAc], [EMIm][OBz] and [EMTr][OBz] with experimental density data, as well as with the previously published results for [EMIm][OAc].

Figure 4 .
Figure 4. Comparison of the RDFs of the hydrogen bond between the marked oxygen atoms of [OAc] − and the ring protons of [EMTr]+ calculated from a reference AIMD and force field MD simulations using adapted literature force field parameters[37][38][39] and BILFF.The RDFs are averaged over both ring protons.

Figure 5 .
Figure 5.Comparison of the RDFs of the hydrogen bond between the marked oxygen atom of water and the ring protons of [EMTr] + calculated from a reference AIMD and force field MD simulations of aqueous [EMTr][OBz] using adapted literature force field parameters [34-39] and BILFF.The RDFs are averaged over both ring protons.

Figure 6 .
Figure 6.Comparison of the RDFs of the hydrogen bond between the marked oxygen atoms of [OAc] − /water and the ring protons of [EMIm] + /[EMTr] + calculated from force field MD simulations of aqueous [EMTr][OAc] and [EMIm][OAc] using BILFF.According to the height of the RDFs (Figures 5 and 6), the hydrogen bond between [EMTr] + and the oxygen of water is comparatively weaker than the hydrogen bond between [EMTr] + and [OAc] − .A comparison of the hydrogen bond of [EMTr] + and H 2 O in aqueous [EMTr][OBz] between an AIMD and force field MD simulation with BILFF is shown in Figure 5.It should be noted that the force field parameters were specifically adjusted for [EMTr] + in [EMTr][OAc] and [OBz] − in [EMIm][OBz], but not for the combination of both ions as [EMTr][OBz].However, the results of the RDF calculated from AIMD simulations

Figure 7 .
Figure 7.Comparison of the RDFs of the hydrogen bond between the marked oxygen atoms of [OBz] − and the ring protons of [EMIm]+ calculated from a reference AIMD and force field MD simulations using adapted literature force field parameters[34][35][36] and BILFF.The RDFs are averaged over all three ring protons.

Figure 8 .
Figure 8.Comparison of the RDFs of the hydrogen bond between the marked oxygen atoms of [OBz] − and the protons of water calculated from a reference AIMD and force field MD simulations using adapted literature force field parameters [34-36] and BILFF.

Figure 9 .
Figure 9.Comparison of the RDFs of the hydrogen bond between the marked oxygen atoms of [OBz] − as well as [OAc] − and the ring protons of [EMIm] + and [EMTr] + calculated from force field MD simulations using BILFF.The RDFs are averaged over the marked ring protons.

Figure 10 .
Figure 10.Distance-angle distribution function between an example ring proton of [EMIm] + and the oxygen atoms of [OBz] − in pure [EMIm][OBz] as a result of a reference AIMD simulation (top) and a force field MD simulation with BILFF (bottom).Color code in arbitrary units.The black rectangle demonstrates the geometric criterion for the existence of a hydrogen bond, as used in the lifetime calculation.

Figure 11 .
Figure 11.Distance-angle distribution function between two [EMTr] + (left) and two [OBz] − (right) ring centers in aqueous [EMTr][OBz] as a result of a reference AIMD simulation (top) and a force field MD simulation with BILFF (bottom).code in arbitrary units.In contrast to this, the corresponding CDF of [OBz] − • • • [OBz] − (cf. Figure 11; right panel) shows two maxima at 500 pm for angles of 0 • and 180 • , along with a broad distribution in the center region (500 pm, 90 • ).The unusually large distance between the molecules combined with the higher probability of α = 90 • shows that the coplanar [OBz] − • • • [OBz] − motif is no longer the dominant structure (as opposed to either cations).Instead, they are twisted relative to each other, adopting a T-shaped arrangement rather than a stacked configuration.In both the anion• • • anion and cation• • • cation CDFs, the results obtained from FFMD simulations exhibit a good agreement with AIMD simulations.It is important to note here that the force field parameters were developed for [EMTr] + in [EMTr][OAc] and [OBz] − in [EMIm][OBz] but not for the combination [EMTr][OBz].The fact that the microstructure of the AIMD simulation can nevertheless be reproduced well even in the case of the cross combination [EMTr][OBz] in water shows a good transferability of the force field parameters.This aspect will be discussed in more detail in Section 2.3.3.

Figure 12 .
Figure12.Spatial distribution function of the arrangement of molecules around the two different cations with the protons and oxygen atoms of water (blue, 26 nm −3 ) and the oxygen atom of the anion (green, 7 nm −3 ) as well as the ring center of the cation (red, 7 nm −3 ) and anion (cyan, 10 nm −3 ) in all four systems calculated with BILFF.

Figure 13 .
Figure 13.Visualization of the competing hydrogen bonds in aqueous [EMTr][OBz] as Sankey diagram calculated from a force field MD simulation with BILFF.The inserted numbers represent coordination numbers and are shown comparatively to the AIMD simulation.The triazolium cation exhibits a coordination number of 1.0 in aqueous [EMTr][OBz] (cf. Figure 13), and 1.21 in aqueous [EMTr][OAc] (cf. Figure S11 in the Supplementary Materials).The imidazolium cation has similar coordination numbers of 0.95 in aqueous [EMIm][OBz] and 1.15 in aqueous [EMTr][OAc] (cf.Figures S12 and S13 in the Supplementary Materials).Considering the anions, benzoate is surrounded by about 4.5 other molecules of the cation or water in the corresponding ILs (cf.Figures 12 and S12 in the Supplementary Materials).The smaller acetate, on the other hand, has a slightly higher coordination number of about 5 (cf.Figures 12 and S12 in the Supplementary Materials).Compared to the results for [EMTr][OBz] from the AIMD simulation, the coordination numbers of BILFF can be well reproduced, especially for the cation.It should be noted, that the force field parameters were not explicitly optimized for this combination, but for [EMTr][OAc] and [EMIm][OBz].

Figure 14 .
Figure 14.Radial distribution functions (RDFs) of the hydrogen bonds between the marked oxygen atoms of [OBz] − and the ring protons of [EMTr] + in aqueous [EMTr][OBz] at different temperatures, simulated by AIMD simulation and force field MD simulation using BILFF.
[EMTr][OAc], [EMIm][OBz] and [EMTr][OBz] is shown.In addition, a comparison with the IL [EMIm][OAc] is made.The lifetime of the hydrogen bonds is studied in the presence and absence of water.
Pure [EMTr][OBz] has the highest density and (aqueous) [EMIm][OAc] the lowest.This difference can be attributed to the larger molecular volume occupied by [EMTr] + and [OBz] − .The agreement of the density of pure [EMTr][OAc] as well as of the combined IL [EMTr][OBz] for which the force field parameters were not specifically optimized, is outstanding.The deviation is only 0.9% from the experimental measurement.The experimental density of [EMIm][OBz] is reproduced without any deviation at all.
. The focus of this work is the generalization of the [EMIm][OAc] force field to additional ions ([EMTr] + , [OBz] − ) including all four combinations.This article thus presents the first force field specifically optimized for [EMTr][OAc] and [EMIm][OBz].The objective was to ensure an accurate reproduction of the quantum chemical microstructure as well as dynamic properties such as hydrogen bond lifetimes and diffusion in these ILs.To achieve this, an iterative adjustment of the force field parameters to the results of reference AIMD simulations and experimental system densities was performed.In addition, to verify the transferability of the optimized force field parameters for [EMIm][OBz] and [EMTr][OAc], the results of (aqueous) [EMTr][OBz] were evaluated.By comparing the simulation results of the four ILs [EMTr][OAc], [EMIm][OAc], [EMTr][OAc] and [EMTr][OBz], the hydrogen bonding strength of benzoate was found to

Table 1 .
Nomenclature of the atom types and atom classes of [EMTr] + and [OBz] − in BILFF.The atom types are applied for the non-bonded interactions (see Table

Table 6 .
Overview of the lifetime τ of the hydrogen bonds in all four systems comparing the results of the reference AIMD simulation and the force field MD simulation using BILFF at the given temperatures (C + = Cation, A − = Anion).(No AIMD simulations of anhydrous [EMTr][OBz] as well as [EMTr][OAc]/H 2 O and [EMIm][OBz]/H 2 O have been performed, so no data are available for these).

Table 7 .
[16]arison of the system densities from force field MD simulations using BILFF at 350 K.The molar fraction of the ionic liquids in the aqueous systems is χ IL = 0.25.Measurements at 323.15 K from Ref.[15].bCalculatedfrom the MD simulations of our already published article[27].cExtrapolatedvalues of temperature dependent measurements from Ref.[41] (x IL = 0.252 in the aqueous system).d/eMeasurements at 358.15 K from Ref.[16]. a

Table 8 .
[27]-diffusion coefficients D from force field MD simulations using BILFF.For comparison, the diffusion coefficients of[EMIm][OAc] are also shown from our previously published article[27].

Table 9 .
Simulation parameters and physical simulation time (sim.time) of the final equilibrated ab initio and force field MD simulations of [EMIm][OBz] and aqueous [EMTr][OBz] for the development and validation of the force field.