Nanoscale domains in ionic liquids: A statistical mechanics definition for molecular dynamics studies

TL;DR

This paper introduces a first-principles statistical mechanics definition of nanoscale domains in ionic liquids, using adaptive molecular resolution to identify atomistically self-contained nanodroplets, revealing heterogeneity scales.

Contribution

It provides a rigorous, physics-based method to define and analyze nanoscale domains in ionic liquids, advancing beyond arbitrary criteria.

Findings

Identification of atomistic nanodroplets in ILs

Quantitative measure of heterogeneity length scale

Application to four different imidazolium-based ILs

Abstract

One of the many open questions concerning Ionic Liquids (ILs) is the existence of nanoscale supramolecular domains which characterize the bulk. The hypothesis of their existence does not meet a general consensus since their definition seems to be based on ad hoc arbitrary criteria rather than on general and solid first principles of physics. In this work, we propose a suitable definition of supramolecular domains based on first principles of statistical mechanics. Such principles can be realized through the application of a recently developed computational tool which employs adaptive molecular resolution. The method can identify the smallest region of a liquid for which the atomistic details are strictly required, while the exterior plays the role of a generic structureless thermodynamic reservoir. We consider four different imidazolium-based ILs and show that indeed one can quantitatively represent the liquid as a collection of atomistically self-contained nanodroplets embedded in a generic thermodynamic bath. Such nanodroplets express a characteristic length scale for heterogeneity in ILs.

Tables

Table 1. Table 1: Radius of the minimal atomistic region of AdResS (size of the nanodroplet) for each system considered.

System	Radius (nm)
$\left[DMIM\right]\left[Cl\right]$	3.0
$\left[EMIM\right]\left[Cl\right]$	3.0
$\left[BMIM\right]\left[Cl\right]$	2.0
$\left[BMIM\right]\left[B{F}_4\right]$	3.0

Table 2. Table 2: Total dipole moment of the nanodroplet for each IL considered. It is reported also the systematic error obtained from AdResS and full atomistic simulations. The unit of μ 𝜇 \mu is in Debye.

System	AdResS	full atomistic
$\left[DMIM\right]\left[Cl\right]$	0.56$\pm$0.005	0.55$\pm$0.005
$\left[EMIM\right]\left[Cl\right]$	0.84$\pm$0.007	0.81$\pm$0.005
$\left[BMIM\right]\left[Cl\right]$	0.83$\pm$0.005	0.86$\pm$0.005
$\left[BMIM\right]\left[B{F}_4\right]$	0.98$\pm$0.008	1.02$\pm$0.009

Table 3. Table 3: Average number of ion pairs computed in the subregion of the AdResS and the same quantity computed in the equivalent subdomain of the full atomistic simulations

System	full atomistic	AdResS
$\left[DMIM\right]\left[Cl\right]$	582$\pm$2	594$\pm$2
$\left[EMIM\right]\left[Cl\right]$	565$\pm$2	554$\pm$2
$\left[BMIM\right]\left[Cl\right]$	137$\pm$1	130$\pm$1
$\left[BMIM\right]\left[B{F}_4\right]$	336$\pm$2	343$\pm$1

Equations

$$P({\mathbf{r}}_1........{\mathbf{r}}_N)\approx {\Pi }_{i,j}g(r_{ij})$$

$$F_{\alpha \beta }=w(X_{\alpha })w(X_{\alpha })F_{\alpha \beta }^{AT}+[1-w(X_{\alpha })w(X_{\alpha })]F_{\alpha \beta }^{CG}$$

Full text

{tocentry}

Dipolar nanodroplet model of ionic liquid embedded in a generic fluid

