Perspectives in the Computational Modeling of New Generation, Biocompatible Ionic Liquids.

In this Perspective, I review the current state of computational simulations on ionic liquids with an emphasis on the recent biocompatible variants. These materials are used here as an example of relatively complex systems that highlights the limits of some of the approaches commonly used to study their structure and dynamics. The source of these limits consists of the coexistence of nontrivial electrostatic, many-body quantum effects, strong hydrogen bonds, and chemical processes affecting the mutual protonation state of the constituent molecular ions. I also provide examples on how it is possible to overcome these problems using suitable simulation paradigms and recently improved techniques that, I expect, will be gradually introduced in the state-of-the-art of computational simulations of ionic liquids.

Introduction

Few inventions have impacted society and especially societal relations as the computer has done. The recent diffusion of easily portable computers in the form of smartphones and other devices has altered the way in which people communicate and has made possible the emergence of a new form of society (unthinkable just few decades ago) where the dissemination of information is immediate and global. The same revolution has impacted science and the way in which it is communicated or divulged. A less evident (at least outside the scientific community), but likewise innovative revolution has impacted the way in which science is done: a revolution boosted by the impressive surge in computer power and by the ever-increasing availability of computers.

The importance of computational approaches in science can be traced back few decades by just highlighting some of the past Nobel prizes. It starts with Fukui and Mulliken and the work of Hauptman and Karle and continues in 1998 with the prize to Kohn and Pople and, even more recently, with the recognition of the pioneering work of Karplus, Levitt, and Warshel in 2013. The 2021 Nobel prize has been awarded to Manabe and Hasselmann for “for the physical modelling of Earth’s climate” and to Parisi, who was involved in the birth of APE supercomputers in early nineties, for his work on complex systems. Other important achievements in human knowledge have been reached through the intensive use of computer and computations, just to name two: the measurement of the rate of the expansion of the Universe and the mapping of the human genome.

As implicitly recognized by the short and incomplete list above, chemistry is one of the fields that has benefitted most from the introduction of computers and their use. The application of computations to chemistry has spawned an entire new discipline aptly named “computational chemistry”. Under this umbrella we find the widest field of activities ranging from biology to molecular physics. This means that computational chemistry concerns a huge range of different systems, from the micro to the nano scale, from millions of atoms (e.g., proteins) to few atom systems (isolated molecules). The study of molecular systems from a fundamental point of view through computer models, has consequences in high societal impact areas such as climate changes, green economy, medicine, etc.

Despite the huge improvements in computer performances and in methodologies, the modeling and simulation of complex and heterogeneous molecular systems is still a challenging, often frustrating, task in chemistry. Notwithstanding, the opportunity to provide reliable data using a computer model has been and will be a great advantage for an efficient design of laboratory activities. It is sufficient to think that a well-designed set of computational simulations can avoid tedious and repetitive experimental procedures, thereby allowing the research to focus on the most promising materials or processes for a given aim. The realization of this scenario obviously depends on the reliability of the computer models, a question that is now central to the computational chemistry community and that is the subject of this article.

Computational chemistry consists of a wide array of techniques and methods which have their roots and foundation in physics, in the molecular Schrodinger equation (SE, mainly in its nonrelativistic form). In principle, the SE tells us how a given physical quantum system evolves in time. As is well-known, solving the SE for the quantum evolution of a molecular system of a decent complexity (let us say more than five atoms) is impractical or simply impossible. Computational chemistry, using Dirac’s words, consists of the development and use of “approximate practical methods of applying quantum mechanics [···] which can lead to an explanation of the main features of complex atomic systems without too much computation”.

Two of these approximate methods are the cornerstone of modern computational chemistry: the Born–Oppenheimer (BO) approximation and the classical motion of the nuclei. The former allows the calculation of electronic energies and the use of these as the source of the forces that act on the nuclei, the latter permits to blend nuclei into “classical” atoms for the implementation of simple and inexpensive methods to compute the evolution of complex chemical systems.

The BO approximation can be formulated for all molecular systems and states that the nuclei move under the effect of forces that are the gradient of the electronic energies computed at fixed nuclear geometries. This, effectively, decouples the solution of the nuclear and electronic SE in a two-step procedure: (i) first, one computes the electronic energy as a function of the nuclear geometries, and (ii) then, one solves the nuclear motion using the electronic energy as the potential. The computational techniques for the first step are well-known and are based on either wave function-based methods (Hartree–Fock or post-Hartree–Fock calculations)[1,2] or Density Functional Theory (DFT).[3] Both approaches are often called ab initio methods because they do not employ empirical parameters.

