Ionic Liquids Treated within the Grand Canonical Adaptive Resolution Molecular Dynamics Technique

Abstract: We use the Grand Canonical Adaptive Resolution Molecular Dynamics Technique (GC-AdResS) to examine the essential degrees of freedom necessary for reproducing the structural properties of the imidazolium class of ionic liquids. In this technique, the atomistic details are treated as an open sub-region of the system while the surrounding environment is modelled as a generic coarse-grained model. We systematically characterize the spatial quantities such as intramolecular, intermolecular radial distribution functions, other structural and orientational properties of ILs. The spatial quantities computed in an open sub-region of the system are in excellent agreement with the equivalent quantities calculated in a full atomistic simulation, suggesting that the atomistic degrees of freedom outside the sub-region are negligible. The size of the sub-region considered in this study is 2 nm, which is essentially the size of a few ions. Insight from the study suggests that a higher degree of spatial locality seems to play a crucial role in characterizing the properties of imidazolium based ionic liquids.


Introduction
Ionic liquids are a new class of chemical compounds that are highly tunable in nature with unusual properties such as high electrical conductivities, improved electrical and thermal stabilities, good solvation ability, negligible vapor pressures, non-flammability, to mention a few [1][2][3].By chemical modification of the polar or apolar components or charge of the ions, numerous ILs with a range of physicochemical properties such as viscosity, catalytic activity, or solvation can be generated.We now have task specific ionic liquids that can be applied for a range of designated applications [4,5]: ionic liquids as drugs [6,7], as batteries and super-capacitors for solving energy demands of the environment [8][9][10][11], to capture greenhouse gases, to dissolve proteins [12,13] and cellulose [14,15].
The physicochemical properties of ionic liquids can be understood by studying the morphological structure of these novel materials.Therefore, a plethora of experiments and simulations have been performed to probe the structural morphology of ionic liquids with an emphasis on pair correlation functions (PCFs) [16][17][18][19].The conventional approach to probe the structure of ionic liquids using scattering experiments involves determination of structure factor that can be Fourier transformed to obtain pair correlation functions (PCFs) or radial distribution functions (RDFs) [20].Atomic-level organization of ILs can be obtained by properly analyzing the structure factor and corresponding pair correlation functions including large scale correlation functions [21] of various components of ionic liquids.However, in our previous studies on ILs, we have shown that a large class of ILs can be characterized by a local scale where important properties that may be expected to be dominated by large scale pair correlations, are instead highly localized and depend only on the immediate neighboring molecules [22][23][24].Motivated by these results, the adaptive resolution scheme (AdResS) [25,26] has been extended to grand canonical-like version [27][28][29][30] for the study of ILs.
In fact, AdResS simulation allows one to obtain full atomistic details of a region of interest and couple it with a coarse-grained region modeled without any atomistic details.The generic coarse-grained region acts as a reservoir of particles and energy, which suggests that GC-AdResS can be used as a tool to identify the essential atomistic degrees of freedom required to have a certain property.This idea already been used in previous AdResS studies to infer the locality of ionic liquids [31,32], IR spectrum in water [33] and solvation properties for hydrophobic molecules in water [34].In this perspective, the aim of this work is to show the technical feasibility of an AdResS study of ILs consisting of cations with longer alkyl chains and relatively larger anions compared to our earlier study [31,32].In this work, we have chosen [BMI M][Cl] and [BMI M][BF 4 ] as the test systems for our simulation not only for its inherent complexity but because these ILs have been shown to be useful in the development of breakthrough application in bio-nanotechnology [35,36].

Results
The analysis of GC-AdResS simulation firstly involves whether the necessary thermodynamic conditions for physical consistency such as the structural and orientational properties of the ILs are fulfilled.

Structural Properties of ILs
In Figure 1 we show the density of ion pairs for the [BMIM][Cl] and [BMIM][BF 4 ] system as a function of the distance from the center of the atomistic region.As expected, the largest deviation in density occurs in the hybrid region with an upper bound of about 10% and 7% for [BMIM][BF 4 ] and [BMIM][Cl] respectively.However, within the atomistic region, which is the region of interest in our study, the largest deviation observed is ≤5 for both the systems.Given the critical technical conditions, in addition to the inherent complexity of the ILs (for example: viscosity, η of [BMIM][BF 4 ] and [BMIM][Cl] at 298 K are 101 cP and 29 cP respectively [37]) we consider the density difference to be satisfactory for studying the structural properties of ILs.The fluctuations in the density can be further reduced by improving the convergence of the thermodynamic force.
In order to understand the behavior of the particle number fluctuation (ie) the exchange of matter with the large reservoir, we plot the diffusive number density profile for the ion pairs located in the AT, HY and CG regions at various time frames for the [BMI M][Cl] and [BMI M][BF 4 ] system (see Figures 2 and 3).The results show that there is diffusion of particles from each region to the other suggesting that the GC-AdResS method behaves reasonably well.