Ionic liquids (ILs) are an unusual class of liquids where the complex interplay of van der Waals, electrostatic and short-ranged solvophobic interactions determines the formation of alternating nano-metric structures in the bulk or at interfaces. The non-polar alkyl chains form non-polar domains by excluding themselves solvophobically from the ionic polar domains, in turn, this ability to spontaneously form self-assembled nano-structures without the need of external factors, allows them to play a key role in the evolution of modern nanotechnology. The list of properties of ILs which are used in current (nano)-technology is very broad, interested readers can consult the referenced literature 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12. Some relevant examples regarding Imidazolium based ILs, i.e. the class of systems considered in this work, are medicinal applications for the treatment of gram-negative and gram-positive bacteria’s, fungi, algae, human tumor cell lines13, 14. Importantly, it turns out that all these technological applications are often limited by the stability or dimensions of the nanoscale domains in bulk or at the interface. In fact, it has been suggested that the optimization of macroscopic properties such as heat capacity and compressibility is directly linked to the specific structure of such nanodomains. The natural consequence is that a rational design of the specific molecular chemistry for the optimization of a macroscopic property needs the knowledge about the ability of chemistry-specific ion pairs in building nanodomains with certain characteristics 15, 16.

As a consequence, the key question that chemical engineers must face concerns the determination of the optimum number of ions that account for the formation of stable apolar/polar nanoscale domains according to the specific chemical structure of the ions. The question can be actually reverted by asking how to specifically design a molecule (ion pair) so that a desired size of the nanoscale domain is achieved. Molecular design via computer simulation can play a major role in addressing such a question, in fact the chemical structure of a molecule can be modified in silico and then its liquid properties can be simulated at the thermodynamic conditions desired. The success of molecular simulation in explaining the important properties of ionic liquids and in promoting computational design of new species has been proven by a large number of publications (see Refs.9, 11, 8, 10, 17 and references therein). In literature, the definition of nanodomains in ILs is analyzed with a severe critical eye and some of the related results, such as the critical number of ions, are considered “magic” data, because associated to an arbitrary length scale which does not have proven physical or chemical basis 18, 16. In general one can define nanoscale domains such domains that arise due to solvophobic self-assembly and arising because of the amphiphilic nature of the ions; these can be thought of having an experimental meaning since their characteristic size is measured by x-ray scattering. In addition, there are small regions or aggregates which behave as if they hold one unit of charge (and thus giving a low ion concentration overall). This is usually considered only a physical model which is used to explain the observed screening length; these aggregates are not directly measured (or measurable). In this context, this work aims at contributing with an original analysis of possible nanodomains for a class of imidazolium-based ILs. Specifically, we aim at identifying the smallest region of the liquid whose properties strictly require the explicit consideration of atomistic degrees of freedom. In relation to the discussion above, both kinds of aggregates, i.e. those produced by the mechanism of solvophobic self-assembly and those which behave as a unit of charge, require (locally in space) the explicit chemical structure of the molecules. The extension in space for which the explicit consideration of the atomistic degrees of freedom (and thus charges) is mandatory would give information on the interplay between the two type of aggregates in terms of the interplay between entropy (molecular packing) and energy (direct atom-atom interactions). Such an interplay, in turn, can lead to a physical model of nanodomain that contains both kinds of aggregates, each with its specific signature within the definition of a nanodomain given by our analysis. It must be added that not all ILs are characterized by the heterogeneity of the first kind, only ILs with alkyl chain length $n\geq 4$ 16. For this reason we will consider ILs with alkyl chain length up to 4, so that we can treat both kinds of aggregates. The proposed analysis implies that in a subregion of the system the explicit chemical structure of the molecules is necessary for any calculated property, instead the rest of the system plays the role of a generic thermodynamic environment where explicit chemical details do not play any role.

Such an analysis is possible with a simulation tool developed in our group over the last few years, that is the Grand Canonical Adaptive Resolution Simulation method (GC-AdResS) 19, 20, 21, 22, 23, 24. Such a method allows for treating open subregions of a system at different level of resolution, specifically a high resolution region (atomistic) embedded in a large uncharged coarse-grained mean-field environment (see Fig1).