Essentially, operating under the BO approximations consists of assuming that the electrons are instantaneously in a stationary state at a given nuclear geometry and that the couplings due to the simultaneous motion of nuclei and electrons are negligible. In practice, this is achieved by initially fixing the position of the nuclei (an approach known as clamped nuclei) and solving the resulting electronic SE where the nuclear kinetic energy has been neglected. Such SE is solved through ab initio approaches to find an approximation to the ground state wave function (or energy) that governs the electronic motion subjected to the specific, geometry-dependent nuclear electrostatic potential. The ground state electronic energy (i.e., the eigenvalue of the electronic SE) depends on the nuclear geometry and is a function of the nuclear coordinates. As such, in the BO approximation, the electronic energy is interpreted as the potential (i.e., the source of the force) that governs the nuclear motion.

Chemistry is essentially the study of transformations that take place over a given time. It therefore follows that dynamic, i.e. the nuclear (atomic) motion, is the crucial issue. Electronic structure calculations such as DFT are inherently static, that is they provide an instantaneous “snapshot” of a given chemical system but yield only limited information about its chemical evolution. But it is precisely the evolution of the molecular motion that gives rise to chemical processes. Very often, these processes, especially when complex materials are involved, stem from emerging properties which do not simply appear to arise from the combination of those of the molecular constituents.[4] Emerging patterns and complex properties surface because the simpler constituents behave collectively over a given time span in certain ways and not in others. This means that statistics and probability (i.e., entropy) play a crucial role in determining the outcome of a given chemical phenomena.

The problem with dynamic is that the quantum evolution of given chemical system of certain complexity is impossible to calculate. The trick that computational chemists use consists of exploiting the BO approximation and blend the nuclei into ideal classical objects whose quantum features can be neglected. In most applications, these objects are meant to be “atoms”, but, in principle, they can also represent other entities.[5] This approach is broadly known as Molecular Dynamics (MD).[6] MD is indeed based on treating the evolution of the atoms using Newton’s laws of motion where the forces are determined by the electronic energies. In this context, disregarding the quantum nature of the heavy particles produces a set of systematic errors that lie in neglecting zero-point-energies and tunnelling, neglecting coherence, symmetry, and exchange effects, and ignoring the discrete nature of the energies of the molecular motions. Even though some of these issues can actually be accounted for in some way,[7] their inclusion is still missing in most traditional MD techniques.

When it comes to MD, another crucial issue is the origin and the accuracy of the forces that acts on the atoms. In principle one should use the gradient of the electronic energies obtained by solving the electronic SE, but often, in practice, one is limited by size and must resort to a rough parametrization known as the “force field”. The former approach is known as Born–Oppenheimer MD (BOMD) also called ab initio MD (AIMD).[8] The latter consists of the so-called “classical” MD.[9] The forces fields are built using simple analytical expressions that mimic the true quantum electronic potential. Generally, in these force fields, the intramolecular interactions are treated separately from the intermolecular ones. For example, a typical intermolecular atom–atom interaction is modeled using an expression that is the sum of a repulsive term (thus reproducing the exchange interaction at short nuclear distance) and an attractive one (that accounts for the long-range dispersive forces between atoms). The intramolecular forces are treated using elaborate analytical expressions that describe the covalently bound groups of atoms and give shape to the molecular structure.

The choice between BO and classical (force field) MD is dictated by the size of the systems and by the time scale along which the chemical events take place. For systems comprising thousands or more atoms and for time scales larger than a few hundreds of picoseconds, classical MD is the only option. For smaller systems and shorter time scales AIMD becomes feasible. Methods that require a less expensive way to solve the electronic Schrodinger equation, such as semiempirical approaches,[10,11] can also be used to alleviate the computational costs associated with AIMD.

The Materials

As the title says, we focus the following discussion on biocompatible ionic liquids. These liquids are built using relatively simple molecular ions, but due to a surprisingly cooperative dynamics, show a range of complex phenomena that represent the optimal playground to illustrate some of the difficulties and challenges that theoretical and computational simulations will have to confront in the immediate future; especially if we desire to apply them to the increasingly complex problems that emerge from new technologies of societal relevance (new batteries, sustainable materials, new biomaterials, nanomedicine, etc.).

Ionic liquids (ILs) represent a significant opportunity for replacing, at least partly, the solvents currently used in various fields in industry, chemical synthesis, preservation, analysis, etc. To name a few, ILs have been proposed as solvents for electrochemical applications,[12,13] corrosion inhibitors,[14] catalysts,[15] removal agents of polluting gases,[16−21] and dissolving agents for biomasses.[22]