Orientational Properties of ILs
In addition to the structural properties, we compute the orientational alignment of the molecules to ensure the possibility of any artificial molecular orientation due to the interface introduced by the different spacial resolution of the coarse grained models.The orientational order parameter is defined as |cosθ| = | v |v| .x|,where v/|v| is the unit vector along the principal axis of symmetry of the cation and x is the unit vector defined between the cation and the center of the atomistic region [31].For a random rotation of the molecules, the value of cosθ must be 0.5.The probability distribution of the average order parameter, |cosθ| measured within the sub-regions of the AdResS simulations are compared with the corresponding quantity in the equivalent sub-region of the full atomistic simulation (see Figure 4).The average value of cosθ obtained is 0.5 for the AdResS simulations compared to the reference system, suggesting that the rotation of the molecules are random.full atomistic AdResS:q(CG)=+/-0.8 AdResS:q(CG)=+/-0

Locality of Structural Properties
After assuring the thermodynamic and statistical behavior of the atomistic region, we studied the essential degrees of freedom necessary for reproducing the local structural properties of ILs.The intermolecular (anion-anion, anion-cation, cation-cation) and intramolecular (H-H, H-Cl, C-C) radial distribution functions are shown in Figure 5 for [BMI M][Cl] and in Figure 6 for [BMI M][BF 4 ] system respectively.The center of mass of the cations and anions are represented as CAT, ANN or Cl respectively.CAA is the carbon, of the methyl group attached to the imidazolium ring and H is the hydrogen that is present in between the 1-methyl and 3-butyl chain of the cations.The intermolecular and intramolecular RDFs are calculated in the atomistic region (2 nm) and compared with the equivalent functions in the corresponding sub-region of the full atomistic system.The radial distribution functions represent a second order probability distribution function of the system, as it is reduced from a 3N-probability distribution function to a factorized function of independent two-body terms (see also [28]).Therefore, reproducing the radial distribution functions of the atomistic region of GC-AdResS reproduce the behavior of the same quantities calculated in the sub region of a full atomistic simulation, suggesting that the two sub-regions are statistically equivalent at least up to the second order.Our results for the atom-atom radial distribution functions in the AT region of AdResS reproduce the results of the reference full atomistic simulation, which implies that the only relevant atomistic degrees of freedom are those of the AT region of 2 nm radius.Thus from these local radial distribution functions we can further infer that the ensemble average of quantities, depending on spatial coordinates, can be localized in the sub-region with an accuracy at least up to the second order in the probability distribution.

Method
In this study we employ the Grand Canonical Adaptive Resolution Simulation Scheme (GC-AdResS) for molecular dynamics [25,26,[28][29][30]38].Two regions of a simulation box having different resolution regions: 1. atomistic (AT) and 2. a coarse-grained (CG) system can be linked using this multi-resolution simulation method.Free exchange of particles between the two resolution regions are ensured through a coupling transition region (HY) where molecules have space-dependent hybrid atomistic/coarse-grained resolution.The force between two molecules α and β present in the hybrid region is computed via a space-dependent interpolation formula and can be written as here, the force between the particles derived using the atomistic and CG potentials are F AT αβ and F CG αβ respectively.The interpolating switch w(x) is defined as: In the above equation, the spatial extent of the atomistic and hybrid regions are represented as d AT and d ∆ respectively (see Figure 7).The yellow line represents the weighting function that allows the coarse grained particles to change their resolution into an atomistic molecule as it varies smoothly from 0 to 1 and vice versa.Earlier studies have shown that this simulation setup is sufficient for performing accurate simulations in both the atomistic and coarse-grained region [39].
In order to ensure proper thermodynamic equilibrium and correct particles' exchange between he AT and CG region [40,41], thermodynamic force is introduced in the transition region.Studies have shown that models that can reproduce macroscopic thermodynamic quantities such as density and temperature can be used as coarse-grained models, so that AdResS is re-framed within the Grand Ensemble idea of open boundary systems [29,30,38].In this work, we use a spherical atomistic region of radius d at that is surrounded by a hybrid region of width d hy .The remainder of the box is treated as the CG region that serves as a particle reservoir for the AT region.The force parameters of the atomistic and coarse-grained model are discussed below in Sections 3.2 and 3.3.with the radius of the atomistic region set to 2 nm.The CG, HY and AT represents the coarse-grained region, hybrid region and the atomistic region respectively.The forces in the AT and CG regions are derived from the respective potentials.In the HY region, a space dependent interpolation scheme with a weighting function w(x) (shown in yellow) is used for obtaining the forces between two particles.The cation is modelled as two coarse grained beads.One sphere for the imidazolium ring and the other sphere for the alkyl side chain.Similarly, the anion is modelled as one coarse grained bead (bottom column).