The method is based on a rigorous numerical implementation of Grand Canonical principles (for further details see Refs. 25, 26 and references therein); a relevant consequence for the current work is that the smallest atomistic region in GC-AdResS that can reproduce certain quantities of a full atomistic simulation, unambiguously defines a self-contained atomistic (sub)region (atomistic nanodroplet, w.r.t such quantities). In this perspective we will consider a series of imidazolium-based ILs, 1,3-dimethylimidazolium chloride (DMIM Cl), 1-ethyl-3-methylimidazolium chloride, (EMIM Cl), 1-butyl-3-methylimidazolium chloride (BMIM Cl) and 1-butyl-3-methylimidazolium tetrafluoroborate (BMIM BF4) (see Fig.2) and apply the GC-AdResS analysis to define the nanodroplet domain and analyze its physical and chemical characteristics.

Using state of the art force fields available in literature 27, 28, we will show that one can quantitatively describe the liquid as a collection of supra-molecular domains embedded in a chargeless particle-based mean-field. As a matter of fact, we provide an objective route to the definition of nanodomains with their associated physical properties. A relevant consequence of our analysis is the direct link between the microscopic scale of specific chemistry and mesoscopic structure of the liquid; in perspective, this link offers a possible route to a systematic rational design of large scale liquid properties.

0.1 Determination of the size of nanodroplets via GC-AdResS analysis

In order to determine the size of the atomistic nanodroplets we apply the GC-AdResS tool to the liquid. The application of such a tool implies that quantities calculated in the atomistically resolved region carry the full chemical characteristics of the molecules of such a region, while the influence of the exterior (coarse-grained region) enters only in the form of an action of a generic bath that thermodynamically equilibrates the atomistic region. The implication is that if a certain quantity is calculated in the atomistic region of GC-AdResS and displays the same behavior when calculated in the equivalent subregion of a full atomistic simulation, then the explicit role of atomistic degrees of freedom outside such region is negligible for the quantity considered. As a consequence, with respect to such a quantity, the atomistic region of GC-AdResS can be considered as an atomistically self-contained nanodroplet embedded in a liquid of uncharged coarse-grained spheres( one for the anion and one, or more bonded, for the cation). In essence, the exterior represents a generic fluid, i.e. without any information regarding the chemical structure of the molecules, that operates as a large thermodynamic bath for the atomistically resolved droplet. The quantification of the minimal size of such subdomain is determined by running a sequence of GC-AdResS simulations with decreasing size of the atomistic region. The smallest atomistic region for which the target property has the same behavior as in a full atomistic simulation (within maximum $5\%$ of discrepancy) defines the size of the nanodroplet. The choice of the quantity to use as a target for the definition of an atomistic nanodroplet plays a key role; it must be such that on the one hand it carries a specific signature of the system and on the other it should be a quantity that can be directly used to determine a broad range of physical properties. In previous work 23 we have proposed that atomistic radial distribution functions are ideal quantities to use as target. In fact they uniquely characterize the liquid structure of a substance and in addition, at statistical mechanics level, the probability distribution function of an atomistic system can be approximated as a factorization of (all possible) atom-atom radial distribution functions (two-body correlation functions):

[TABLE]

where ${\bf r}_{1}........{\bf r}_{N}$ is the spatial configuration of the $N$ molecules of the liquid, $i$ and $j$ labels atoms belonging to different molecules at a distance $r_{ij}$. Since the knowledge of $P({\bf r}_{1}........{\bf r}_{N})$ implies the knowledge of every average property/quantity of the system, its approximation via the factorization proposed above implies the complete statistical knowledge of the system up to the two-body level (second order). The determination of the size of the nanodroplets using GC-AdResS is done by determining the smallest atomistic region for which every possible atom-atom radial distribution function has the same behavior as the corresponding quantity calculated in the equivalent subregion of a full atomistic simulation. In addition, being the atomistic region an open system governed by Grand Canonical principles, it is mandatory to require that the molecular number probability, $P(N)$, must have the same behavior (within a certain accuracy) of the , $P(N)$, calculated in the equivalent subregion of a full atomistic simulation of reference. In molecular simulation, for a system at uniform resolution (i.e. full atomistic), it is trivial to obtain all the atomistic radial distribution functions, provided that the size is sufficiently large to justify the statistical average. In our specific case the situation is very different: outside the atomistic region the atomistic structure of the liquid does not exist, thus it is not obvious to reproduce the atom-atom radial distribution functions in the atomistic region of AdResS. The agreement is possible only if the atomistic information of the bulk is not required for the local structure-based microscopic properties. We must clarify that here the dynamics is not treated, and that we focus only on static properties. This criterion has been proven to be unambiguous and effective in identifying relevant length scale at nanoscopic level in previous studies of relevant systems 21, 29, 30, and in preliminary work on ILs 31, 32, 33. At this point a warning must be added, that is, Eq.1 implicitly means that properties where three-body (and higher) accuracy is required cannot be treated with the current approach, unless the corresponding n-th-body correlation function is used as a target. This means that the extension to higher accuracy and further statistical properties, is technically possible