The impressive increase in recent years of the studies concerning ILs in the most varied field of research and technology, is also due to the fact that, for a long time, they have been considered inherently green replacements for traditional chemical solvents (often notoriously harmful to environment and living species). Their extremely low volatility and high chemical stability are indeed huge advantages in case of leaks or laboratory operations. However, recent studies have revisited the presupposed green nature of ILs and discovered that most of them are actually toxic toward organisms and less environmentally benign than previously assumed.[23−26] Despite the large number of studies about the toxicity of traditional ILs, a molecular understanding of their action toward biomolecules is still largely missing.[23]

In the past few years, part of the research effort has steered toward the quest for truly biocompatible ILs. Typically, these liquids are obtained using the cholinium cation, a metabolic harmless substance, as a replacement of the imidazolium of traditional ILs.[27−29] Deprotonated amino acids (see Scheme and Figure ) or simple organic acids (see Figure ) are used as anions giving rise to the subgroups of ILs known respectively as amino acid based ILs (AAILs) and acid–based ILs (AILs). Formally, (A)AILs are obtained using a simple acid/base reaction, hence they pertain to the class of ILs called Protic ILs (PILs) where ionization occurs because of the following equilibriumwhere, to form a PIL, the reaction must be completely shifted toward the products, when all neutral components (the reagents) disappear.

Scheme 1

Isomers of the [Asp]− Anion That May Be in Equilibrium inside the Bulk Phase of an AAILs

Figure 5

Lowest energy, optimized structures of the [Cys]22– (left) and [Glu]22– (right) dianionic dimers as obtained from ab initio calculations in the bulk phase. Reprinted with permission from ref (41). Copyright 2018 American Chemical Society.

Figure 6

Left: Time evolution of the N–H distances of the triethylammonium ions in the bulk phase when coupled to acetate. The initial stage of proton migration from ammonium to acetate is illustrated by the two red trajectories. Right: The final state of the liquid—in red the neutral phase; in blue, the ionic one.

ILs, in general, show a rather complicated pattern of structural features at the nanoscopic scale.[30] This includes segregation phenomena,[31−33] self-assembly,[34,35] and unusual HB networks.[36−42] PILs, in particular, are characterized by the existence of a pervasive network of hydrogen bonds (HBs) that makes them unique in the vast population of possible ILs. The presence of these HBs is the source of their peculiar solvation properties that has spawned a concrete interest from the biomedical community,[43,44] for possible applications in circular economy,[45−47] and in electrochemistry.[12,48,49] This very same network of HB spans a range of nanoscopic phenomena that makes the study of these substances particularly intriguing from the computational point of view, but it also represents a challenge for current state-of-the-art techniques.

Methods: State-of-the-Art

The computer simulations of ILs are performed through a variety of techniques. These are generally chosen depending on the size of the simulated system. Roughly speaking one can divide the continuous size scale (from hundreds of picometers of isolate molecules well into the nanometer range of macromolecules) into three different realms (Figure ): isolated molecules or dimers, clusters composed by a small number of molecules, and bulk. As shown in Figure , as the size of the system (i.e., the number of atoms or entities making up the simulation) increases, the level of detail which we obtain from the simulations tends to decrease. Although there might be significant exceptions, Figure also illustrates the typical appropriate computational method for each size.

Figure 1

Sketch indicates the three system sizes typically explored by computational chemists in the simulations of ILs. A reasonable indication of the appropriate method is also reported.

Typically, when moving from isolated molecules to clusters to bulk, one has to compromise (due to the computational cost) and loose a certain level of details. While ab initio electronic structure calculation such as those based on electronic wave functions or DFT are available at the molecular and cluster level, they rapidly become impractical when simulating the bulk phase that is the realm of methods based on the adoption of parametric forces (the force fields) where the details of the electronic structure and connected phenomena are lost. It is worth pointing out that the implicit hierarchy in Figure is based exclusively on considerations of system size and details, but it does not necessarily entail any assumption about the accuracy of the methods involved. Often a force field-based MD can better match experiments than ab initio depending on how well the force field has been parametrized.