Models
The all atom forcefield parameters to model imidazolium ionic liquids that has long alkyl chains attached to the cation ([BMIM]) and relatively hydrophobic anions (Cl − , BF − 4 ) were taken from the Dommert et al. [42] and Kirchner, B. et al. [43].The model considers reduced/scaled charges of +0.8 e and −0.8 e instead of the full formal charges +1.0 e and −1.0 e.Within the CG region, ionic liquids are modelled as charged/neutral spherical atoms.We use two spherical beads to model the cation (in which the imidazolium ring and the alkyl chain are modeled using separate spherical beads) and one bead to represent the anion.We use Inverse Boltzmann Iterative procedure (IBI) to develop the coarse-grained interaction parameters.This procedure takes into account the corresponding charge in the coarse-grained model during the iterative process.

Technical Details
[BMI M][Cl] − : We performed all the MD simulations using GROMACS package [44].We set up two systems, a relatively smaller system consisting of 400 ion pairs and a larger system with 19,968 ion pairs.The full atomistic 400 ion pairs simulation was used to generate coarse-grained models, which was used in GC-AdResS simulations.However, the full atomistic 19,968 ion pairs simulation was used as a reference system to compare the statistical quantities computed within the sub-region of the GC-AdResS simulations with that computed in the equivalent subregion of the atomistic simulations.Both the systems were optimized using full atomistic NPT calculations.The temperature of the system was set to 400 K with a simulation time step of 2 fs.The particle mesh Ewald (PME) technique was used to calculate the electrostatic interactions between ions.Berendsen barostat [45] was used for the first 5 ns and then Parrinello-Rahman [46] barostat was used for another 5 ns.The box size is considered to be converged, when the changes in the box length were of the order of 0.0001 nm.For 400 and 19,968 ion pairs we obtained a cubic box with 4.31080 nm and 15.88645 nm respectively.The production runs are performed using NVT ensemble.The radial distribution function was obtained after 5 ns full atomistic NVT simulations at 380 K and with 2 fs timesteps.The two coarse-grained models were derived from the 400 ion pairs configuration.Inverse Boltzmann iteration (IBI) procedure [47] was used to reproduce the radial distribution functions of the full atomistic target systems (for neutral and charged model see Figure 8).We used the tabulated coarse grained potentials to set up 19,968 ion pairs GC-AdResS system.We fix the resolution of the atomistic region in the GC-AdResS simulation to 2 nm, which is sufficient to capture the structural properties of the system similar to the reference atomistic system.The width of the hybrid region that borders the atomistic region is set to 2 nm, and the remaining length corresponds to the CG region.A Langevin thermostat is used with Γ = 5 ps −1 .Generalized reaction field method with a self-consistent dielectric constant was used to compute electrostatic interactions between ions [48].In the case of charged coarse-grained model, the reaction field method is applied in the coarse-graining region similar to the atomistic region.The thermodynamic force through the iterative procedure was considered converged after 100 iterations, sufficient to reach an accuracy of <5 particle density and for the ion-ion radial distribution functions ( compared to the full atomistic of the reference).
[BMI M][BF 4 ] − : We set up two systems, a relatively smaller system consisting of 313 ion pairs and a larger system with 939 ion pairs.Both the systems were optimized using full atomistic NPT calculations.The temperature of the system was set to 400 K with a simulation time step of 2 fs.The particle mesh Ewald (PME) technique was used to calculate the electrostatic interactions between ions.Berendsen barostat [45] was used for the first 1 ns and then Parrinello-Rahman [46] barostat was used for another 3 ns.The box size is considered to be converged, when the changes in the box length were of the order of 0.0001 nm.For 313 and 939 ion pairs, we obtained a cubic and orthogonal box of 4.72228 nm and 14.16684 nm × 4.72228 nm × 4.72228 nm respectively.The production runs are performed using NVT ensemble.The radial distribution function was obtained after 10 ns full atomistic NVT simulations at 400 K and with 2 fs timesteps.The two coarse-grained models were derived from the 313 ion pairs configuration.Inverse Boltzmann iteration, IBI, procedure [47] was used to reproduce the radial distribution functions of the full atomistic target systems (for charged model see Figure 9).We used the tabulated coarse grained potentials to set up 939 ion pairs GC-AdResS system.We fix the resolution of the atomistic region in the GC-AdResS simulation to 2 nm, which is sufficient to capture the structural properties of the system similar to the reference atomistic system.The width of the hybrid region that borders the atomistic region is set to 2 nm, and the remaining length corresponds to the CG region.A Langevin thermostat is used with Γ = 5 ps −1 .Generalized reaction field method with a self-consistent dielectric constant was used to compute electrostatic interactions between ions [48].In the case of charged coarse-grained model, the reaction field method is applied in the coarse-graining region similar to the atomistic region.The thermodynamic force through the iterative procedure was considered converged after 30 iterations, sufficient to reach an accuracy of <5% particle density and for the ion-ion radial distribution functions (compared to the full atomistic of the reference).