In the next section we report the calculations done for each system and we quantitatively define the nano-metric characteristics of each liquid. Furthermore, as a relevant physical property of such nanodomains, we calculate the total dipole moment of the atomistic region of GC-AdResS and compare it with the equivalent quantity calculated in the same sub-domain of a full atomistic simulation and prove that indeed the values agree within a maximum discrepancy of $5\%$. This calculation is of particular interest, in fact implies the existence, within the defined nanodomains, of smaller domains with a predominance of positive and negative ions.

1 Results

The technical details of the simulation are reported in the methods section (for further details, see 31, 32, 33). Here we report the results of the atom-atom radial distribution functions and molecular number probability. It must be underlined once again that we consider all possible atom-atom radial distribution functions (according to the symmetry of the molecules), however for clarity of presentation and in order to avoid an excess of equivalent information Figs.(3-6) report only some representative examples of the atom-atom radial distribution functions calculated in GC-AdResS and compared with the same quantities calculated in the equivalent subregion of a full atomistic simulation. All the other functions are not reported explicitly, but they are characterized by the same accuracy regarding the agreement between the AdResS simulations and the full atomistic simulation of reference.

The size of the atomistic domain varies between $2.0$ and $3.0$ nm for the different ILs considered (see table 1; as made clear before, outside the atom-atom radial distribution functions do not exist, thus if the atomistic data of the exterior were important for the radial distributions in the atomistic region, then one must have seen a clear perturbation of the functions.

As anticipated above, the atomistic region of GC-AdResS is an open system, thus the definition of a statistically self-contained atomistic region would be justified if and only if the exchange of particles with the exterior occurs in the same way as in an equivalent subregion of a full atomistic simulation. Such a behavior is expressed by $P(N)$ i.e. the particle number probability,

Fig.7 shows that for every system considered the atomistic domain identified with GC-AdResS display a very good agreement with the full atomistic simulation of reference regarding $P(N)$ (within $3\%$ accuracy). In the light of the results reported above, one can be certain that the representation of the exact many-body probability distribution function of the atomistic domain is correct up to the second order. As a consequence any quantity calculated as an ensemble average in the atomistic region has the same value when calculated in the equivalent subregion of a full atomistic simulation of reference, i.e. we can identify local characteristics of the liquid in space by considering only one atomistic droplet defined according to the criteria listed above. In this perspective, our proposed definition of a nanodomain can be considered very solid from the point of view of statistical mechanics.

1.1 Dipole Moment of the Nanodroplet

We can now proceed to further explore the physical characteristics of the nanodroplet defined by GC-AdResS. Being ILs strongly characterized as a mixture of positive and negative ions, a particularly interesting quantity to look at is the (ensemble) average dipole moment of the droplet. Although ILs display spatial heterogeneity at certain scale (see e.g. 34, 35) in general one expects that on larger spatial scales there is a certain homogeneity, that is the nanodroplets are expected to be neutral (in average). However, within the single nanodroplet there may exist spatial heterogeneity, e.g. regions with higher concentration of anions (cations) thus the nanodroplets are characterized by an electric dipole moment. Indeed we have found that a picture of this kind can be drawn from our results. Figs.8, show the comparison between the AdResS results and the same quantity calculated in the corresponding sphere in a full atomistic simulation. The agreement is highly satisfactory and is within a minimum accuracy of $5\%$ (see also table 2).

Interestingly the existence of anion-rich/cation-rich domains at smaller scale, is consistent with results and conclusion of previous work 9, 11, 8, 10 and once again it strengthen the idea of ILs as substances characterized by scales with different spatial heterogeneity at mesoscopic level 34, 35. In perspective, one may think of a model for this class of ILs as a collections of dipolar nanodroplets embedded in a generic (mean-field) fluid, as sketched in Fig.9, however, while it may represent an interesting research plan, at this stage it goes beyond the scope of the current work.

Another interesting point is that the picture of the liquid as a collection of dipolar nanodroplets may strengthen the complementarity in the interpretation of relevant experimental results. In fact, Israelachvili et al. 36 have shown that ionic liquids behave as weak diluted electrolytes where the electrostatic energy between ions decay exponentially as a function of distance between them, that is, the Debye screening or decay length monotonically decreases even in the limit of high ion concentration because of the effective screening of charges at shorter separations. In other words, long-range forces are expected to not play a major role in concentrated electrolytes such as ionic liquids (when we use the Debye screening model). At the same time, the group of Perkin 37, 38 have shown that the conventional of strong screening in ILs does not hold in the limit of high charge concentration of ions. Using simple scaling theory they show that there is an anomalous dependence of the decay length on the ion concentration of an order of magnitude larger than the ion diameter in ILs. Our results allow for an interpretation that enforces the complementarity between these two views. In fact, while the local structure of the droplet does not explicitly depend on the electrostatic long range interactions (our coarse-grained model does not have charges), however each nanodroplet is characterized by a dipole, that is the larger spatial scale structure of the liquid (i.e. the structure of the liquid as a collection of nanodroplets) instead will be influenced by the dipole-dipole interaction of the nanodroplets (consistently with other proposals, see Refs. 39, 36, 40, 38, 41). Moreover, it is important to note that the long screening length measured in Ref.38 is in accordance with inflated screening lengths in highly concentrated simple electrolytes such as 3M NaCl in water, thus it demonstrates that the details of molecular/ionic structures are probably not important for determining the long screening length. This conclusion can be drawn also from our study where beyond the local scale of the nanodroplet the specific chemical structure of the molecules does not play anymore a role, in fact we can use generic (and uncharged) spherical coarse-grained models. Finally, another important aspect related to the dipole of the nanodroplet is that the unambiguous locality proven by AdResS is also qualitatively supported by electronic structure-based related research. In fact the vibrational spectra calculated from quantum simulations of small samples are in very good agreement with experiments. One can qualitatively conclude that indeed there may exist subregions of the liquid which do not need the explicit environment for the calculation of certain microscopic properties42.

2 Conclusion

We have employed molecular dynamics within the GC-AdResS scheme and used a statistical mechanics/probabilistic criterion to define and to quantify the size of nanodomain in a class of imidazolium-based ILs. The main conclusion concerns the characterization of the ILs as liquid formed by mesoscopic nanodroplets, a picture consistent and complementary to the results obtained in other studies. Thus, in conclusion, our aim is to enforce the idea, based on previous excellent work of other groups, that molecular simulation can indeed be employed to defined nanodomains in ILs. As a practical consequence, such studies can certainly help with a detailed in silico rational design of ions for macroscopic properties on the liquid on demand.

3 Methods

Molecular Dynamics has been used for the simulation of the systems considered. In particular, within such an approach, we have used the Adaptive resolution simulation method (AdResS) 43, 22 in its Grand Canonical (GC) set up 44, 23, 24, 26, 45. The standard Molecular Dynamics method is well established and does not need further explanation, instead the AdResS method is not as popular and thus we report here its basic features that complement the description provided before in the text. A detailed description can be found in the references cited above.

In the Adaptive Resolution Simulation technique, space is divided in three regions, CG, HY and AT, characterized by coarse-grained, hybrid and atomistic resolution respectively and molecules adapt their structural resolution according to the region in which they are instantaneously located. The algorithm is based on a space dependent interpolation for the force between two molecules, $\alpha,\beta$ (see Fig.1):

[TABLE]

$F_{\alpha\beta}^{AT}$ is the atomistic force while $F_{\alpha\beta}^{CG}$ is the coarse-grained force. $w(x)$ is a smooth function that goes from 0 to 1 in a transition region $\Delta$ (with $x$ center of mass of the molecule), a thermostat assures a target thermodynamic equilibrium. An additional force, ${\bf F}_{th}(x)$ (thermodynamic force), derived from first principles of statistical mechanics, acts on the center of mass of the molecule in $\Delta$. Such a force gives thermodynamic consistency to the method, i.e. assures that the effective chemical potential of the whole system corresponds to that of the (target) atomistic resolution. With this set up, the approach becomes numerically robust and conceptual rigorous for the simulation of system with open boundaries, thus in case the coarse-grained region is large enough, the method describes a Grand Canonical system for the atomistic region (GC-AdResS). In practice, ${\bf F}_{th}(x)$ is calculated during a simulation (equilibration run), as the gradient of the number density of the molecules, in an iterative form 25: $F_{k+1}^{th}(x)=F_{k}^{th}(x)-\frac{M_{\alpha}}{[\rho_{ref}]^{2}\kappa}\nabla\rho_{k}(x)$, $M_{\alpha}$ is the mass of the molecule, $\kappa$ a tunable constant, $\rho_{k}(x)$ is the molecular density in $\Delta$, at the $k$-th iteration and $\rho_{ref}$ is the target density of reference. The criterion of convergence depends on the accuracy required for the simulation, in general, $|\rho_{final}-\rho_{ref}|$ must be always below $10\%$ in $\Delta$. The GC-AdResS method has been implemented in the Molecular Dynamics package GROMACS 46, below are reported the specific details of the simulations of this work.

For each of the ILs, we have set up two simulations: 1. a full atomistic simulation which we use it as a reference to compare the statistical properties calculated in the sub-atomistic region of AdResS with the equivalent subregion of the full atomistic simulations. 2. GC-AdResS simulation with spherical atomistic resolution. All simulations were optimized using NPT ensemble. The temperature and time step corresponds to 400K and 2fs respectively. The particle mesh Ewald (PME) is used to determine the coulombic interactions. For the first 5ns we use Berendsen barostat47 which is then switched to Parrinello-Rahman48 barostat for the remaining length of the simulation. The equilibration is carried over until the fluctuating box length reach an order of 0.0001nm. This configuration is used as the initial configuration for the production runs in NVT ensemble.All the statistical quantities were computed during the production run. The production run length was 10ns for [DMIM][Cl] and 25-40ns for [BMIM][BF4]. The force field parameters, the step-by-step procedure to obtain the neutral coarse grained potentials and the thermodynamic forces for the coarse grained beads for the four different IL systems are discussed elsewhere 31, 32, 33. We use the tabulated CG potentials and thermodynamic force data to setup the GC-AdResS simulations. The GC-AdResS simulation with spherical resolution of the atomistic region for two of the IL systems ([DMIM][Cl] and [BMIM][BF4]) are discussed here and the reminder systems are discussed elsewhere31, 32. The radius of the atomistic region in the GC-AdResS simulation is fixed at 3.0nm for both [DMIM][Cl] and [BMIM][BF4] (see table 1), which turned to be the minimal size of the atomistic region that reproduces all the possible atom-atom radial distribution functions. The hybrid region that borders the sub atomic region has a width of 2nm for each system. The coarse grained region spreads out throughout the reminder of the cubic box of length 15.69333Å (20000 ion pairs [DMIM][Cl]) and 14.17269(8451 ion pairs [BMIM][BF4]) respectively. An optimized $\Gamma=5ps^{-1}$ parameter is used with Langevin thermostat. Reaction field method is used to calculate the long ranged electrostatic interaction between particles49. The quality of the thermodynamic force is such that the density in the atomistic subregion of the GC-AdResS is comparable to the bulk atomistic system with a maximum fluctuations $<5\%$ for both the systems.

4 Acknowledgment

This research has been funded by Deutsche Forschungsgemeinschaft (DFG) through the grant CRC 1114: “Scaling Cascades in Complex Systems”, project C01. We would like to thank Prof.Barbara Kirchner and Prof.Susan Perkin for relevant suggestions that shaped the work.