The advantage of an approach based on the study of isolated molecular entities or dimers is that all the details of the electronic and nuclear motions are accessible and that the accuracy of the calculation can be pushed to the boundary of what is possible today. The main drawbacks of such approach consist in the lack of many-body effects and in the impossibility of modeling all those patterns and phenomena that clearly emerge because of collective or cooperative molecular behaviors (segregation, assembly etc.). When instead one works at the bulk scale (disregarding the case of perfect crystals where periodicity translates into a simplified computational model), one loses the electronic details, but there is a chance of obtaining information on collective phenomena.

Between these two regimes, there is space for an approach based on the analysis of small aggregates of molecules or ions. Although this is a level only seldom explored, it represents the ideal approach for looking at how the molecular properties are gradually affected by an increased presence of a surrounding. The advantage of an approach based on a cluster is that highly detailed calculations are still doable together with the possibility of implementing nontrivial dynamic analysis as well. Its drawbacks substantially lie in the inevitable presence of border effects due to the finite size of the clusters.

It is rather evident from the above discussion that to sample the behavior of relatively complex materials such as ILs, considering the “dynamic” factor is essential to understand the origin of their properties and to provide predictive and reliable models on which the planning of new technologies could be implemented. For this reason, I will concentrate mainly on dynamics in the rest of the paper.

The aims of performing MD on ILs is manyfold: first and foremost, it is needed to accumulate a set of configurations that should span the dominant regions of the conformational (phase) space of the system; this, in turn, allows the calculation of thermodynamic functions (energy, enthalpy, free energy, molar volume, etc.) using time or ensemble averages using statistical mechanics. Second, MD can be used to study the structural features of the system, albeit in an averaged fashion (X-ray diffraction patterns, for example). Finally, it provides the frictional properties of the liquid (viscosities, conductivities, diffusion etc.) which are linked to the kinetic of the molecular motion inside fluid itself, with these typically obtained using Green–Kubo relations.[50−52]

MD can be performed in various ways, the most common being ab initio MD, classical or force-field based MD and coarse-grained MD. The latter[53,54] remains outside the scope of the present paper and we defer the reader to the specialized literature. The remaining two variants of MD have both advantages and drawbacks. We summarize them with the scheme of Figure .

Figure 2

This sketch illustrates the main features of the two different ways of implementing MD in the typical setups of ILs theorical studies.

AIMD is based on the simultaneous treatment of electrons and nuclei. The forces acting upon the nuclei are calculated via a solution of the electronic Schrödinger equation (typically with DFT) and differentiating the resulting electronic energy. The chemical topology (i.e., the presence of certain chemical bonds and the coordination of atoms) emerges naturally from the electronic structure and is not fixed. It is ideal to study the dynamic of chemical reactions, or any other bond breaking/forming process and it includes all nuclear many body potentials. It is extremely computationally expensive, and its applications are limited to a system with thousand atoms or less and for processes acting on a scale of picoseconds.

Classical MD, on the other hand, is much cheaper computationally, but the electrons disappear, the nuclei are blended into atoms, and the potential acting between them is parametric and (for the large part) made only by two-body terms. Such a simplified expression of the interatomic potential completely bypasses the need to evaluate forces using explicit electrons and results in a performance gain of several orders of magnitude with respect to ab initio methods. Typical systems for classical MD are well into the nanometric scale and in the microsecond regime.

Introducing (partially) many-body effects in classical MD is possible by using the so-called “polarizable” force fields (vide infra) that include electrostatic effects due to induced dipoles.[55−58] In classical MD, the chemical topology of bonds is fixed and decided a priori. This last condition can be relaxed using special (“reactive”) force fields, that allow a change in atomic coordination, but have only seldom been applied to ILs.[59]

In order to realistically simulate a chemical event using MD one should, in principle, run the simulation long enough to see that event occurring on a statistically significant scale. In the case of ILs the sampling of the phase space for obtaining meaningful thermodynamics averages is a challenge per se. ILs are extremely viscous[60] due to strong electrostatic forces between the molecular components, and collecting data with MD can prove frustrating. An idea of the sluggish motion of the centers of mass of molecular ions can be had by looking at Figure that contains, for the sake of providing an example, data extracted from a simulation done with classical MD on a traditional IL. There I report 1 ns of motions of the centers of mass of 3 selected 1,3-dimethylimidazolium cations (blue) and 3 Tf2N anions (red) against the background (in shaded colors) of the other ions. In 1 ns, each ion does little more than oscillating around its initial position (the side of the cube is 60 Å wide), and it is clear that describing ionic motions needs very long time scales.

Figure 3

Exemplar results from a classical MD simulation of [MMIm][Tf2N] that illustrate the sluggish motion of molecular ions on a timespan of 1 ns. Against a shaded background, I have highlighted six ions (three cations in blue and three anions in red) and shown their centers of mass motion as a series of dots.