Conclusions
We have tested the spatial locality of complex ILs consisting of cations with long alkyl chains and hydrophobic anions using the Grand Canonical Adaptive Resolution Molecular Dynamics technique.By calculating various radial distribution functions between anions and cations within an atomistic region of radius 2 nm and comparing them with the equivalent quantities in a full atomistic simulation, we could conclude that even for a spherical region of radius 2 nm, the atomistic degrees of freedom of the bulk do not play a relevant role.Although, the spatial localization has been observed in complex liquids characterized by strong local anisotropic interactions such as: class of imidazolium based ionic liquids [31,32], water [33], the spherical region of radius 2 nm is system specific and could vary depending upon the complexity of the system of study.The insight from this study can be used for designing ILs at the molecular level in the framework of structure-function relations.

Figure 2 .Figure 3 .
Figure 2. Evolution in time of an instantaneous diffusive number density profile along the trajectory for the ion pairs that are at time, t = 0, located in the atomistic region (top panel), hybrid region (middle panel) and in the coarse-grained region (bottom panel).GC-AdResS simulation with the charged and uncharged coarse-grained model of the [BMIM][Cl] system is shown in left and right respectively.Exchange of ion pairs among different regions are shown.

Figure 4 .
Figure 4.The probability distribution of the averaged order parameter |cosθ| in the atomistic region of the 2 nm AdResS simulations compared with the corresponding quantity in the equivalent sub-region of the full atomistic simulation for the [BMI M][Cl] −1 (left) and [BMI M][BF 4 ] −1 (right) system.

Figure 5 .Figure 6 .
Figure5.The variation of intermolecular (left) and intramolecular (right) radial distribution profile calculated in the atomistic region of AdResS for (charged/neutral) coarse-grained models and compared with the corresponding quantity calculated in the equivalent region of a full atomistic simulation the [BMI M][Cl] − system.The center of mass of the cations and anions are represented as CAT, ANN or Cl respectively.CAA is the carbon, of the methyl group attached to the imidazolium ring and H is the hydrogen that is present in between the 1-methyl and 3-butyl chain of the cations.

Figure 7 .
Figure 7. GC-AdResS scheme is shown in the (top column); Simulation snapshot of the [BMIM][Cl] systemwith the radius of the atomistic region set to 2 nm.The CG, HY and AT represents the coarse-grained region, hybrid region and the atomistic region respectively.The forces in the AT and CG regions are derived from the respective potentials.In the HY region, a space dependent interpolation scheme with a weighting function w(x) (shown in yellow) is used for obtaining the forces between two particles.The cation is modelled as two coarse grained beads.One sphere for the imidazolium ring and the other sphere for the alkyl side chain.Similarly, the anion is modelled as one coarse grained bead (bottom column).

1 Figure 8 .
Figure 8.The radial distribution function obtained from reference atomistic and coarse-grained simulation with charged (bottom) and neutral model as a function of distance r for the [BMI M][Cl] − .The CG simulation employs the IBI potential, V CG (r)(kJ/mol −1 ) for the charged and uncharged CG model (inset).

8 Figure 9 .
Figure 9.The radial distribution function obtained from reference atomistic and coarse-grained with charged (bottom) and neutral model (top) as a function of distance r for the [BMI M][BF 4 ] − .The CG simulation employs the IBI potential, V CG (r)(kJ/mol −1 ) for the charged and neutral CG models.(inset).