Such example should illustrate how properties that are related to the ionic motion (ionic diffusion, ligand exchange, rearrangements of solvation shells etc.) requires extreme simulation times that lie well beyond the present limits of methods based on a detailed calculation of the forces (ab initio MD, Figure , left). In other words, the more accurate the interactions, the less likely is to sample the time scales necessary to converge the desired property. To access long time scales, one has necessarily to resort to a simplified form of the interaction potentials, hence, to force fields-based methods. As I will detail below, the main challenge in performing simulations of ionic liquids stems from this apparent irremediable juxtaposition of accuracy vs time scale.

The problem lies in the fact that the growing use of ILs in complex technological setups that are beyond their simple use as bulk solvents (protein extraction, antimicrobic agents, energy storage devices, catalysis, advanced materials, etc.) already requires, and will do so even more in the future, an increasingly detailed description of their interaction and dynamic that goes well beyond classical MD models. The knowledge of these details has a cost that is linked to the long simulation times and that, at the moment, has to be paid in terms of computational resources. A great effort has been produced in the past and is presently put forward by the computational community in terms of optimizing and improving current algorithms and procedure to achieve the objective of allowing more and more detailed and reliable simulation techniques to become the de facto standard for the description of these systems.

Increasing Details in Classical MD: Polarization and Charge Transfer

It is well-known to the ILs community that, due to ionic nature of the materials, polarization is a crucial phenomenon that should be considered in the simulations. In nonpolarizable MD simulations, the atoms of a molecule interact with the atoms of another one via a potential made by three pieces: a repulsive, short-range part that defines the atomic volume; an attractive part that takes into account the van der Waals cohesive energy; and an electrostatic contribution due to the presence of a partial charge on each atom. In nonpolarizable MD, the partial atomic charges are fixed and depend on the nature of the atom and on the functional group to which it belongs. The atoms interact with each other solely via this potential; hence, traditional MD is entirely based on intermolecular forces that emerge from the sum of two body (atom–atom) interactions. When dealing with systems with partial charges significantly different from zero, the above approach begins to show some limitation. In the real system, induction effects becomes important, and the atomic charges fluctuate depending on the chemical surrounding of the atom. As an example, consider the solvation of a bivalent ion in water: the charge density of the ion induces dipoles on the water molecules that alter their charge distributions and that interact with each other. These interactions between the induced dipoles are many-body effects due to the simultaneous presence of at least three bodies (one inducer and two induced molecules). The potential used in traditional MD must be modified to account for the presence of these effects. This is the purpose of “polarizable” force fields. Using two different algorithmic procedures,[57,58] they introduce a way to allow for charge fluctuations and to account for polarization and induction energy.

In practice, polarization effects are particularly important in fully ionized systems as ILs[55] where these many-body interactions between the induced dipoles tend to reduce the electrostatic cohesive energy, thus explaining why polarizable force field lead to more reliable calculations of dynamical quantities such as fluidity, ionic mobility, and conductivities, all of which are heavily underestimated by fixed charge models.[55,61]

A part of the community is already involved in the setup of polarizable force fields which are presently recognized as an essential addition to traditional MD. A summary of these efforts has been recently reviewed[55] and we only mention the recent availability of the CL&Pol[56,57,62] and Amoeba-IL[58] force fields. The two force fields are based on two entirely different algorithmic approaches, but both aims at becoming a sort of “standard” for polarizable simulations of ILs. It would be auspicious that in the future a wider dissemination of these force fields in the relevant community could provide the necessary validation to assess their potential and reliability.

As we mentioned above, the introduction of polarization is crucial whenever dynamics quantities are involved in the calculations. While structural and static properties can be easily obtained with a high degree of confidence also from nonpolarizable force fields, all those physical properties and phenomena that are related to the reorganization of the ILs molecules, hence friction, mobility, solvation, and transport, require the use of more sophisticated electrostatic models that go beyond the traditional fixed-charge schemes. Unfortunately, introducing polarization has a cost and the performance of MD degrades. An improvement of the algorithms is therefore mandatory to allow polarizable simulations to become the de facto standard method for producing reliable simulations in the ILs field.

Closely connected to polarization there is another peculiar phenomenon that affects the electrostatic of ILs and that is generally referred under the name of charge transfer. This is a quantum effects due to electronic delocalization across molecules. It consists of the fact that, when modeling ILs, the solution of the electronic SE (hence an ab initio approach) clearly indicates that the charge on each molecular ion is not unitary. In practice this means that the overall charge of cations and anions is lower than what is expected from the chemical structure. This effect is linked to the tight packing of the molecular ions caused by the strong electrostatic potential and to the transfer of significant portions of the electronic wave function across them.[63]

One example of this phenomenon is illustrated in Figure , where I report a calculation of the charge distributions of the two molecular ions of triethylammonium mesylate obtained through the approximate resolution of the electronic SE and sampled along 10 ps evolution. The net charge of the two ions (mesylate, negative, and triethylammonium, positive) is never ±1, but averages at ±0.93.

Figure 4

Charge distribution of the molecular ions in triethylammonium mesylate according to a semiempirical treatment of the electronic degree of freedom (DFTB) in a bulk simulation of 10 ps.

The charge transfer phenomenon is another important factor in explaining why the frictional properties of an ILs are overestimated by fixed charge models that assume an integer charge on the molecular ions. This effect is inherently quantum mechanical and cannot be easily taken into account using a model potential. Reducing the overall ionic charges by uniform scaling or introducing a dielectric screening are known and broadly used workarounds,[64] but they are unsatisfactory from various points of view: the scaling factor is peculiar for each liquid and has to be known a priori from ab initio data or from reverse fitting experimental data; in addition, a uniform scaling does not produce a charge distribution such as that in Figure but, rather, a single value of the ionic charges. More elaborate schemes to account for these effects are under study in force fields such as Amoeba+[65,66] and in methods such as those based on the effective fragment potential.[67]

Presently, the best way to accurately take into account the phenomena of charge transfer is to perform AIMD, where the resolution of the SE at each step of the time-evolution ensures the calculation of forces that include polarization, many-body effects and, precisely, charge transfer. It is well-known, however, that AIMD require a huge computational effort that renders it, at the moment, impractical for the study of ILs over a time scale longer than few hundreds of picoseconds. We know that the bottleneck of the AIMD simulations is the resolution of the electronic SE and that any way of reducing the cost of this step opens the way for longer time scales and larger systems. One of the most promising way of implementing a reduction of the computational cost of AIMD is to resort to semiempirical approaches to solve the SE, hence drastically reducing the amount of linear algebra and iterations required at each step of the dynamics.

One of the approach that has been recently utilized in the ILs context is the density functional tight binding method (DFTB)[42,68−70] that has shown good accuracy and performances. In its most recent implementation, the method is based on an expansion of the energy functional up to third order, it includes dispersion energy corrections and naturally produces realistic fluctuating atomic charges. It relies on the a priori knowledge of a set of parameters appearing in the energy expression that are calculated from atomic properties. This method is the one that produced the data reported in Figure .

It seems that semiempirical methods have been generally overlooked by the community as a possible tool to approach an almost first-principle description of ILs. While performance still is an issue for the wide applicability of these methods, there is a wide space for them in the near future especially in light of their ability to incorporate the electronic degrees of freedom in specific situations[71] and to overcome several severe limit of classical MD (see Section ) without being excessively expensive.

Modeling Collective Effects: Like-Charge Aggregation and Other Phenomena

It has been apparent from the beginning that, due to the variety of shapes and functional groups in the molecular ions, ILs can present a wide range of nanoscale structures whose nature is still not entirely understood. From the computational point of view, the emergence of these elusive, transient structures, which nevertheless leave their mark on the bulk liquid properties, is a fascinating, but difficult challenge. The nature of the conundrum is the following. Capturing the nature of these features at the molecular level using atomistic simulations involves operating with sizes large enough to accommodate these (often nanosized) structures for very long time scales.[72] The problem is that to reach such regimes, computational efficiency requires the most simplified potential models (typically nonpolarizable force fields), but the lack of details in the potential might jeopardize the discovery of the driving force of the phenomena.

The presence of such nanostructures in the imidazolium-based generation of ILs, is clearly indicated by the appearance of low-Q peaks in the diffraction patterns.[34,73−75] The aggregation and coalescence of alkyl groups in apolar domains is very likely at the origin of such spectral evidence. Computational studies, mostly based on classical MD in this context, have only provided indirect evidence about the nature of these domains, since the scales at which these aggregations occur are large both in terms of sizes and time. Only recently, more rigorous analysis have shed light on how these segregated domains of alkyl groups have a consequences on bulk properties such as structural relaxation and viscosities.[33,76]

In PILs the network of HBs represents an additional source of disorder and the formation of segregated apolar domains is often suppressed (for peculiar cases see refs (34 and 35)). However, even in the presence of conspicuous HBs, the formation of peculiar structures is not impossible. Our recent works on AAILs (rich in HBs)[77,78] and the work of Ludwig and co-workers[36,37,39,40] have shown that the presence of HBs can promote a different form of aggregation that defies our chemical intuition. This aggregation phenomena is due to like-charge interactions and seems to be peculiar of liquids rich in HBs. It consists of the formation of dimers or oligomers of ions of the same charge that are connected by strong HBs and where the expected Coulomb repulsion is weakened. The repulsive potential turns out to be mitigated by the charge delocalization in certain functional groups (for example carboxylates or sulfonates), by charge transfer which reduce the overall charge of the ions and by the electrostatic screening of the surrounding medium (polarization).

For example, the anionic dimers [Cys]22– and [Glu]22–, shown in Figure , are surprisingly stable when immersed in the bulk and their binding energies are comparable to those of the cation/anion counterpart.[79] This is completely counterintuitive.

To correctly account for such elusive aggregation phenomena the computational model must include all the effects that are responsible for them. A purely fixed charge model would eventually lead to the dissociation of like-charge aggregates under their mutual repulsion. It is evident from the cited literature, that an accurate modeling of polarization and charge transfer are both prerequisites to describe the transient existence of these like-charge aggregated domains and the peculiar role of strong HBs in this context.

Accounting for Variable Chemical Topology: Ionicity and Reactions

PILs, as discussed in Section , arise formally from and acid–base equilibrium reaction. When the equilibrium is not entirely shifted toward the ionic pair, the ensuing liquid is a mixture of a neutral (polar) phase and an ionic one. In this case the system can undergo phase separation to various extent and to evaporation of the polar volatile component.[80,81]

With the help of recent works,[82−84,69] a more complex picture than expected is emerging. Several mixtures of acids and bases that were expected to form fully ionic liquids are instead partially neutral even when the difference of pKa between the two components is large. In addition, the extent of ionization does not seem to be trivially connected to this difference, but rather to the self-solvation abilities of the liquid toward its own ion. This issue is known as the “ionicity problem”. Predicting from molecular properties the extent of equilibrium in the bulk fluid is nontrivial, and further studies and simulations are needed in the near future to fully address the ionicity of PILs.

The coexistence of a neutral and ionic phase makes the computational treatment of these liquids lying at the border of the ionic-to-polar transition extremely difficult. This is not an issue that can be simply solved by focusing on acid/base pairs that give rise to a full ionization. There are, in fact, many instances where PILs are chosen because of their ability to sustain proton conductivity.[85−89] Proton conduction is achieved by using specific mixtures with partially dissociated acids or by precisely exploiting equilibrium . In all these cases the fate of proton transfer reactions between the molecular ions or within a mixture must be properly described by the computational model.

The problem is that traditional MD (even the polarizable variant) are based on a fixed chemical topology where the H atoms (protons) are initially bound to one molecular entity and remain bound to it during the entire simulations. Hence, traditional MD is unable to account for the proton transfer process, proton conductivity, tautomerization reaction, etc. To explore the dynamics of mobile protons the only option is to abandon fixed topology. This can be achieved in two ways: the first is to adopt first-principles MD either ab initio or semiempirical, and the second is to use a force field where the coordination number of each atom is not fixed but changes in response to the chemical environment.

The former approach has been described above and, since chemical topology emerges naturally from the electronic SE, such method is, in principle, able to describe the evolution of any chemical reaction including proton transfer. To provide an example, we show in Figure (left) the evolution of the N–H distances in the liquid triethylammonium acetate calculated using the DFTB method along 350 ps. The system at t = 0 is prepared as a fully ionic system (protons on the ammonium) and all N–H distances are around 1.1 Å. After 125 ps, a first proton hops onto one of the carboxylates neutralizing a pair of molecules (red trajectory). After 300 ps another one does the same thing. Along the dynamics one can also notice how several ammonium protons temporarily break the N–H bond (blue trajectories).

At the end of the proton migration process, the equilibrium state of the fluid should resemble the situation depicted in Figure on the right. The red regions are the volumetric representation of the neutral phase, the blue one the aggregated region where the ionic species reside. It is clear from this simple illustration how chemical reactions and complex electrostatic forces do contribute to shape the unconventional nanoscopic structure of the fluid and of the aggregation/segregation phenomena within.

Most of the biocompatible, new generation ILs are obtained using amino acids anions and the cholinium cation. Choline is a very poor base and equilibrium is entirely shifted to the right. This however does not mean that the problem of proton transfer is irrelevant for these liquids. For example, several amino acids have protic side chains (−SH, −PO3H2 or −COOH) that may lose the proton and lead to isomerization reactions that cause the appearance of zwitterionic-anionic forms of the amino acid. In addition, the spatial vicinity of anionic species due to anionic clustering (Section ) makes the proton exchange between these protic groups and −NH2 very efficient. Hence, like in the case of incomplete proton transfer, the presence of additional (mobile) protons in AAILs could induce intra- and intermolecular isomerization reactions. For example, the singly deprotonated aspartate anion can exist in a set of almost isoenergetic isomers shown in Scheme . These three isomers are presumably in equilibrium through continuous and reversible proton transfer processes and their relative abundances can have a sizable effect on the dynamic and thermodynamic properties of the bulk phase.

Accounting for the presence of these phenomena in a simulation is not easy. The problem is that the typical time scales for establishing proton transfer equilibria can be very long, especially in highly viscous systems as PILs. While semiempirical and ab initio MD, are viable options to describe such equilibria, their typical time scales might lie well outside the reach of the current implementations of these methods. In other words, due the high viscosities of the fluids, simulations long enough to reach a steady state for the isomerization reactions are impractical and massively time-consuming. This means that, at present, the only methods that are appropriate to the study of proton conduction or partial ionization are still far from being able to provide a complete description of these processes simply because they cannot be easily extended to the time scales necessary to reach thermodynamic equilibrium.

When chemical reactions are involved, a more efficient approach is available and consists of using force fields where the atomic coordination numbers are variable and depend on the chemical surrounding. Such possibility is, however, still unexplored in the realm of ILs, and its use is certainly complicated by the fact that the transfer of protons also changes the charge and geometry of the molecular ion. While such force fields exists in other contexts,[90,91] as far as I know, they have never been applied to proton transfer in ILs. The only case that I am aware of is the simulation of CO2 absorption by the AAIL [P4444][Gly].[59] The problems of such an approach are that one must know in advance all possible reactive events to provide a realistic parametrization and that a large number of ab initio calculations are necessary for the characterization of these events.

Perspectives for the Near Future

Computational modeling of complex materials has become a key ingredient of research. The possibility of screening compounds and substance for a given desired property beforehand experimental investigation is an obvious benefit of the use of simulations. Beyond this, the growth in the availability of computer simulations coupled with improved algorithms and with the surge of hardware technical advances have allowed in the past few decades an unprecedented number of researchers to dedicate their work to the fundamental understanding of the phenomena that characterize new materials such as ILs. Despite these advances, though, many challenges still must be faced, some of which I have reviewed in this paper. The most prominent are those connected to the reliability of the models of the intermolecular potentials and to the applications of well-established techniques in situations where additional complications due to chemistry arise. In this context, I have laid out the status of computational research. What is to come is more difficult to foresee, but a glimpse of future development can be, nevertheless, attempted.

Accessing to reliable methods for the modeling of proton transfer dynamics in ILs is an important question: not only because ILs are currently widely employed as prototypical conductive medium in energy devices but also because their interaction with other molecules would necessarily depend on their acidity/basicity properties. I am especially thinking about the interaction of ILs (and their subgroups based on biocompatible molecular ions) with biomolecules and, ultimately, with living organisms. If the use of ILs should become so pervasive as the current research may induce one to foresee, it is evident that a proper modeling of their toxicity, biocompatibility, and biodegradability will be mandatory in the near future. While experimental evidence is accumulating,[23,92,93] simulations and computational approaches about this subject are still scarce.[93−95]

In the context of proton mobility, a substantial advance in the current simulation techniques is necessary and 2-fold: on the one hand, the development of reliable potential models (semiempirical and reactive force fields) is stringent; on the other, the need of accounting for quantum nuclear effects might become an additional challenge.[7] Tunnelling and other phenomena connected to the quantum nature of the light nuclei (protons) can be relevant to understand the biochemical activity of HB-dominated liquids such as AAILs.

For what concerns the outstanding problems of the computational cost, apart from improving algorithms and specific codes, the research could take advantage of relatively new procedures that exploit recent application of the chemoinformatic models to ILs. These statistical models allow one to trace the structure–activity relationship (SAR, QSAR, etc.)[96−100] avoiding altogether the need to perform actual time-consuming MD and use relatively cheap ab initio calculations on isolated constituents.[101]