Wasserscheid P. & Welton T. (2008). *Ionic liquid in synthesis*, 2nd ed., Wiley-CH, ISBN 978-3-

Welton T. (1999). Room-Temperature Ionic Liquids. Solvents for Synthesis and Catalysis,

Wilkes J.S. (2004). Properties of ionic liquid solvent for catalysis*, J. Mol. Catal. A: Chem.*,

Xu J.M., Qian C., Liu B.K., Wu Q. & Lin X.F. (2007). A fast and highly efficient protocol for

Yang L., Xu L.W. & Xia C.G. (2005). Highly efficient KF/Al2O3-catalyzed versatile hetero-Michael

Yoshizawa M. & Ohno H. (2001). Synthesis of molten salt-type polymer brush and effect of

Yu B., Zhou F., Wang C. & W. Liu (2007). A novel gel polymer electrolyte based on poly

Zaderenko P., Soledad Gil M. & Ballesteros P. (1994). Synthesis and regioslective hydrolysis

Zilkova N., Zukal A. & Cejka J. (2006). Synthesis of organized mesoporous alumina

*Tetrahedron Letters*, vol.46, n°19, pp.3279-3282, (May 2005), ISSN 0040-4039 Yang L., Xu L.W., Zhou W., Li L. & Xia C.G. (2006). Highly efficient aza-Michael reactions of

Michael addition of N-heterocycles to α,β-unsaturated compound using a basic ionic liquid [bmIm]OH as catalyst and green solvent. *Tetrahedron*, vol.63, n°4,

addition of nitrogen, oxygen, and sulfur nucleophiles to α,β-ethylenic compounds.

aromatic amines and N-heterocyles catalysed by basic ionic liquid under solventfree conditions. *Tetrahedron Letters*, vol.47, n°44, pp.7723-7726, (October 2006), ISSN

brush structure on the ionic conductivity, *Electrochimica Acta*, vol.46, n°10-11,

ionic liquid 1-ethyl 3-(2-methacryloyloxy ethyl) imidazolium iodide, *Eur. Polym. J*.,

of 2-imidazol-1-ylsuccinic ester, *J. Org. Chem,* vol.59, n°21, pp.6268-6273, (October

templated with ionic liquids. *Microporous Mesoporous Mater.*, vol.95, n°1-3, pp.176–

*Chem. Rev.*, vol.99, n°8, pp.2071-2083 (July 1999), ISSN 0009-2665

vol.214, n°1, pp.11-17, (May 2004), ISSN 1381-1169

pp.986-990, (January 2007), ISSN 0040-4020

pp.1723-1728, (March 2001), ISSN 0013-4686

vol.43, n°6, pp.2699-2707 (June 2007), ISSN 0014-3057

vol.45, n°5, pp. 1577–1582 (March 2004), ISSN 0032-3861

0033-4545

1521-3773

0040-4039

1994), ISSN 0022-3263

179, (October 2006), ISSN 1387-1811

(August 2003), ISSN 1359-7345

527-31239-9, Weinheim.

household microwave oven , *Pure Appl. Chem.,* Vol. 73, No. 8, pp. 1309-1313, ISSN

supported ionic liquid catalysts (silca) in the hydrogenation of α,β-unsaturated aldehydes. *Chem. Eng. Sci.*, vol.62, n°14, pp.3660–3671, (July 2007), ISSN 0009-2509 Washiro S., Yoshizawa M., Nakajima H. and Ohno H. (2004). Highly ion conductive flexible

films composed of network polymers based on polymerizable ionic liquids, *Polymer*

Catalysis, *Ang. Chem. Int. Ed*., vol.39, pp.3772-3789, (November 2000), Online ISSN:

New functionalised ionic liquids from Michael-type reactions—a chance for combinatorial ionic liquid development, *Chem. Comm*., vol. 39, n°16, pp.2038-2039,

## **Theoretical Description of Ionic Liquids**

Enrico Bodo and Valentina Migliorati

*University of Rome "La Sapienza", Chemistry Department Italy*

#### **1. Introduction**

Ionic liquids are universally considered to be materials of the future. Their peculiar properties appeal to the most diverse technological areas such as chemical industry, electrochemistry, optics, environmental chemistry, medicine and nanotechnology. Their extraordinary macroscopic properties originate from the peculiar ionic nanoscopic structure and the possibility of modulating these properties depends ultimately on small modifications of the material at the molecular level.

Theoretical simulations of ionic liquids are becoming increasingly important in predicting the properties of these new materials and in helping to interpret the experimental determinations. From the theoretical point of view, ionic liquids are peculiar materials because they are entirely composed of ions. They can be defined as pure electrolytes. All the technologically interesting properties that are proper to this class of compounds are due to the strong electrostatic interactions that govern the fluid behavior. Though these interactions are quite strong, the size and nature of the cationic partners (more often) are such that the fluid remains liquid even at room temperature. Its behavior, from a microscopic point of view, is fairly similar to the more conventional high temperature molten salts.

The theoretical framework in which it is possible to provide high quality studies of the microscopic structure of the ionic liquid is mainly represented by classical molecular mechanics and, only very recently, by ab-initio molecular dynamics. While the employed theoretical techniques are not very different from those used for conventional fluids, many difficulties arise because of the microscopic nature of ionic liquids. In particular these substances are extremely viscous and simulation times become quickly prohibitive if one wants to describe dynamical properties, even as simple as diffusion coefficients. Recent technological advances such as the introduction of GPU clusters might allow unprecedented possibilities in the simulation of these material opening the route to the simulation of rare events and long time scale phenomena.

The choice and validation of the many force fields which have been established in the last twenty years, some of which are explicitly tailored for a specific class of ionic liquid, represent a crucial step in performing classical trajectory studies. This task has to be performed by taking into account many different sources of experimental data: on the one hand structure studies with X-ray or neutron diffraction provide indications on the short-range structure, on the other hand measurements of bulk properties such as evaporation enthalpies point to the overall energetic behavior of the fluid.

In this chapter we would like to summarize and review the recent developments and possible future opportunities of theoretical simulations of ionic liquids. In particular we will review some of the issues connected to the theoretical simulations of ionic liquids, focusing on our recent work on a few compounds. We shall also review the theoretical methods in general and the peculiarities in their application to ionic liquids.

## **2. Ionic liquids**

Ionic liquids are one of the most promising classes of new materials investigated in the last decade. They do not easily fit the conventional description of molecular fluids therefore promoting a necessary exploration of their physical properties at a microscopic level. Conventionally, ionic liquids are chemicals entirely made by ions (Freemantle, 2009; Gaune-Escard & Seddon, 2010; Rogers & Seddon, 2005a) that show a melting point lower than 100◦C, therefore they are liquid under conventional ambient conditions. Their negligible vapor pressure, high thermal stability, properties tunability upon slight changes in the chemical architecture (including polarity, hydrophobicity, density, solvating activity etc.) have made these materials tailored for a constantly increasing range of applications (Plechkova & Seddon, 2008a; Rogers et al., 2009; Rogers & Seddon, 2005b). Among these applications we find: catalysis, synthesis, sensoristics, medicine, electrochemistry and green chemistry in general (Fischer et al., 1999; Galinski et al., 2006; Garcia et al., 2004; Hough & Rogers, 2007; Hough et al., 2007; Welton & Wasserscheid, 2008). Their versatility has opened the possibility of their employment as functional advanced materials, mediums for materials production and components for preparing highly functional materials. Even if ionic liquids are commonly utilized as solvents, their specific composition where ions solely exist provides amazing functional properties such as dissolution of bio-related materials that never dissolve in conventional solvents. The dissolution of biopolymers as cellulose (Fujita et al., 2009) is one of such examples. Chloride based ionic liquids dissolve cellulose better than other solvents because of the hydrogen bonding between chloride anions with the hydroxyl groups of the polymer. The use of ionic liquids allows a simple, benign system for the processing of cellulose and has potential environmental and cost advantages over current processing methodologies. Ionic liquids are made entirely by ionic couples. The anion is generally inorganic as for example [PF6]−, [NTf2]−, [BF4]−, Br−, Cl−. The cation, instead, is an asymmetric organic cation such as alkylammonium, alkylphosphonium, N,N-dialkylimidazolium, and N-alkylpyridinium. Due to the freedom in designing the organic cation (for example by variation of the side chain length or by varying the substituents on the ring and/or on the chain) and to the possible combination of cation and anion, one can generate a huge number of different ionic liquids. Unlike the about six hundred conventional solvents that are extensively used in industrial and synthetic processes, ionic liquids exists in at least a million pure forms and a trillion ternary mixtures. This liberty allows one to design these materials to optimize a specific task such as a solvent for a reaction. This is the reason why these liquids have been termed "designer solvents" (Fei et al., 2006; Freemantle, 1998; Rogers & Seddon, 2003).

In recent years, ionic liquids have emerged as possible "green solvents", that is environmentally benign substances mainly because they have negligible vapor pressure (Liu et al., 2005). One of the primary driving forces behind research into ionic liquids is the perceived benefit of substituting traditional industrial solvents, most of which are volatile organic compounds (VOCs), with nonvolatile ionic liquids. Replacement of conventional solvents by ionic liquids would prevent the emission of VOCs, a major source of environmental pollution (Polshettiwar & Varma, 2008).

One of the fields in which ionic liquids have more to offer in terms of new opportunities is electrochemistry. While there are many properties that individual members of the ionic liquids family exhibit, the single ubiquitous property that we expect to find in all of them is ionic conductivity (Galinski et al., 2006; Garcia et al., 2004). Unfortunately while the 2 Will-be-set-by-IN-TECH

recent work on a few compounds. We shall also review the theoretical methods in general and

Ionic liquids are one of the most promising classes of new materials investigated in the last decade. They do not easily fit the conventional description of molecular fluids therefore promoting a necessary exploration of their physical properties at a microscopic level. Conventionally, ionic liquids are chemicals entirely made by ions (Freemantle, 2009; Gaune-Escard & Seddon, 2010; Rogers & Seddon, 2005a) that show a melting point lower than 100◦C, therefore they are liquid under conventional ambient conditions. Their negligible vapor pressure, high thermal stability, properties tunability upon slight changes in the chemical architecture (including polarity, hydrophobicity, density, solvating activity etc.) have made these materials tailored for a constantly increasing range of applications (Plechkova & Seddon, 2008a; Rogers et al., 2009; Rogers & Seddon, 2005b). Among these applications we find: catalysis, synthesis, sensoristics, medicine, electrochemistry and green chemistry in general (Fischer et al., 1999; Galinski et al., 2006; Garcia et al., 2004; Hough & Rogers, 2007; Hough et al., 2007; Welton & Wasserscheid, 2008). Their versatility has opened the possibility of their employment as functional advanced materials, mediums for materials production and components for preparing highly functional materials. Even if ionic liquids are commonly utilized as solvents, their specific composition where ions solely exist provides amazing functional properties such as dissolution of bio-related materials that never dissolve in conventional solvents. The dissolution of biopolymers as cellulose (Fujita et al., 2009) is one of such examples. Chloride based ionic liquids dissolve cellulose better than other solvents because of the hydrogen bonding between chloride anions with the hydroxyl groups of the polymer. The use of ionic liquids allows a simple, benign system for the processing of cellulose and has potential environmental and cost advantages over current processing methodologies. Ionic liquids are made entirely by ionic couples. The anion is generally inorganic as for example [PF6]−, [NTf2]−, [BF4]−, Br−, Cl−. The cation, instead, is an asymmetric organic cation such as alkylammonium, alkylphosphonium, N,N-dialkylimidazolium, and N-alkylpyridinium. Due to the freedom in designing the organic cation (for example by variation of the side chain length or by varying the substituents on the ring and/or on the chain) and to the possible combination of cation and anion, one can generate a huge number of different ionic liquids. Unlike the about six hundred conventional solvents that are extensively used in industrial and synthetic processes, ionic liquids exists in at least a million pure forms and a trillion ternary mixtures. This liberty allows one to design these materials to optimize a specific task such as a solvent for a reaction. This is the reason why these liquids have been termed "designer solvents" (Fei et al., 2006; Freemantle, 1998; Rogers & Seddon, 2003). In recent years, ionic liquids have emerged as possible "green solvents", that is environmentally benign substances mainly because they have negligible vapor pressure (Liu et al., 2005). One of the primary driving forces behind research into ionic liquids is the perceived benefit of substituting traditional industrial solvents, most of which are volatile organic compounds (VOCs), with nonvolatile ionic liquids. Replacement of conventional solvents by ionic liquids would prevent the emission of VOCs, a major source of

the peculiarities in their application to ionic liquids.

environmental pollution (Polshettiwar & Varma, 2008).

One of the fields in which ionic liquids have more to offer in terms of new opportunities is electrochemistry. While there are many properties that individual members of the ionic liquids family exhibit, the single ubiquitous property that we expect to find in all of them is ionic conductivity (Galinski et al., 2006; Garcia et al., 2004). Unfortunately while the

**2. Ionic liquids**

ionic liquids are intrinsic ion conductors, their ion conductivity often falls short of that of solvent-based electrolytes because of the high viscosities (Hu & Margulis, 2007). Of course viscosity is the result of the electrostatic interactions that are intrinsic to the nature of these substances. Untangling this "catch 22" situation has been the leading theme of research in electrochemistry (Wilkes, 2002). Many of the structural features that make this substances attractive for electrochemistry are also at the origin of the difficulties in using them in the electrochemical processes (Forsyth et al., 2004; MacFarlane et al., 2007; Xu & Angell, 2003). One potential application of ionic liquids in the electrochemical field concerns their use in supercapacitors also known as electrochemical double-layer capacitors (EDLC). Another intriguing application of ionic liquids concerns the construction of dye-sensitized solar cells (DSSCs) (Wang et al., 2002; Wolfbauer et al., 2001). These are novel photovoltaic devices, also known as Gratzel cells, consisting of a layer of a nanoporous semiconductor, often TiO2, impregnated with a light sensitive dye and confined between two metallic electrodes. The microporous structure of the semiconductor is embedded in a liquid electrolyte that contains a redox couple, for instance I−/*I*3 −. Different organic solvents have been tested in the manufacture of DSSCs but, as a result of evaporation and leaking, most of them deteriorate the performance of the cell after some time of operation. Ionic liquids have been proposed as an alternative in the design of DSSCs because of their high viscosity and low vapor pressure that minimize solvent leaking and losses.

Molecular simulations are playing an increasing role in developing an understanding of the physical chemistry of ionic liquids. Unlike conventional organic liquids that have been studied for decades, the ionic liquid research field is still young and much is still unknown about their properties and behavior, and, in this respect, simulations are on an equal footing with experimental investigations (Lipkowitz et al., 2009). This means that molecular simulation will be used both for property prediction and to provide qualitative insight into the nature of these substances. Property predictions is of paramount importance under conditions where experiments are difficult to be conducted and thus cannot determine the structural and dynamic properties of the system under investigation.

#### **3. A brief history**

The history of ionic liquids is a recent one even if one can trace back the origins of the materials to the early 19th century, by examining the molten salt research (Wilkes, 2002). An ionic liquid can be loosely defined, indeed, as a molten salt with a melting temperature below the boiling point of water. In general an ionic liquid differs from a traditional molten salt because is made by an organic cation and an inorganic anionic partner. Although ionic liquids can be handled as ordinary solvents, they present many specific features that are not seen in ordinary organic solvents or even in high temperature molten salts. These features are mainly due to the strong ion-ion interactions that are peculiar to these materials. One of the most important consequences of the network of strong ionic interactions is the negligible vapor pressure and the huge liquidus range i.e. the span in temperatures between the freezing and boiling point. The consequence for green chemistry is that the ionic liquids represent the ultimate non-volatile solvents (Wilkes, 2002).

As stated before, the roots of nowadays ionic liquids research is firmly planted in the traditional high temperature molten salts development and use. Molten salts have long been established as solvents for particular applications, mainly electrochemical: their main advantages are that they present a high electrical conductivity due to ion mobilities, they have a wide liquidus range due to their thermal stability and they are relatively cheap because they are derived from naturally occurring minerals. Their main drawback is that they operate as liquids at very high temperatures so that their use as a reaction media can be considered unknown to most chemists. Ionic liquids present the advantages of molten salts but avoid the problem of the high temperatures.

The first report for the existence of room temperature ionic liquids dates back to 1914 (Walden, 1914a; Walden., 1914b) when Walden investigated ethylammonium nitrate. After this work, many years passed without a substantial understanding of the considerable potential of this class of materials and only in the 1970's the first generation of ionic liquids (composed of water-unstable organic chloro-aluminates) was studied. As soon as air and water stable ionic liquids have been synthesized in the early nineties (Wilkes & Zaworotko, 1992), ionic liquids have attracted the interest of a considerable part of the material science, electrochemical and organic chemistry community. At the present moment, more than 8000 papers have been published in the last decade and review papers are appearing every two or three days (Plechkova & Seddon, 2008b).

From the theoretical point of view the research on ionic liquids lays its foundations in the simulations of molten salts of the late '60. Although between 1990 and 1997 there have been various theoretical studies on ammonium salt derived ionic liquids, the first study of a dimethy-limidazolium compound traces back to 2001 when Hanke et al. (Hanke et al., 2001) studied the [C1mim][Cl] and [C1mim][PF6] at high temperature (mostly because of computational limitations). Another study followed with the dissolution in the same compounds of water, methanol and propane (Lynden-Bell et al., 2002). Many other theoretical investigations of ionic liquids then followed and helped elucidating some of the key properties and features of this class of compounds. These studies have also highlighted the difficulties of obtaining reliable results on ionic liquids from molecular simulations and the challenges of the next generation of computational simulations. One major problem is that the ionic liquid dynamics is slow, much slower than in conventional liquids and this makes computational method more difficult to be adopted unless very high temperatures or very long simulation times are used. Another common feature of ionic liquids, that can be studied through the use of simulations and that was recognized in these systems from the beginning (Morrow & Maginn, 2002), is the spatial etherogenity due to a long range alternating ordering effect between cations and anions which persists up to 20 Å around each ion center. As an example we report here a comparison we have made in ref. (Bodo et al., 2010) between the measured total distribution functions and those predicted by means of molecular dynamics simulations for the series of compounds [C*n*mim][Tf2N] with *n* =1, 4, 6 and 8. Long range ordering effects can easily be seen from Figure 1 up to 20-25 Å. It has to be pointed out that the theoretical predictions have been obtained by using two different force field: the one from Canongia-Lopes (Lopes & Pádua, 2004) and another one from Köddermann (Ködderman et al., 2007). Differences are very small between the two simulations and both force fields are able to reliably reproduce the structural features of the liquid. These study showed that relatively simple force fields as those adopted in our work that are based on the OPLS one (Jorgensen et al., 1996) are reasonably accurate in describing the structural features of ionic liquids. It turns out to be much more difficult to reliably predict dynamical features as we will show below.

#### **4. Present challenges**

A bulk ionic liquid can be considered as a pure electrolyte (Rogers & Seddon, 2005a). It is made entirely of molecular ions and besides its actual composition or the actual nature of the molecular ions that it is made of, many of the appealing properties which have been exploited in the technological applications are due to this feature.

4 Will-be-set-by-IN-TECH

as liquids at very high temperatures so that their use as a reaction media can be considered unknown to most chemists. Ionic liquids present the advantages of molten salts but avoid the

The first report for the existence of room temperature ionic liquids dates back to 1914 (Walden, 1914a; Walden., 1914b) when Walden investigated ethylammonium nitrate. After this work, many years passed without a substantial understanding of the considerable potential of this class of materials and only in the 1970's the first generation of ionic liquids (composed of water-unstable organic chloro-aluminates) was studied. As soon as air and water stable ionic liquids have been synthesized in the early nineties (Wilkes & Zaworotko, 1992), ionic liquids have attracted the interest of a considerable part of the material science, electrochemical and organic chemistry community. At the present moment, more than 8000 papers have been published in the last decade and review papers are appearing every two or three days

From the theoretical point of view the research on ionic liquids lays its foundations in the simulations of molten salts of the late '60. Although between 1990 and 1997 there have been various theoretical studies on ammonium salt derived ionic liquids, the first study of a dimethy-limidazolium compound traces back to 2001 when Hanke et al. (Hanke et al., 2001) studied the [C1mim][Cl] and [C1mim][PF6] at high temperature (mostly because of computational limitations). Another study followed with the dissolution in the same compounds of water, methanol and propane (Lynden-Bell et al., 2002). Many other theoretical investigations of ionic liquids then followed and helped elucidating some of the key properties and features of this class of compounds. These studies have also highlighted the difficulties of obtaining reliable results on ionic liquids from molecular simulations and the challenges of the next generation of computational simulations. One major problem is that the ionic liquid dynamics is slow, much slower than in conventional liquids and this makes computational method more difficult to be adopted unless very high temperatures or very long simulation times are used. Another common feature of ionic liquids, that can be studied through the use of simulations and that was recognized in these systems from the beginning (Morrow & Maginn, 2002), is the spatial etherogenity due to a long range alternating ordering effect between cations and anions which persists up to 20 Å around each ion center. As an example we report here a comparison we have made in ref. (Bodo et al., 2010) between the measured total distribution functions and those predicted by means of molecular dynamics simulations for the series of compounds [C*n*mim][Tf2N] with *n* =1, 4, 6 and 8. Long range ordering effects can easily be seen from Figure 1 up to 20-25 Å. It has to be pointed out that the theoretical predictions have been obtained by using two different force field: the one from Canongia-Lopes (Lopes & Pádua, 2004) and another one from Köddermann (Ködderman et al., 2007). Differences are very small between the two simulations and both force fields are able to reliably reproduce the structural features of the liquid. These study showed that relatively simple force fields as those adopted in our work that are based on the OPLS one (Jorgensen et al., 1996) are reasonably accurate in describing the structural features of ionic liquids. It turns out to be much more difficult to reliably predict dynamical features as we will

A bulk ionic liquid can be considered as a pure electrolyte (Rogers & Seddon, 2005a). It is made entirely of molecular ions and besides its actual composition or the actual nature of the molecular ions that it is made of, many of the appealing properties which have been exploited

problem of the high temperatures.

(Plechkova & Seddon, 2008b).

show below.

**4. Present challenges**

in the technological applications are due to this feature.

Fig. 1. Experimental and theoretical (MD) *G*(*r*) for [C*n*mim][NTf2] with *n*=1 (upper-left corner), *n*=4 (upper-right), *n*=6 (lower-left), *n*=8 (lower right).

The structure of an ionic liquid at the molecular level is therefore characterized by very strong electrostatic interactions which give rise to many of its macroscopic properties: viscosity, negligible vapor pressure, poor conductivity1. The relatively low melting points, on the other hand, are due to the complexity of the molecular components and to their asymmetry which defy the natural tendency of ionic compounds to form crystals.

Another aspect that is important to keep in mind is that, unless the halogen anions are used, the charge is not localized on an atom but is distributed all over the molecular structure, therefore making even more difficult to describe the ionic liquid as made by point charges.

From the theoretical point of view, the description of an ionic liquid in its bulk liquid state still represents a formidable challenge under many point of view. Molecular dynamics simulations are the main (if not the only possible) choice for the molecular, atomistic description of an ionic liquid in the bulk or nanoconfined state. As a first point, consider that if one wants to perform a statistically meaningful simulation of a liquid, one has to provide at least few ns of dynamics using a cell large enough to capture the long range structural correlations of the system under investigation. This means that cells of the order of 50-70 Å have often to be used which generally contain few hundreds of ionic couples. This substantially rules out the possibility of

1 Actually all the ionic liquids are intrinsically conductors, but since the ion drift or diffusion is very low it turns out that only in a limited number of mixture, it has been possible to improve the intrinsic ionic conductivity(Forsyth et al., 2004)

using ab-initio molecular dynamics although it may represent a promising route for the near future or for the simplest (in terms of molecular structure) inhabitants of the vast zoology of ionic liquids. Despite the fact that ab-initio molecular dynamics on such systems is slowly getting within the grasp of modern supercomputer, it still can not represent a viable choice for general applications. This limitation is due to the high computational cost of methods such as Car-Parrinello or even Born-Oppenheimer techniques. The reason is simple: most of the ionic liquids present a very high viscosity, due to the strong electrostatic interactions, and therefore an ab-initio molecular dynamics simulation requires both equilibration and production simulation times that are unaffordable in terms of computational costs. The application of ab-initio techniques, i.e. quantum mechanical methods (QM) is limited, in most cases, to the study of small portion of materials such as clusters or even gas phase isolated molecules. Other issues may arise if one wants to calculate dynamical properties, such as the self diffusion coefficients either from the velocity correlation functions or mean square displacements. In this case the simulation time should be in the hundreds of nanosecond scale and this can be achieved only for small simulation cells or for the simplest compounds. Molecular dynamics simulations can exploit the fact that they bypass the solution of the electronic quantum problem and use instead a very crude representation of the force acting on nuclei generally obtained by means of mixed QM and empirical data. The choice of the classical force field is of crucial importance for the reliability of the theoretical results. Even if the force fields have been validated against various experimental quantities (viscosity, densities, evaporation enthalpy, diffraction data), their transferability to other systems, not

directly involved in the validation stage, can often represent a problem. Given the huge number of molecular variations that are possible with ionic liquids, the use of a given force field has to be experimentally verified anyway, even if the same force field has been successfully used elsewhere. Another challenge that awaits the molecular dynamics simulations is the lack of validated force fields for many of the compounds of interest for experimental applications. For example, while imidazolium based compounds have been extensively investigated and the series of force fields provided by Canongia-Lopes and coworkers (Canongia Lopes & Pádua, 2006b) have been used successfully in predicting the physical properties of the relevant compounds, much less has been done on other compounds such as ammonium based salts although they represent the prototypical ionic liquid. The task of developing and validating a force field is tedious, but it is essential for the theoretical community because, especially if the validation is performed by means of different experimental sources, it provides the basis on which one is able to obtain reliable theoretical data.

It is our strong belief that in the next years theory will not only serve as a post-experiment interpretation tool or a simple "black box" able to predict experimental quantities (even though these already represent highly complex duties), but it will provide a sophisticated tool able to obtain complementary information with respect to experimental quantities and to see many things that are precluded to experiments. In order to do so theory must rely on a well established approximation to the true forces governing the nuclear motion. In this sense, the validation of force fields is crucial.

## **5. The force field**

The problem of choosing, creating or improving a force field for molecular dynamics calculations is obviously non trivial and requires a careful validation of the field with the existing experimental observables. There are even many different paradigms which one may choose from in order to build a force field: for example the "atomistic" level and the 6 Will-be-set-by-IN-TECH

using ab-initio molecular dynamics although it may represent a promising route for the near future or for the simplest (in terms of molecular structure) inhabitants of the vast zoology of ionic liquids. Despite the fact that ab-initio molecular dynamics on such systems is slowly getting within the grasp of modern supercomputer, it still can not represent a viable choice for general applications. This limitation is due to the high computational cost of methods such as Car-Parrinello or even Born-Oppenheimer techniques. The reason is simple: most of the ionic liquids present a very high viscosity, due to the strong electrostatic interactions, and therefore an ab-initio molecular dynamics simulation requires both equilibration and production simulation times that are unaffordable in terms of computational costs. The application of ab-initio techniques, i.e. quantum mechanical methods (QM) is limited, in most cases, to the study of small portion of materials such as clusters or even gas phase isolated molecules. Other issues may arise if one wants to calculate dynamical properties, such as the self diffusion coefficients either from the velocity correlation functions or mean square displacements. In this case the simulation time should be in the hundreds of nanosecond scale and this can be achieved only for small simulation cells or for the simplest compounds. Molecular dynamics simulations can exploit the fact that they bypass the solution of the electronic quantum problem and use instead a very crude representation of the force acting on nuclei generally obtained by means of mixed QM and empirical data. The choice of the classical force field is of crucial importance for the reliability of the theoretical results. Even if the force fields have been validated against various experimental quantities (viscosity, densities, evaporation enthalpy, diffraction data), their transferability to other systems, not directly involved in the validation stage, can often represent a problem. Given the huge number of molecular variations that are possible with ionic liquids, the use of a given force field has to be experimentally verified anyway, even if the same force field has been successfully used elsewhere. Another challenge that awaits the molecular dynamics simulations is the lack of validated force fields for many of the compounds of interest for experimental applications. For example, while imidazolium based compounds have been extensively investigated and the series of force fields provided by Canongia-Lopes and coworkers (Canongia Lopes & Pádua, 2006b) have been used successfully in predicting the physical properties of the relevant compounds, much less has been done on other compounds such as ammonium based salts although they represent the prototypical ionic liquid. The task of developing and validating a force field is tedious, but it is essential for the theoretical community because, especially if the validation is performed by means of different experimental sources, it provides the basis on which one is able to obtain reliable theoretical

It is our strong belief that in the next years theory will not only serve as a post-experiment interpretation tool or a simple "black box" able to predict experimental quantities (even though these already represent highly complex duties), but it will provide a sophisticated tool able to obtain complementary information with respect to experimental quantities and to see many things that are precluded to experiments. In order to do so theory must rely on a well established approximation to the true forces governing the nuclear motion. In this sense,

The problem of choosing, creating or improving a force field for molecular dynamics calculations is obviously non trivial and requires a careful validation of the field with the existing experimental observables. There are even many different paradigms which one may choose from in order to build a force field: for example the "atomistic" level and the

data.

the validation of force fields is crucial.

**5. The force field**

electrostatic model. The former can be assessed in the following terms: in principle the most accurate simulations can be obtained from an all-atom calculation in which, as the name says, all of the atoms in the system are explicitly treated. Moreover, all of the bonds among atoms in the molecules can be allowed to vibrate. However, it is well known that some of the vibrational motions of the involved bonds cannot be correctly described by classical mechanics at least in the room temperature range. The frequencies associated to vibrational motions are often too high to allow for the classical equipartition theorem to work correctly and quantization effects should be considered or included. This is especially valid for the stretching of bonds involving the hydrogen atom. As a first step, it is therefore a common practice to freeze the C-H bonds (or the O-H and N-H ones). For this reason it may become computationally convenient to fuse, for example, the three/four Lennard-Jones entities of a CH2/CH3 group into a new L-J entity that behaves as a rigid unit. This leads to the construction of United Atoms (UA) force fields that reduce the computational complexity by adopting the above scheme. This approximation can be moved one step forward into considering entire subunits (more or less rigid) of the molecular system as unique entity: the resulting coarse-graining models eliminate out those degrees of freedom that enter into macroscopic properties only through their cooperative effect and substitute them by effective degrees of freedom. Building a coarse grained model implies a rather complex fitting procedure which is based on the knowledge of some previous molecular dynamics run. In a first step one has to decide how to map an all-atom system into a reduced complexity system which has a different structure. One then has to build a reliable force field acting between the newly defined subunits. Needless to say, the second step is the most complex one (Wang & Voth, 2005).

Coarse grained models have been successfully employed in the description of ionic liquids (Wang & Voth, 2005; 2006). It has been shown that for sufficiently long side chains in imidazolium based liquids, the apolar alkali chains tend to aggregate. This important result is in agreement with many experimental observations.

From the point of view of electrostatic interactions (that are obviously of crucial importance for the description of ionic liquids), there are two main choices: fixed partial charges and polarizable models.

Among the many possible variants of force fields the OPLS has been the most popular in this field and has been suitably taylored for the description of a large number of classes of ionic liquids by A. Padua and coworkers (Canongia Lopes & Pádua, 2006a; Canongia Lopes et al., 2005; Canongia Lopes & Pádua, 2006b; Lopes et al., 2004; Lopes & Pádua, 2004)

The non polarizable OPLS force fields can be written by using the following functional form (Jorgensen et al., 1996).

$$V = V\_{\text{bonded}} + V\_{\text{non-bonded}} = \sum\_{ij} k\_{ij} (r\_{ij} - r\_{ij}^0)^2 + \sum\_{ij\bar{k}} a\_{i\bar{j}\bar{k}} (\theta\_{i\bar{j}\bar{k}} - \theta\_{i\bar{j}\bar{k}}^0)^2 + \dots$$

$$\sum\_{ij\bar{k}} \sum\_{n=1}^3 \frac{1}{2} V\_n \cos(n\phi\_{i\bar{j}l\bar{l}}) + \sum\_{ij} \left[ \left(\frac{\sigma\_{ij}}{r\_{ij}}\right)^{12} - \left(\frac{\sigma\_{ij}}{r\_{i\bar{j}}}\right)^6 \right] + \frac{q\_{\bar{l}} q\_{\bar{j}}}{4\pi \epsilon\_0 r\_{i\bar{j}}} \tag{1}$$

where the first term is the harmonic vibrational energy of the bonded pairs, the second term is the bending energy, the third term is due to 1-4 torsional interactions and the last two terms are the Lennard-Jones short range term and the electrostatic interaction energy. The latters are normally referred to as non-bonded interactions. Other force fields use different forms of each specific term. For example the MM3 one, has a more complex expression for the bonded terms that includes higher order terms in addition to the quadratic one.

The electrostatic term is very simple in the OPLS form but may be more complex in other force fields. This term is actually crucial in the simulation especially if one wants to calculate energetic quantities such as vaporization enthalpies or diffusion coefficients (see below). If we add a polarization term, this has usually the following form:

$$V^{\rm pol} = -\frac{1}{2}\mu\_i \cdot \mathbf{E}\_1^0 \tag{2}$$

where μ*i* = *αi*E*i* is the induced dipole at force center *i*, with polarizability *αi* (isotropic) and E0 is the force field due only to the partial atomic charges. Atomic polarizabilities are found by fitting the gas-phase polarizability as determined by quantum mechanics calculations. In general it is found that polarizable force fields, when accurately parametrized, (Borodin, 2009) lead to a realistic description of many important properties of ionic liquids such as densities and transport properties.

In any case a force field has to be "validated" against experimental data. The best way to obtain the parameters for the bonded terms such as bond and bending is that of producing fairly high quality ab-initio calculations in the gas phase for the isolated cation and anion and to obtain equilibrium distances and force constants from an optimization followed by a normal mode analysis. The dihedral angle potential is normally obtained by calculating the torsional profile and fitting the energy expression to the raw points. Electrostatic is generally obtained by assigning a partial charge to each atom (see below). Lennard-Jones parameters are much more difficult to obtain from a quantum calculation and generally the parameters adopted in widely available force fields (AMBER, OPLS) are collected and used.

#### **6. A combined approach**

The theoretical simulation of ionic liquids requires a combined approach that involves the use of more than one techniques. It is customary to perform a series of preliminary ab-initio calculations before moving to molecular dynamics simulations. The ab-initio simulations are dictated by the necessity of obtaining the "true" charge distribution in order to check whether it matches the partial atomic charges of the force field and, if necessary, to provide new partial atomic charges.

The first step in the study of a molecular ionic liquid is the modeling of the gas-phase structure of its ionic constituents. This approach is often necessary because the partial atomic charges are not known for many particular structures. Two approaches can be followed: the first one which is more apt to allow for the construction of a general force field is to optimize the isolated cationic or anionic moieties and to calculate the partial atomic charges. The second one is to consider the ionic couple in the gas phase and optimize the two partners altogether. In the former case the net charge on each of the partner will be entire while in the latter case the net charge of the partner can be fractional.

#### **6.1 Charges determination**

The simplest way to calculate partial atomic charges is that of fitting the electric potential generated by the electronic and nuclear charge distribution by point charges placed on the various atomic centers. The charge distribution can be easily obtained by ab-initio calculations using well tested DFT or MP2 methods. The more traditional way of producing atomic charges based on wavefunction atomic orbital analysis (Mulliken or Lodwin population) should be avoided because they are not reliable especially when an augmented basis has 8 Will-be-set-by-IN-TECH

each specific term. For example the MM3 one, has a more complex expression for the bonded

The electrostatic term is very simple in the OPLS form but may be more complex in other force fields. This term is actually crucial in the simulation especially if one wants to calculate energetic quantities such as vaporization enthalpies or diffusion coefficients (see below). If we

> 2 μ*i* · E0

where μ*i* = *αi*E*i* is the induced dipole at force center *i*, with polarizability *αi* (isotropic) and E0 is the force field due only to the partial atomic charges. Atomic polarizabilities are found by fitting the gas-phase polarizability as determined by quantum mechanics calculations. In general it is found that polarizable force fields, when accurately parametrized, (Borodin, 2009) lead to a realistic description of many important properties of ionic liquids such as densities

In any case a force field has to be "validated" against experimental data. The best way to obtain the parameters for the bonded terms such as bond and bending is that of producing fairly high quality ab-initio calculations in the gas phase for the isolated cation and anion and to obtain equilibrium distances and force constants from an optimization followed by a normal mode analysis. The dihedral angle potential is normally obtained by calculating the torsional profile and fitting the energy expression to the raw points. Electrostatic is generally obtained by assigning a partial charge to each atom (see below). Lennard-Jones parameters are much more difficult to obtain from a quantum calculation and generally the parameters

The theoretical simulation of ionic liquids requires a combined approach that involves the use of more than one techniques. It is customary to perform a series of preliminary ab-initio calculations before moving to molecular dynamics simulations. The ab-initio simulations are dictated by the necessity of obtaining the "true" charge distribution in order to check whether it matches the partial atomic charges of the force field and, if necessary, to provide new partial

The first step in the study of a molecular ionic liquid is the modeling of the gas-phase structure of its ionic constituents. This approach is often necessary because the partial atomic charges are not known for many particular structures. Two approaches can be followed: the first one which is more apt to allow for the construction of a general force field is to optimize the isolated cationic or anionic moieties and to calculate the partial atomic charges. The second one is to consider the ionic couple in the gas phase and optimize the two partners altogether. In the former case the net charge on each of the partner will be entire while in the latter case

The simplest way to calculate partial atomic charges is that of fitting the electric potential generated by the electronic and nuclear charge distribution by point charges placed on the various atomic centers. The charge distribution can be easily obtained by ab-initio calculations using well tested DFT or MP2 methods. The more traditional way of producing atomic charges based on wavefunction atomic orbital analysis (Mulliken or Lodwin population) should be avoided because they are not reliable especially when an augmented basis has

adopted in widely available force fields (AMBER, OPLS) are collected and used.

1 (2)

*V*pol = −1

terms that includes higher order terms in addition to the quadratic one.

add a polarization term, this has usually the following form:

and transport properties.

**6. A combined approach**

**6.1 Charges determination**

the net charge of the partner can be fractional.

atomic charges.

been used (which is the normal case for the anionic partner). There are various fitting recipes that can be found implemented in the most common ab-initio packages, but the CHELPG procedure (Frisch et al., 2009) has become a de facto standard for force field calculation when coupled to B3LYP calculations.

A charge distribution and the resulting atomic charges are obviously dependent upon the particular geometry of the molecular system. This is a subtle point that may not be apparent when dealing with rigid ionic partners such as 1,3-dimethyl-imidazolium where the conformational mobility is limited to the CH3 rotation. There are however examples where the conformational mobility is very large and it turns out not to be obvious which gas-phase geometry is representative of the liquid phase. We will make more comments to assess this specific point with an example in the final section of this chapter. It might be not obvious, in a system with a high conformational mobility, to choose one particular geometry in order to calculate charges. A polarizable force field partially solves and accounts for this problem, but it still relies upon a certain specific geometry to provide some important quantities. The example reviewed below represents an extreme case of conformational mobility. The geometry of the isolated cationic partner is completely different from that it assumes when in contact with the anion. The resulting charge distribution is even more diverse.

Another puzzling question which has to do with the determination of partial atomic charges is whether the ionic couple has to be considered as a real entity that survives in the ionic liquid or not. Having a net entire charge on the two ionic partners is obviously a great advantage if the main objective is that of creating a force field as transferable as possible to similar compounds. It has, however, been shown that fractional charges obtained by an optimization of the gas-phase ionic pair, might improve the reliability of the following molecular dynamics runs especially in terms of diffusion coefficients. It is interesting to note that less charged partners lead to less viscous fluids and to higher diffusion coefficients. This fact has a very simple explanation: due to the relatively high density of the liquid phase of these compounds, the ionic partner are considerably in close contact during the liquid dynamics. It is also easy to predict that in an ionic pair there will always be a certain amount of charge transfer (especially in protic ionic liquids). Therefore it might be concluded that a simulation in which the two ionic partners do not have a net entire charge might be more realistic. See (Borodin & Smith, 2006) for an example where a non integer charge on the cation and the anion were matched to an ab-initio simulation.

#### **6.2 The role of dispersion energy**

Given the strong electrostatic nature of the interaction energy it is clear that the contribution of dispersion forces is of secondary importance and limited to short range interactions. There is however a secondary, but still relevant possible contribution. The presence of *π* electron clouds on aromatic rings on the cation, and of long alkyl chains can make the importance of London interactions more important than one would initially consider. For example in ref. (Zahn & Kirchner, 2008; Zahn et al., 2008), the authors assess the importance of London forces on the structure of selected gas phase structures of typical ionic liquids pairs. They conclude that dispersion forces can play a role in determining the structure of ionic liquids. They also make a comparison of the performance of various GGA, meta-GGA and dispersion corrected density functionals concluding that in order to obtain reliable dissociation energies one has to use dispersion corrected functionals.

#### **6.3 Molecular dynamics**

A recent edition of Account of Chemical Research (see (Maginn, 2007)) contains several reviews focussing on the theoretical problem of modeling ionic liquids. Maginn's chapter in Reviews in Computational Chemistry is also a good starting point for a general discussion of molecular modeling of ionic liquids (Maginn, 2009).

Once a force field has been chosen, a sampling of the liquid phase space has to be performed. This can be achieved both using Monte Carlo and Molecular Dynamics. In the former a series of configurations of the system are generated along a Markov Chain and a transition probability from an element to the other has to be defined in order to allow the best exploration of the phase space. Statistical averages can then be run over each of the Markov Chain elements. In a molecular dynamics trajectory calculation, one solves the Newton equations using a finite step approximation (usually employing highly efficient simplectic integrators such as Verlet, velocity Verlet or Beeman). A snapshot of the system geometry, velocity and forces is then saved at a given time-step (much larger than the integrator time-step) creating what is generally called the "trajectory". The sequence of snapshots (trajectory) is then used to calculate statistical averages.

Here follows a series of examples for various properties that can be computed from the analysis of the trajectory.

In particular, one can calculate the static structure factor either of neutrons or of X-ray (Keen, 2000) from the trajectory of the liquid phase simulation and directly compare it with the experimental one. For X-ray scattering data the radial distribution function between particle type *i* and *j* can be used, which are defined as:

$$g\_{ij}(r) = \frac{n\_{ij}(r)}{4\pi r^2 dr \rho\_j} \tag{3}$$

where *ρj* is the number density of particles *j*, and *nij*(*r*) is the number of particles of type *j* between distance *r* and *r* + *dr* from a particle of type *i*. Given *cj* the molar fraction of particles *j*, *ρj* = *cj* · *ρ*0 where *ρ*0 is the total number density of the sample. The total scattering for X-ray can be then written as:

$$F(Q) = \rho\_0 \sum\_{i=1}^{n} c\_i c\_j f\_{ij}(Q) \int\_0^{\infty} 4\pi r^2 [g\_{ij}(r) - 1] \frac{\sin(Qr)}{Qr} \tag{4}$$

10 Will-be-set-by-IN-TECH

A recent edition of Account of Chemical Research (see (Maginn, 2007)) contains several reviews focussing on the theoretical problem of modeling ionic liquids. Maginn's chapter in Reviews in Computational Chemistry is also a good starting point for a general discussion

Once a force field has been chosen, a sampling of the liquid phase space has to be performed. This can be achieved both using Monte Carlo and Molecular Dynamics. In the former a series of configurations of the system are generated along a Markov Chain and a transition probability from an element to the other has to be defined in order to allow the best exploration of the phase space. Statistical averages can then be run over each of the Markov Chain elements. In a molecular dynamics trajectory calculation, one solves the Newton equations using a finite step approximation (usually employing highly efficient simplectic integrators such as Verlet, velocity Verlet or Beeman). A snapshot of the system geometry, velocity and forces is then saved at a given time-step (much larger than the integrator time-step) creating what is generally called the "trajectory". The sequence of snapshots (trajectory) is then used to

Here follows a series of examples for various properties that can be computed from the

• Density: the easiest to calculate quantity and the more commonly experimentally available one. It is always a good starting point for validating a force field even if it is well known that an incorrect force field can still give rise to a correct density. So in general it is possible to accurately reproduce density while other properties are only poorly reproduced. In general OPLS based force fields underestimate slightly the experimental densities. • Liquid structure: a very strict test of the quality of the force field used in a simulation is represented by the comparison between the structural results obtained from the trajectory and X-ray or neutron diffraction experimental data. When a good agreement between theory and experiment is found, one can be sure that the molecular dynamics simulation provides a reliable and correct description of the structural properties of the system under

In particular, one can calculate the static structure factor either of neutrons or of X-ray (Keen, 2000) from the trajectory of the liquid phase simulation and directly compare it with the experimental one. For X-ray scattering data the radial distribution function between

where *ρj* is the number density of particles *j*, and *nij*(*r*) is the number of particles of type *j* between distance *r* and *r* + *dr* from a particle of type *i*. Given *cj* the molar fraction of particles *j*, *ρj* = *cj* · *ρ*0 where *ρ*0 is the total number density of the sample. The total

> ∞ 0 4*πr*

*nij*(*r*) 4*πr*2*drρj*

2[*gij*(*r*) − 1]

sin(*Qr*)

*Qr* (4)

(3)

*gij*(*r*) =

*cicj fij*(*Q*)

**6.3 Molecular dynamics**

calculate statistical averages.

analysis of the trajectory.

investigation.

of molecular modeling of ionic liquids (Maginn, 2009).

particle type *i* and *j* can be used, which are defined as:

scattering for X-ray can be then written as:

*F*(*Q*) = *ρ*0

*n* ∑ *i*=1 where *Q* id the scattering wave-vector modulus and *fi* are the atomic Q-dependent X-ray scattering factors and where

$$f\_{\vec{i}\vec{j}}(\mathbb{Q}) = \frac{f\_{\vec{i}}(\mathbb{Q})f\_{\vec{j}}(\mathbb{Q})}{[\sum\_{i=1}^{n}c\_{i}f\_{\vec{i}}(\mathbb{Q})]^{2}} \tag{5}$$

$$\mathcal{C}\_p \sim \left[ \frac{H(T + \Delta T) - H(T)}{\Delta T} \right]\_P \tag{6}$$

$$D\_i = \frac{1}{6} \lim\_{t \to \infty} \frac{\left\langle \left|| \mathbf{r}\_i(t) - \mathbf{r}\_i(0) \right||^2 \right\rangle}{t} \tag{7}$$

where **r***i*(*t*) is the position of the ion at time *t* while **r***i*(0) is the ion initial position.

Being a dynamic property, the diffusion coefficient can only be computed by molecular dynamics methods and not by Monte Carlo techniques. As we shall see, the difficult part for an ionic liquid is that the long enough simulation time can be of the order of 10-100 ns which is hardly attainable for current hardware and sophisticated molecular entities. Diffusive regimes where

$$\beta = \frac{\Delta \log(\Delta r^2)}{\Delta \log(t)} \sim 1 \tag{8}$$

can be obtained with some ionic liquids only after 10 ns of simulation time. In this case only very long trajectories can provide reliable information on the diffusion coefficients of the system. As we have mentioned above, non polarizable force fields, or those force fields which imply an integer charge on the ionic pairs constituents are often under-predicting diffusion coefficients to some extent. Polarizable force fields or charge scaling procedure yielded much better results.

Fig. 2. Schematic structure of the [C*n*dmim][NTf2] cations: *n* can be 3, 6, 9 and 12.

## **7. An example: dicationic imidazolium**

A special class of ionic liquids have recently been synthetized using geminal imidazolium dications (Anderson et al., 2005; Payagala et al., 2007; Sun et al., 2009). This class of compounds represents an interesting variation of the cationic partner and may present several advantages over the traditional monocationic ionic liquid in their application as lubricants (Jin et al., 2006), solvents (Han & Armstrong, 2005; Xiao & Shreeve, 2005) and separation media (Anderson et al., 2005).

We have recently performed a theoretical study of the gas phase structure of these ionic liquids (Bodo & Caminiti, 2010) because it provides:

A schematic view of the molecular structures that we have examined is reported in figure 2: we have a linkage chain (whose length can be 3, 6, 9 or 12 for the compounds analyzed) that connects two imidazolium rings with a net positive charge on them and that are substituted with a methyl group. Various ionic liquids have been obtained by mixing these dications with Br−, Tf2N−, BF− 4 and PF− 6 . These ionic liquids have thermal stabilities which are larger than their mono-cationic counterpart, their density decreases when increasing the linkage chain length (Anderson et al., 2005) while the viscosity increases. Except for those with the Tf2N anions, most of these kind of dicationic ionic liquids (especially those with short linkage chains) are solid at room temperature and melt at temperature higher than 100◦C (Anderson et al., 2005). Increasing the linkage chain length decreases the melting temperature and the C12 is a liquid at RT with all anion counterparts.

#### **7.1 The minimum energy structure problem**

When the number of atoms in a certain complex molecule or molecular system increases, the number of local minima of the potential energy surface (PES) quickly becomes so large that a systematic search is impossible. The study of stable geometries of large molecular systems characterized by a high conformational mobility, has to locate possible good candidates 12 Will-be-set-by-IN-TECH

+ +

+ + ní2

A special class of ionic liquids have recently been synthetized using geminal imidazolium dications (Anderson et al., 2005; Payagala et al., 2007; Sun et al., 2009). This class of compounds represents an interesting variation of the cationic partner and may present several advantages over the traditional monocationic ionic liquid in their application as lubricants (Jin et al., 2006), solvents (Han & Armstrong, 2005; Xiao & Shreeve, 2005) and separation media (Anderson

We have recently performed a theoretical study of the gas phase structure of these ionic liquids

• a series of model structures that may be present also in the liquid phase. The ionic interactions in these complexes are often strong enough to preserve the gas-phase structure also in the liquid phase(Logotheti et al., 2009) or at least make it one of the dominant

• the electronic density and the molecular electrostatic potential which in turn can be fed into the force fields for the molecular dynamics simulations of the liquid phase through a

• the charge distribution which is of great importance for the description the structural properties of the liquid phase (see for example refs. (Li & Kobrak, 2009; Spohr & Patey,

A schematic view of the molecular structures that we have examined is reported in figure 2: we have a linkage chain (whose length can be 3, 6, 9 or 12 for the compounds analyzed) that connects two imidazolium rings with a net positive charge on them and that are substituted with a methyl group. Various ionic liquids have been obtained by mixing these dications with

their mono-cationic counterpart, their density decreases when increasing the linkage chain length (Anderson et al., 2005) while the viscosity increases. Except for those with the Tf2N anions, most of these kind of dicationic ionic liquids (especially those with short linkage chains) are solid at room temperature and melt at temperature higher than 100◦C (Anderson et al., 2005). Increasing the linkage chain length decreases the melting temperature and the

When the number of atoms in a certain complex molecule or molecular system increases, the number of local minima of the potential energy surface (PES) quickly becomes so large that a systematic search is impossible. The study of stable geometries of large molecular systems characterized by a high conformational mobility, has to locate possible good candidates

6 . These ionic liquids have thermal stabilities which are larger than

Fig. 2. Schematic structure of the [C*n*dmim][NTf2] cations: *n* can be 3, 6, 9 and 12.

C1 C1 C1 C1 CS

ní2

2 1 3

45

CR NA NA

CW CW

**7. An example: dicationic imidazolium**

(Bodo & Caminiti, 2010) because it provides:

4 and PF−

C12 is a liquid at RT with all anion counterparts.

**7.1 The minimum energy structure problem**

fitting procedure generating partial atomic charges;

et al., 2005).

conformers;

2009)).

Br−, Tf2N−, BF−

for the time-consuming, *ab-initio* optimizations. One possible approach, in relatively small systems, is that of generating randomly (or by chemical intuition) many different starting structures and optimizing them separately. However, this route is impracticable in the present compounds because of the large number of possible initial realistic geometric guesses which are determined by the relative position of the chain and of the two anions.

We have employed a well tested force field for imidazolium based ionic liquids (Lopes et al., 2004) and we have used the *tinker* (Ren & Ponder, 2002) package in order to locate various minima on the PES. Our recipe is the following: we start with an isolated ionic complex made by the di-cation and two Tf2N anions in close contact but in a random orientation. We run a 100 ns molecular dynamics trajectory at 150-200 K with 1 fs timestep saving a snapshot of the system each 1000 steps. We then minimize each snapshot and locate the various minimum energy structures. We select two or three structures among the lowest energy ones and we repeat the dynamics for each of these structures but with a lower temperature. After a few cycles, we are generally able to locate a recurring structure that, very likely, is the global minima for the adopted force field.

The global minima determined above has then been optimized at the CAM-B3LYP/6-31G(d) (Peach et al., 2006; Yanai et al., 2004) level using the Gaussian program (Frisch et al., 2009). Obviously, given the high number of local minima, we cannot be sure whether we are describing the global minima over the potential energy surface of these complexes or a local one. However, the molecular dynamics-based annealing procedure described above should be certainly able to produce reasonable structures so that we can consider the ab-initio analysis to describe realistically at least the key structural feature of the complexes. One should however always keep in mind that small energetic variations are certainly possible when considering the various position that the linkage chain can assume especially when it is long enough to acquire a certain conformational mobility as for the C9 and C12 compounds.

#### **7.2 The gas phase structures**

Possible lowest energy structures for the linkage chains with 3, 6, 9 and 12 carbon atoms are reported in figure 3. The structures show clearly that the linkage chain "entangles" in order to maximize the interactions of the two imidazolium rings with the two anions. Crystals of [C2dmim][PF6] have been obtained (Lee et al., 2010) and the reported structure showed a parallel arrangement of the imidazolium rings. As far as we know no solid state structures have been reported for longer linkage chains.

In general, in all the four structures, there are at least two close contacts between the "acidic" C(2)-H hydrogens and the oxygens of the Tf2N anion (Kempter & Kirchner, 2010). For the C3 and for the C9 the distances we have found are 2.16 and 2.19 Å for the former and 2.08 and 2.16 Å for the latter. In the case of the C6 and of the C12 instead, we have 2.32, 2.77 Å and 2.22, 2.88 Å. Moreover in the C3 and C9 we have that the two imidazolium rings are parallel while this is not the case for the other two compounds. Whether these geometric variations are relevant or not in determining the local structure in the liquid phase is difficult to say (see below).

Another important indicator that can help in the interpretation of the liquid phase properties is the charge distribution in the ionic complex. In the entangled configurations the distribution is rather different to what is found, for example, in other mono-imidazolium based liquids. In those systems one finds a highly polar head (the ring) connected through electrostatic interactions (and possibly through a weak, special hydrogen bond) to the anion and a non polar alkyl chain which is free to move creating an apolar tail with high rotational mobility. Here we have a rather spherical and compact structure with no or very low conformational

Fig. 3. Pictorial representation of the lowest energy structures associated with the 4 compounds. From the upper left corner, we have: C3, C6, C9 and C12. Hydrogens and Fluorine atoms have been omitted for clarity

mobility in which the charge is distributed in a quadrupolar way: two negatively charged region are counterbalanced by the two imidazolium rings and a rather apolar domain is created by the bent alkyl chain. We present the maps of the electrostatic potential for the C9 compounds in 4. The positive charge is localized on the rings and the negative one on the oxygens of the Tf2N.

We have run molecular dynamics liquid phase simulations using the above force field. This has been done in order to check if the compact structure survives in the liquid phase. What we have found is that the structure actually survives for the C3 linkage chain, but is almost lost for the C12 one. The other two compounds present a somewhat mixed behavior. Anyway, given the lack of experimental structural studies, we have no validation of the performed dynamics and we cannot say if the adopted force field (especially the atomic charges) would lead to a realistic local liquid structure. Here we report just few of the results we have gathered on such systems. The molecular dynamics protocol we have followed is very similar to that used in our previous publication (Bodo et al., 2010). All-atom molecular dynamics simulations have been carried out on pure [C*n*dmim][NTf2] ionic liquids with *n* =3, 6, 9 or 12 using the DLPOLY2 14 Will-be-set-by-IN-TECH

Fig. 3. Pictorial representation of the lowest energy structures associated with the 4 compounds. From the upper left corner, we have: C3, C6, C9 and C12. Hydrogens and

mobility in which the charge is distributed in a quadrupolar way: two negatively charged region are counterbalanced by the two imidazolium rings and a rather apolar domain is created by the bent alkyl chain. We present the maps of the electrostatic potential for the C9 compounds in 4. The positive charge is localized on the rings and the negative one on the

We have run molecular dynamics liquid phase simulations using the above force field. This has been done in order to check if the compact structure survives in the liquid phase. What we have found is that the structure actually survives for the C3 linkage chain, but is almost lost for the C12 one. The other two compounds present a somewhat mixed behavior. Anyway, given the lack of experimental structural studies, we have no validation of the performed dynamics and we cannot say if the adopted force field (especially the atomic charges) would lead to a realistic local liquid structure. Here we report just few of the results we have gathered on such systems. The molecular dynamics protocol we have followed is very similar to that used in our previous publication (Bodo et al., 2010). All-atom molecular dynamics simulations have been carried out on pure [C*n*dmim][NTf2] ionic liquids with *n* =3, 6, 9 or 12 using the DLPOLY2

Fluorine atoms have been omitted for clarity

oxygens of the Tf2N.

Fig. 4. Charge distribution mapped on the electronic density for the structure depicted on the right for the C9dmim(Tf2N)2 ionic complex. Red regions are charged negatively, blue ones are charged positively.

(Smith et al., 2002) package. We have employed the force field from Canongia-Lopes and Pádua (Canongia Lopes et al., 2005; Lopes et al., 2004; Lopes & Pádua, 2004).

The initial configurations have been generated by randomly distributing 300 cations and 600 anions in very large simulation cells. A short (100 ps) isobaric-isothermal (NPT) run at 500 K and with 103 atmospheres pressure was carried out in order to provide densely-packed, mixed initial configurations. These initial configurations were left to relax during another NPT run of 300-500 ps at 1 atm and 300 K up to the point in which an approximately constant density was obtained. The densities we have obtained are within few percent from the experimental ones (Anderson et al., 2005). The volume was then fixed and an additional mixing was provided by a short (300 ps) NVT run at 300 K where the electrostatic interactions were turned off. The resulting configurations were finally left to evolve in a constant volume NVT production runs of about 5 ns with T = 300 K after an equilibration period.

Apart from the densities the other experimental quantities available are solid liquid transition temperatures and viscosities which are extremely expensive from the computational point of view and therefore we have not computed them.

However the agreement on density allows us to trace a few geometric considerations. One important feature of the liquid state is obviously represented by the position of the linkage chain: we have calculated the *gij*(*r*) for the two lateral "C1" carbon atoms and we report it in figure 5 for the [C3dmim][Tf2N] and [C9dmim][Tf2N] compounds. We report in the same figure also the distances between the two lateral C1 carbon atoms in the ionic complex with two anions (smaller value) and in the isolated cation (larger value). As can be seen in both moieties the linkage in the liquid phase (at least in our simulations) seems to present a somewhat bent configuration which is not as compact as in the isolated ionic complex, nor extended as in the isolated cation. In any case the most abundant geometric conformation is the one in which the linkage chain is almost extended. Polarization effect might play an

Fig. 5. Lateral C1 self radial distribution for the [C3dmim][Tf2N] and [C9dmim][Tf2N] compounds. The vertical lines identify the C1-C1 distances attained in the gas by the ionic complex with two Tf2N or in the isolated cation.

important role in determining a different degree of entanglement of the complexes, but we were unable to verify it.

Another interesting feature of these ionic liquids is the relative position of the two imidazolium rings of the same cation. By looking at the gas phase structures above, one can see that some of the complexes show a parallel arrangement of the two imidazolium rings. We have calculated the angular distribution between imidazolium rings of the same molecule and we report it in Figure 6 for the [C3dmim][Tf2N] and [C9dmim][Tf2N] compounds. For the C3 compound we find that the two imidazolium rings prefer a non parallel configuration in which the angle between them is around 50 degrees. We see, however, that a sizable fraction of molecules presents an angle near to zero, i.e. a parallel arrangement of the two rings.

#### **8. Conclusions**

In this Chapter we have reviewed the typical theoretical approaches that are commonly used to describe ionic liquids using *in silico* simulations. The study of these materials requires a combined approach made by ab-initio calculations and MD simulations. These techniques require a special attention in the application to ionic liquids because they of the intrinsic difficulties due to the extremely slow dynamics in the real systems. We have concluded our chapter by reporting few highlights from a recent work that we are carrying out on a special class of these compounds and we have pointed out a few interesting features we have discovered in order to provide an example of how powerful theoretical predictions can be for molecular systems which are still relatively unexplored experimentally.

16 Will-be-set-by-IN-TECH

[C9

[C3

mim][Tf2

mim][Tf2

N]

N]

024 6 8 10 12 14 16 18 20 R / Å

Fig. 5. Lateral C1 self radial distribution for the [C3dmim][Tf2N] and [C9dmim][Tf2N] compounds. The vertical lines identify the C1-C1 distances attained in the gas by the ionic

important role in determining a different degree of entanglement of the complexes, but we

Another interesting feature of these ionic liquids is the relative position of the two imidazolium rings of the same cation. By looking at the gas phase structures above, one can see that some of the complexes show a parallel arrangement of the two imidazolium rings. We have calculated the angular distribution between imidazolium rings of the same molecule and we report it in Figure 6 for the [C3dmim][Tf2N] and [C9dmim][Tf2N] compounds. For the C3 compound we find that the two imidazolium rings prefer a non parallel configuration in which the angle between them is around 50 degrees. We see, however, that a sizable fraction of molecules presents an angle near to zero, i.e. a parallel arrangement of the two rings.

In this Chapter we have reviewed the typical theoretical approaches that are commonly used to describe ionic liquids using *in silico* simulations. The study of these materials requires a combined approach made by ab-initio calculations and MD simulations. These techniques require a special attention in the application to ionic liquids because they of the intrinsic difficulties due to the extremely slow dynamics in the real systems. We have concluded our chapter by reporting few highlights from a recent work that we are carrying out on a special class of these compounds and we have pointed out a few interesting features we have discovered in order to provide an example of how powerful theoretical predictions can be for

molecular systems which are still relatively unexplored experimentally.

0

were unable to verify it.

**8. Conclusions**

complex with two Tf2N or in the isolated cation.

0.005

0.01

0.015

g(r) / arb. units

0.02

0.025

0.03

Fig. 6. Angle distribution for the two imidazolium rings of the same cation for the [C3dmim][Tf2N] and [C9dmim][Tf2N] compounds

#### **9. References**

18 Will-be-set-by-IN-TECH

Fei, Z., Geldbach, T. J., Zhao, D. & Dyson, P. J. (2006). From dysfunction to bis-function:

Fischer, T., Sethi, A., Welton, T. & Woolf, J. (1999). Diels-alder reactions in room-temperature

Forsyth, S. A., Pringle, J. M. & MacFarlane, D. R. (2004). Ionic liquids-an overview., *Aust. J.*

Frisch, M. J., Trucks, G. W., Schlegel, H. B., Scuseria, G. E., Robb, M. A., Cheeseman, J. R.,

Fujita, K., Nakamura, N., Igarashi, K., Samejima, M. & Ohno, H. (2009). Biocatalytic oxidation

Galinski, M., Lewandowski, A. & Stepniak, I. (2006). Ionic liquids as electrolytes., *Electrochim.*

Garcia, B., Lavallee, S., Perron, G., Michot, C. & Armand, M. (2004). Room temperature molten

Gaune-Escard, M. & Seddon, K. R. (eds) (2010). *Molten Salts and Ionic Liquids: Never the Twain?*,

Han, X. & Armstrong, D. W. (2005). Using geminal dicationic ionic liquids as solvents for

Hough, W. L. & Rogers, R. D. (2007). Ionic liquids then and now: from solvents to materials to active pharmaceutical ingredients., *Bull. Chem. Soc. Jpn.* 80: 2262–2269. Hough, W. L., Smiglak, M., Rodriguez, H., Swatloski, R. P., Spear, S. K., Daly, D. T., Pernak,

Hu, Z. & Margulis, C. J. (2007). Room-temperature ionic liquids: Slow dynamics, viscosity,

Jin, C.-M., Ye, C., Phillips, B. S., Zabinski, J. S., Liu, X., Liu, W. & Shreeve, J. M.

Jorgensen, W. L., Maxwell, D. S. & Tirado-Rives, J. (1996). Development and Testing of the

J., Grisel, J. E., Carliss, R. D., Soutullo, M. D., Davis, Jr., J. H. & Rogers, R. D. (2007). The third evolution of ionic liquids: Active pharmaceutical ingredients., *New J. Chem.*

(2006). Polyethylene glycol functionalized dicationic ionic liquids with alkyl or polyfluoroalkyl substituents as high temperature lubricants, *J. Mater. Chem.* 16: 1529.

OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic

of cellobiose in an hydrated ionic liquid, *Green Chem.* 11: 351–354.

salts as lithium battery electrolyte., *Electrochim. Acta* 49: 4583–4588.

high-temperature organic reactions, *Org. Lett.* 7: 4205.

and the red edge effect., *Acc. Chem. Res.* 40: 1097–1105.

Liquids, *J. Am. Chem. Soc.* 118: 11225–11236.

Hanke, C., Price, S. & Lynden-Bell, R. (2001). *Mol. Phys.* 99: 801.

12: 2122–2130.

*Chem.* 57: 113–119.

Freemantle, M. (1998). *Chem . Eng. News.* 76: 32.

Gaussian Inc. Wallingford CT 2009.

*Acta* 51: 5567–5580.

Wiley.

31: 1429–1436.

ionic liquids., *Tetrahedron Lett.* 40: 793–796.

Freemantle, M. (2009). *An Introduction to Ionic Liquids*, RSC Publishing.

## **Classical Density Functional Theory of Ionic Liquids**

Jan Forsman1, Ryan Szparaga1, Sture Nordholm2,

Clifford E.Woodward3 and Robert Penfold4 *1Lund University 2Gothenburg University 3University of New South Wales, ADFA 4Institute of Food Research 1,2Sweden 3Australia 4UK*

#### **1. Introduction**

20 Will-be-set-by-IN-TECH

126 Ionic Liquids – Classes and Properties

Rogers, R. D., Plechkova, N. V. & Seddon, K. R. (eds) (2009). *Ionic Liquids: From Knowledge*

Rogers, R. D. & Seddon, K. R. (2003). Ionic liquids–solvents of the future?, *Science*

Rogers, R. & Seddon, K. (eds) (2005a). *Ionic Liquids IIIA: Fundamentals, Progress, Challenges, and*

Rogers, R. & Seddon, K. (eds) (2005b). *Ionic Liquids IIIB: Fundamentals, Progress, Challenges,*

Smith, W., Yong, C. W. & Rodger, P. M. (2002). Dl poly: Application to molecular simulation,

Spohr, H. V. & Patey, G. N. (2009). Structural and dynamical properties of ionic liquids: The

Sun, H., Zhang, D., Liu, C. & Zhang, C. (2009). Geometrical and electronic

Wang, P., Zakeeruddin, S. M., Exnar, I. & Gratzel, M. (2002). High efficiency dye-sensitized

Wang, Y. & Voth, G. A. (2005). Unique spatial heterogeneity in ionic liquids., *J. Am. Chem. Soc.*

Wang, Y. & Voth, G. A. (2006). Tail aggregation and domain diffusion in ionic liquids., *J. Phys.*

Welton, T. & Wasserscheid, P. (eds) (2008). *Ionic Liq. Synth.*, 2nd edn, Wiley VCH, Weinheim. Wilkes, J. S. (2002). A short history of ionic liquids-from molten salts to neoteric solvents,

Wilkes, J. S. & Zaworotko, M. J. (1992). Air and water stable 1-ethyl-3-methylimidazolium

Wolfbauer, G., Bond, A., Eklund, J. C. & Macfarlane, D. (2001). *Sol. Energy Mater. Sol. Cells*

Xiao, J.-C. & Shreeve, J. M. (2005). Synthesis of 2,2'-biimidazolium-based ionic liquids: Use

Xu, W. & Angell, C. A. (2003). Solvent-free electrolytes with aqueous solution-like

Yanai, T., Tew, D. & Handy, N. (2004). A new hybrid exchange-correlation functional using the coulomb-attenuating method (cam-b3lyp), *Chem. Phys. Lett.* 393: 51. Zahn, S. & Kirchner, B. (2008). Validation of dispersion corrected functional theory for ionic

Zahn, S., Uhlig, F., Thar, J., Spickermann, C. & Kirchner, B. (2008). Intermolecular forces in

as a new reaction medium and ligand for palladium-catalyzed suzuki cross-coupling

an ionic liquid ([mmim][cl]) versus those in a typical salt (nacl), *Angewandte Chemie*

based ionic liquids, *J. Chem. Soc., Chem. Commun.* pp. 965–967.

13-bis[3-methylimidazolium]propane bromide, *Theochem* 900: 37.

influence of charge location, *J. Chem. Phys.* 130: 104506.

D.C.

302(5646): 792–793.

Society, Washington D.C.

Walden., P. (1914b). *Chem. Zentralbl.* p. 1800.

*Commun.* pp. 2972–2973.

*Chem. B* 110: 18601–18608.

reactions, *J. Org. Chem.* 70: 3072.

conductivities, *Science* 302(5644): 422–425.

lquid systems, *J. Phys. Chem. A* 112: 8430.

*International Edition* 47(19): 3639–3641.

*Green Chem.* 4: 73–80.

70: 85–101.

127: 12192–12193.

Chemical Society, Washington D.C.

*Molecular Simulation* 28: 892–7022.

Walden, P. (1914a). *Bull. Acad. Imp. Sci. St.-Petersbourg* 8: 405.

*to Application*, Vol. 1030 of *ACS Symp. Ser.*, American Chemical Society, Washington

*Opportunities: Properties and Structure*, Vol. 901 of *ACS Symp. Ser.*, American Chemical

*and Opportunities: Transformations and Processes*, Vol. 902 of *ACS Symp. Ser.*, American

structures of the dication and ion pair in the geminal dicationic ionic liquid

nanocrystalline solar cells based on ionic liquid polymer gel electrolyte, *Chem.*

The current technological revolution based on ionic liquids (ILs) is driven by their unique properties. Though ILs have been known to scientists for close to a century (Walden , 1914), they are now the focus of intense activity mainly because of their promise for "environmentally friendly" applications. Some actual and potential uses of ILs include: specific solvents for heterogeneous and homogeneous catalysis; selective solvents for removal of heavy metal contaminants; electrolytes in various electrochemical processes and devices; and as dispersive agents for stabilization of nanoparticles. In all of these applications, the structural properties of ILs and their mixtures in the bulk phase and at interfaces are crucial to their performance (Maier et al. , 2010). In electrowinning processes, such as aluminium refining, metal ion speciation is a crucial consideration in determining the population of electroactive components (Rocher et al. , 2009) and consequently influences the pertinent dynamical processes that lead to oxidation and reduction. These populations are ultimately governed by the free energy of the IL solution, which in turn depends upon the mutual interactions between the IL molecules and the added components.

It is now clear that fundamental theoretical studies of structure-property relationships in ILs will provide important new insights, which will assist this kind of research. Not only in the choice and design of ILs for a specific application, but also for scenario based modelling of IL processes. For example, theoretical studies of the structure of ILs at charged surfaces are important for understanding the properties that they impart to electrodeposition processes, including factors such as ionic speciation and mass transport of electroactive species (MacFarlane et al. , 2010). As emphasized by Kornyshev in a recent review article (Kornyshev , 2007), applications of ILs at electrified interfaces to energy-storage systems, electrowetting devices or nanojunction gating media will be strongly promoted by a deeper understanding of the structure and properties of the interfacial double layer.

The importance of theoretical investigations of ILs is also manifested in the rapid accumulation of such studies in the scientific literature. This includes a large number of fully atomistic simulations of ILs, which first appeared some ten years ago (Hanke et al. , 2001; Lynden-Bell , 2010; Lynden-Bell et al. , 2007). These reports have provided IL researchers with valuable physical and chemical information with which to rationalize many experimental observations. By deploying numerically intensive Monte Carlo or Molecular Dynamics algorithms, computer simulations circumvent the technical difficulties inherent to dense liquids with strong electrostatic coupling that have hampered the development of a fundamental theory of ILs based upon statistical mechanical formalism. On the other hand, a fundamental theory of this type does exist for traditional aqueous electrolytes wherein ionic concentrations are generally low and electrostatic coupling is much weaker than in ILs. In particular, a simple mean-field theory such as the Poisson-Boltzmann Approximation (PBA) (Gouy , 1910) can be successfully applied. Calculations using the PBA are computationally rapid and reasonably accurate for traditional electrolytes. Even though seminal work by Jönsson and coworkers (Gulbrand et al. , 1984) has shown that the PBA can fail significantly for divalent ions, the PBA nevertheless remains a useful cornerstone theory for aqueous electrolytes and has provided much of the conceptual understanding in those systems.

Conditions that are found in typical ILs would normally cause a simple electrolyte to crystallize. Ionic liquids remain fluid, however, due to steric hindrance accompanying their larger and more complex molecular shapes. Thus, in ILs the steric length scale becomes relevant and there is strong interplay between electrostatic and steric interactions. Under these conditions the PBA fails. This has left something of a conceptual gap in the theoretical treatment of ILs which has been largely filled with detailed all-atom simulations. Unfortunately, these simulations require large computational resources and are not practical for rapid, exploratory calculations. Moreover, the simulation results may be difficult to interpret since important mechanisms may be hidden in the complex numerical calculations. In this chapter we describe an alternative theoretical approach that aims to play a similar role for ILs as the PBA does for traditional electrolytes. As with the PBA, this theory will be mathematically concise (albeit approximate) but contain sufficient corrections to the PBA to render it applicable for ILs. In order to make analytical progress, it will be necessary to use simpler coarse-grained models of ILs than those currently implemented in all-atom simulations. Though this will lead to an inevitable loss of detailed information, we anticipate that this will be outweighed by the value of an approximate but mechanistically more explicit and numerically less demanding theory for ILs. The advantages of such a theory to researchers in the field are both practical and technical:

#### **2. Other approaches**

2 Will-be-set-by-IN-TECH

The importance of theoretical investigations of ILs is also manifested in the rapid accumulation of such studies in the scientific literature. This includes a large number of fully atomistic simulations of ILs, which first appeared some ten years ago (Hanke et al. , 2001; Lynden-Bell , 2010; Lynden-Bell et al. , 2007). These reports have provided IL researchers with valuable physical and chemical information with which to rationalize many experimental observations. By deploying numerically intensive Monte Carlo or Molecular Dynamics algorithms, computer simulations circumvent the technical difficulties inherent to dense liquids with strong electrostatic coupling that have hampered the development of a fundamental theory of ILs based upon statistical mechanical formalism. On the other hand, a fundamental theory of this type does exist for traditional aqueous electrolytes wherein ionic concentrations are generally low and electrostatic coupling is much weaker than in ILs. In particular, a simple mean-field theory such as the Poisson-Boltzmann Approximation (PBA) (Gouy , 1910) can be successfully applied. Calculations using the PBA are computationally rapid and reasonably accurate for traditional electrolytes. Even though seminal work by Jönsson and coworkers (Gulbrand et al. , 1984) has shown that the PBA can fail significantly for divalent ions, the PBA nevertheless remains a useful cornerstone theory for aqueous electrolytes and has provided much of the conceptual understanding in those

Conditions that are found in typical ILs would normally cause a simple electrolyte to crystallize. Ionic liquids remain fluid, however, due to steric hindrance accompanying their larger and more complex molecular shapes. Thus, in ILs the steric length scale becomes relevant and there is strong interplay between electrostatic and steric interactions. Under these conditions the PBA fails. This has left something of a conceptual gap in the theoretical treatment of ILs which has been largely filled with detailed all-atom simulations. Unfortunately, these simulations require large computational resources and are not practical for rapid, exploratory calculations. Moreover, the simulation results may be difficult to interpret since important mechanisms may be hidden in the complex numerical calculations. In this chapter we describe an alternative theoretical approach that aims to play a similar role for ILs as the PBA does for traditional electrolytes. As with the PBA, this theory will be mathematically concise (albeit approximate) but contain sufficient corrections to the PBA to render it applicable for ILs. In order to make analytical progress, it will be necessary to use simpler coarse-grained models of ILs than those currently implemented in all-atom simulations. Though this will lead to an inevitable loss of detailed information, we anticipate that this will be outweighed by the value of an approximate but mechanistically more explicit and numerically less demanding theory for ILs. The advantages of such a theory

• the ability to do rapid calculations will facilitate more ambitious large-scale scenario modelling, which allows us to reasonably cheaply explore a wide range of systems,

• mathematically concise, closed expressions for entropy related quantities (which are difficult to obtain in simulations) will allow identification of general and powerful mechanistic principles, that can be used to describe and organise large classes of

to researchers in the field are both practical and technical:

inaccessible to all-atom simulations;

phenomena in ILs.

systems.

ILs are fascinating from a theoretical standpoint. In order to treat ILs, a theory must account for both electrostatic and steric correlations. Attempts to develop such a theory have appeared in the literature. These include the lattice based model by Kornyshev and co-workers (Fedorov et al. , 2008; Kornyshev , 2007). In another lattice approach by Lauw and coworkers (Lauw et al. , 2009), both cations and anions were modelled by a star-like architecture, with each monomer subject to an independently specified dielectric response. Under certain circumstances, these models were able to qualitatively reproduce the so-called "camel-shaped" differential capacitance, *CD*, profile observed in imidazolium ILs at carbon electrodes (Islam et al. , 2009; Lockett et al. , 2008a). This behaviour is an indication of subtle structuring effects at charged surfaces and is a stern test for theory. Although camel-shaped *CD* profiles are predicted, qualitative discrepancy against experimental results remains. This is not surprising since these crude theories lack a proper description of electrostatic correlations despite accounting for steric contributions. Charge-charge correlations have been shown to be important in interfacial phenomena of highly coupled electrostatic fluids. As with dispersion forces, ionic correlations are always attractive and lower the total energy. Thus, we expect that they will also play a role in establishing the camel shape of the *CD*. That this is the case, is seen from both simulations (Lamperski & Zydor , 2007; Vatamanu et al. , 2010) and modified Poisson-Boltzmann treatments (Lamperski et al. , 2009) of a restricted primitive model (RPM) electrolyte. These investigations show that the RPM displays a camel-shaped *CD* at low temperatures and densities. However, electrostatic correlations do not tell the whole story for ILs. That is, for the RPM the *CD* becomes bell-shaped and exhibits a single maximum at densities and temperatures typical for ILs. This suggests that the simple RPM is an inadequate model for ILs and other factors, such as dispersion forces and molecular structure (Fedorov et al. , 2010; Trulsson et al. , 2010), play an important role in determining their properties. In this work, we will use a model for ILs wherein the molecular ions interact via dispersion, as well as electrostatic, forces and the cation has an oligomeric structure consistent with imidazolium ions. This model has also been adopted by us in previous recent publications (Forsman et al. , 2011; Trulsson et al. , 2010), and employs the methods of polymer free energy density functional theory (DFT) (Woodward , 1991). The result will be a simple, albeit approximate DFT for ILs. We have already introduced the salient features of the theory in a previous publication (Forsman et al. , 2011). Here, we will present a much more thorough treatment, emphasizing the relation to theories of simpler fluids. We will also extend the scope of the study with some emphasis on temperature dependence and the effects of ion-specific adsorption. The IL-DFT that we will discuss here is an off-lattice theory, whose accuracy can be straightforwardly tested against computer simulations of the same molecular model. We will begin with a discussion of relevant earlier work, in order to introduce the essential mechanistic components that will be combined in our density functional theory for ionic liquids.

#### **3. Classical density functional theory (DFT) of simple fluids**

The ideal gas law is known and used by most physicists and chemists but a more complete description of a fluid including its ability to condense into a liquid requires a treatment of the interactions between its particles. The simplest realistic model of a molecular fluid assumes pairwise spherically symmetric interactions by a pair-potential *φ*(*r*) which has a repulsive (positive) inner region characterised by a length scale *σ* and an attractive tail that defines an energy scale *�*. Such models are referred to as simple fluids. The most popular form of simple fluid pair-potential is the Lennard-Jones potential

$$
\phi(r) = 4\epsilon \left( \left( \frac{\sigma}{r} \right)^{12} - \left( \frac{\sigma}{r} \right)^{6} \right) \,\,\,. \tag{1}
$$

We shall let this familiar pair-potential illustrate our simple fluid in the discussion below. Even for a simple fluid, the the analysis quickly becomes quite complicated if we attempt to describe configurational states by explicitly locating all the particles as in simulation methods, or if we try to resolve the coupled hierarchy of correlation functions as in rigorous liquid state theory. Here we shall instead consider an approach based on deliberate and efficient approximations seeking to represent the essential features while hiding the less important details of the fluid properties. The most noteworthy early attempt in this direction was taken by J. van der Waals (van der Waals , 1873) who 140 years ago realized that the interactions in molecular fluids introduced two mechanisms: i) an excluded volume effect due to the size of the molecules and ii) a binding energy due to the attraction between molecules. His equation of state relating pressure *P* to volume per particle *v* and absolute temperature *T* reads

$$P = \frac{k\_B T}{v - v\_0} - \frac{a}{v^2} \ . \tag{2}$$

It incorporates an excluded volume *v*0 per particle and a binding energy per particle *eb* = *a*/*v* due to the attractive interactions. Following his reasoning in the simplest way for the Lennard-Jones potential one would assume that particles do not approach closer than *r* = *σ* and therefore exclude on the average the volume *v*0 = 1 2 × 4*πσ*3/3 = 2*πσ*3/3 per particle as evaluated in the low density limit. The attractive part of the potential is in van der Waals (vdW) theory assumed to be sampled without inducing any correlations, i.e. in the mean-field approximation so that the binding energy parameter becomes:

$$a = \frac{1}{2} \int\_{\sigma}^{\infty} dr \, 4\pi r^2 \phi(r) = \frac{16\pi}{9} \epsilon \sigma^3 \,\,\,\,\,\tag{3}$$

for the Lennard-Jones interaction. The corresponding configurational free energy per particle can be taken to be *fc*(*v*) = −*kBT* ln(*v* − *v*0) − *eb* = −*kBT* ln *v*free − *eb*, where *v*free is the free volume available to each particle after the excluded volume effect has been accounted for, i.e. *v*free = *V*/*N* − *v*0.

The origin of density functional theory can be traced back to the famous failure of the vdW equation of state in its prediction of unphysical "wiggles" in *P*(*v*, *T*) where the pressure could increase with increasing volume and even become negative. The Maxwell reconstruction (Maxwell , 1875) eliminates these wiggles and is based on interpreting the vdW equation of state as the description of a uniform fluid with density *n*(**r**) = *nb* = 1/*v* everywhere in the available volume. The total configurational free energy can then be obtained as

$$\mathcal{F}\_{\mathbf{c}}[n(\mathbf{r})] = \mathcal{F}\_{\mathbf{c}}(n\_{b}) = \int\_{V} d\mathbf{r} \, n(\mathbf{r}) f\_{\mathbf{c}}(\mathbf{r}) = \int\_{V} d\mathbf{r} \, n(\mathbf{r}) \left( k\_{B} T \ln \left( \frac{n(\mathbf{r})}{1 - n(\mathbf{r}) v\_{0}} \right) - a n(\mathbf{r}) \right)$$

$$= N \left( k\_{B} T \ln \left( \frac{n\_{b}}{1 - n\_{b} v\_{0}} \right) - a n\_{b} \right) \,. \tag{4}$$

Now the Maxwell reconstruction can be understood by recognising that the uniform fluid is not the stable state. Instead the fluid separates into a rarefied gas phase (of density *n*(*g*) *b* and 4 Will-be-set-by-IN-TECH

energy scale *�*. Such models are referred to as simple fluids. The most popular form of simple

. (1)

*v*2 . (2)

2 × 4*πσ*3/3 = 2*πσ*3/3 per particle as

9 *�σ*3 , (3)

 *σ r* 12 − *σ r* 6

We shall let this familiar pair-potential illustrate our simple fluid in the discussion below. Even for a simple fluid, the the analysis quickly becomes quite complicated if we attempt to describe configurational states by explicitly locating all the particles as in simulation methods, or if we try to resolve the coupled hierarchy of correlation functions as in rigorous liquid state theory. Here we shall instead consider an approach based on deliberate and efficient approximations seeking to represent the essential features while hiding the less important details of the fluid properties. The most noteworthy early attempt in this direction was taken by J. van der Waals (van der Waals , 1873) who 140 years ago realized that the interactions in molecular fluids introduced two mechanisms: i) an excluded volume effect due to the size of the molecules and ii) a binding energy due to the attraction between molecules. His equation of state relating pressure *P* to volume per particle *v* and absolute temperature *T* reads

> *P* = *kBT v* − *v*0

It incorporates an excluded volume *v*0 per particle and a binding energy per particle *eb* = *a*/*v* due to the attractive interactions. Following his reasoning in the simplest way for the Lennard-Jones potential one would assume that particles do not approach closer than *r* = *σ*

evaluated in the low density limit. The attractive part of the potential is in van der Waals (vdW) theory assumed to be sampled without inducing any correlations, i.e. in the mean-field

for the Lennard-Jones interaction. The corresponding configurational free energy per particle can be taken to be *fc*(*v*) = −*kBT* ln(*v* − *v*0) − *eb* = −*kBT* ln *v*free − *eb*, where *v*free is the free volume available to each particle after the excluded volume effect has been accounted for, i.e.

The origin of density functional theory can be traced back to the famous failure of the vdW equation of state in its prediction of unphysical "wiggles" in *P*(*v*, *T*) where the pressure could increase with increasing volume and even become negative. The Maxwell reconstruction (Maxwell , 1875) eliminates these wiggles and is based on interpreting the vdW equation of state as the description of a uniform fluid with density *n*(**r**) = *nb* = 1/*v* everywhere in the available volume. The total configurational free energy can then be

> *V*

Now the Maxwell reconstruction can be understood by recognising that the uniform fluid is not the stable state. Instead the fluid separates into a rarefied gas phase (of density *n*(*g*)

*d***r** *n*(**r**)

 − *anb*

 *kBT* ln  *n*(**r**) 1 − *n*(**r**)*v*0 . (4)

− *an*(**r**)

*b* and

*dr* 4*πr*

− *a*

2*φ*(*r*) = 16*π*

*φ*(*r*) = 4*�*

fluid pair-potential is the Lennard-Jones potential

and therefore exclude on the average the volume *v*0 = 1

*v*free = *V*/*N* − *v*0.

obtained as

F*c*[*n*(**r**)] = F*c*

 *nb* = *V*

> = *N kBT* ln

approximation so that the binding energy parameter becomes:

*a* = 1 2 ∞ *σ*

*d***r** *n*(**r**)*fc*(**r**) =

 *nb* 1 − *nbv*0 volume *V*(*g*)) coexisting with a condensed liquid phase of density *n*(*l*) *b > n*(*g*) *b* and volume *V*(*l*) *< V*(*g*). Since the total volume and particle number are conserved we have the two relations: *V*(*g*) + *V*(*l*) = *V* and *N*(*g*) + *N*(*l*) = *n*(*g*) *b V*(*g*) + *n*(*l*) *b V*(*l*) = *N*, which constrain the search for the most stable form of the fluid. At fixed *T*, minimizing the configurational free energy with respect to the choice of *n*(*g*) *b* and *n*(*l*) *b* will then produce the familiar Maxwell reconstructed isotherms with flat segments (tielines) replacing the wiggles.

#### **3.1 The generalized van der Waals (GvdW) theory**

The analysis of the traditional van der Waals theory above shows that:

In the treatment of bulk multiphase fluids as above one can assume that the variation of the density is slow on the length scale of either the repulsive or attractive range of interaction. This is certainly true well inside each phase but not neccessarily true in the interface between phases. Such interfaces are generally of minor extent (*N*2/3) compared to the bulk phases (*N*) so that their effect can be neglected in the thermodynamic limit. We may also have slow density variations due to external fields. An example is our atmosphere which is held in place by the gravitational field. We could describe the atmosphere very well with a vdW functional as above parametrized by excluded volume and binding energy.

We now consider two improvements of the vdW functional. First we would like to improve the accuracy of the vdW equation of state for simple fluids over their full density range. In liquids, it turns out that the excluded volume per particle *v*0 is severely overestimated by its low density limit value 2*πσ*3/3 as used above. Consideration of *v*0 as an estimate of the volume per particle in the close packing limit shows that this volume would be just *σ*3 in the simple cubic packing arrangement where each hard sphere of diameter *σ* is assigned a cubic volume in which it fits touching all sides. Thus, setting *v*0 = *σ*3 will improve the ability of the vdW equation of state to describe the full range of fluid densities.

Second, we shall restrict the assumption of local interactions to the short range repulsions (i.e., the excluded volume effect) but give the attractive interactions a nonlocal representation. The GvdW functional for our Lennard-Jones fluid becomes (Nordholm & Haymet , 1980):

$$\mathcal{F}\_{\mathcal{C}}[n(\mathbf{r})] = k\_B T \int\_V d\mathbf{r} \, n(\mathbf{r}) \ln \left( \frac{n(\mathbf{r})}{1 - n(\mathbf{r}) \sigma^3} \right) + \int\_V d\mathbf{r} \, n(\mathbf{r}) \left( \frac{1}{2} \int\_V d\mathbf{r}' \, n(\mathbf{r}') \phi\_{\rm diff}(\mathbf{r}, \mathbf{r}') + V\_{\rm ext}(\mathbf{r}) \right) \, . \tag{5}$$

Here the attractive part of the pair-potential is defined as

$$\phi\_{att}(\mathbf{r}, \mathbf{r}') = \begin{cases} 0 & |\mathbf{r} - \mathbf{r}'| \le \sigma \\ \phi(|\mathbf{r} - \mathbf{r}'|) & |\mathbf{r} - \mathbf{r}'| > \sigma \end{cases} \tag{6}$$

and we have added an energy due to an external field *V*ext(**r**). This functional can treat bulk phase equilibria as well as the main mechanisms of adsorption at a wall and gas-liquid surface tension in a simple fluid of this type (Johnson & Nordholm , 1981; Nordholm & Gibson , 1981).

#### **3.1.1 The fine-grained GvdW theory**

The key to many phenomena in nonuniform fluids is the nonlocality of interactions, that is the nonvanishing range of interaction on the length scale of density variations. From complete neglect of this nonlocality in the vdW analysis of phase separation we have proceeded above to account for the nonlocality of the pairwise attractive interaction leaving only the shorter range repulsive interactions to local approximation. There are, however, many applications where the nonlocality of the repulsive interactions is also important. We know, for example, that at higher densities hard spheres pack into distinct layers at hard walls. Such packing structures are not predicted by theories based on local representation of excluded volume effects. Moreover, it is clear that repulsive nonlocality will enter also into an accurate prediction of gas-liquid surface tension. For this reason fine-grained GvdW theory has been developed to allow both attractive and repulsive interactions to be treated in a nonlocal representation. The key step is to realize that the excluded volume effect is related to an average particle density over a spatial domain determined by the range of the repulsive interactions. If we, as above, treat the particles as interacting hard spheres of diameter *σ* then the simplest coarse-graining average reflecting the domain of the repulsive interactions is

$$
\overline{n}(\mathbf{r}) = \frac{3}{4\pi\sigma^3} \int\_{|\mathbf{r}-\mathbf{r}'| < \sigma} d\mathbf{r}' \, n(\mathbf{r}') \, . \tag{7}
$$

Thus we assume that the excluded volume effect is related to *n*(**r**) while the density itself is *n*(**r**) which may be considerably more finely structured. The free energy density functional is then

$$\mathcal{F}\_{\mathcal{C}}[n(\mathbf{r})] = k\_B T \int\_V d\mathbf{r} \, n(\mathbf{r}) \ln \left( \frac{n(\mathbf{r})}{1 - \overline{n}(\mathbf{r}) \sigma^3} \right) + \int\_V d\mathbf{r} \, n(\mathbf{r}) \left( \frac{1}{2} \int\_V d\mathbf{r}' \, n(\mathbf{r}') \phi\_{\mathrm{d}ll}(\mathbf{r}, \mathbf{r}') + V\_{\mathrm{ext}}(\mathbf{r}) \right) \, . \tag{8}$$

This functional can now describe hard sphere structuring at hard walls and related phenomena (Freasier et al. , 1989; Nordholm et al. , 1980).

Other equations of state, and coarse-graining approaches, have been suggested. For instance, Tarazona and Evans (Tarazona & Evans , 1984) based their free energy functional on the Carnahan-Starling (Carnahan & Starling , 1969) equation of state for hard spheres which yields a configurational free energy per particle

$$f\_{\mathcal{L}} = k\_B T \left[ \ln n\_b - 1 + \frac{\eta (4 - 3\eta)}{\left(1 - \eta\right)^3} \right] \tag{9}$$

in terms of the volume fraction *η* = *πσ*3*nb*/6. The "weighting" of the particle density was done by Tarazona (Tarazona , 1984) according to

$$
\overline{m}(\mathbf{r}) = \int d\mathbf{r}' \, w(|\mathbf{r} - \mathbf{r}'| \because n(\mathbf{r})) n(\mathbf{r}') \, \tag{10}
$$

$$w(r) = w\_0(r) + w\_1(r)n + w\_2(r)n^2 \quad , \tag{11}$$

where the weighting functions *wi* are determined to fit the direct correlation function and the radial distribution function for a uniform hard sphere fluid as described by the Percus-Yevick integral equation theory.

6 Will-be-set-by-IN-TECH

The key to many phenomena in nonuniform fluids is the nonlocality of interactions, that is the nonvanishing range of interaction on the length scale of density variations. From complete neglect of this nonlocality in the vdW analysis of phase separation we have proceeded above to account for the nonlocality of the pairwise attractive interaction leaving only the shorter range repulsive interactions to local approximation. There are, however, many applications where the nonlocality of the repulsive interactions is also important. We know, for example, that at higher densities hard spheres pack into distinct layers at hard walls. Such packing structures are not predicted by theories based on local representation of excluded volume effects. Moreover, it is clear that repulsive nonlocality will enter also into an accurate prediction of gas-liquid surface tension. For this reason fine-grained GvdW theory has been developed to allow both attractive and repulsive interactions to be treated in a nonlocal representation. The key step is to realize that the excluded volume effect is related to an average particle density over a spatial domain determined by the range of the repulsive interactions. If we, as above, treat the particles as interacting hard spheres of diameter *σ* then the simplest coarse-graining average reflecting the domain of the repulsive interactions

*n*(**r**) = 3

 *n*(**r**) 1 − *n*(**r**)*σ*3

*fc* = *kBT*

*n*(**r**) =

 *d***r** � *w* |**r** − **r** � |; *n*(**r**) *n*(**r** �

phenomena (Freasier et al. , 1989; Nordholm et al. , 1980).

4*πσ*3

 + *V*

Thus we assume that the excluded volume effect is related to *n*(**r**) while the density itself is *n*(**r**) which may be considerably more finely structured. The free energy density functional is

This functional can now describe hard sphere structuring at hard walls and related

Other equations of state, and coarse-graining approaches, have been suggested. For instance, Tarazona and Evans (Tarazona & Evans , 1984) based their free energy functional on the Carnahan-Starling (Carnahan & Starling , 1969) equation of state for hard spheres which yields

in terms of the volume fraction *η* = *πσ*3*nb*/6. The "weighting" of the particle density was

where the weighting functions *wi* are determined to fit the direct correlation function and the radial distribution function for a uniform hard sphere fluid as described by the Percus-Yevick

ln *nb* − 1 + *η*(4 − 3*η*)

(1 − *η*) 3

*w*(*r*) = *w*0(*r*) + *w*1(*r*)*n* + *w*2(*r*)*n*2 , (11)

*d***r** � *n*(**r** �

*d***r** *n*(**r**)

1 2 *V d***r** � *n*(**r** � )*φatt*(**r**,**r** �

) . (7)

, (9)

) , (10)

) + *V*ext(**r**)

 . (8)

**3.1.1 The fine-grained GvdW theory**

is

then

F*c*[*n*(**r**)] = *kBT*

 *V*

*d***r** *n*(**r**)ln

a configurational free energy per particle

integral equation theory.

done by Tarazona (Tarazona , 1984) according to

#### **4. Density functional theory (DFT) of simple ionic systems**

The simple fluids considered above are all characterized by short range interactions where the second virial coefficient *B*2(*T*) is finite. In the case of ionic systems the Coulomb energy of two point charges *qi* and *qj* separated by a distance *r* reads

$$\phi\_{ij}^{\mathbb{C}}(r) = \frac{q\_i q\_j}{4\pi\varepsilon\_0 \varepsilon r} \; \text{ \; \, \, \tag{12}$$

where *ε*0 is the vacuum permeability and *ε* is the relative dielectric constant. Hence, the range of interaction is infinite and *B*2(*T*) diverges. This means that interactions have a more profound effect on these systems and the correlations between the ions play a more fundamental role than in the simple fluids discussed above.

#### **4.1 The one-component plasma (OCP)**

The simplest IL is the one-component plasma (OCP) where one of two components is treated as an inactive neutralizing background charge while the other is an active ionic component. An example could be a molten metal where the metal ions are taken to move in a background electron gas which is in its ground state. For active ions of magnitude *q* and density *nb* the bulk charge density is *qnb* while the neutralizing background charge density is −*qnb*. In the absence of correlations, it follows that the system would behave as an ideal gas with zero potential energy on average. The active ions do repel each other, however, so that compact configurations are suppressed, the average potential energy is negative and the OCP is stable. The long range nature of the interaction means each ion is surrounded by a like-charge-depletion zone (a hole) with a shape that exactly neutralizes its charge *q*. Thus correlations play a fundamental role in the OCP. The GvdW free energy density functional for the ion density around a central ion can be written as

$$\mathcal{F}\_{\varepsilon}[n(\mathbf{r})] = k\_{\text{B}}T \int\_{V} d\mathbf{r} \, n(\mathbf{r}) \ln n(\mathbf{r}) + \frac{q^{2}}{4\pi\varepsilon\_{0}\varepsilon} \int\_{V} d\mathbf{r} \, \Delta n(\mathbf{r}) \left(\frac{1}{2} \int\_{V} d\mathbf{r}' \frac{\Delta n(\mathbf{r}')}{|\mathbf{r} - \mathbf{r}'|} + \frac{1}{|\mathbf{r}|}\right) \quad . \tag{13}$$

Here Δ*n*(**r**) = *n*(**r**) − *nb* is the deviation of the ion density from its bulk value. On introducing the Lagrange multiplier *α* to accommodate the canonical constraint of fixed particle number, the variational derivative of this functional yields

$$\frac{\delta \mathcal{F}\_{\mathcal{C}}}{\delta n} = k\_B T \{ \ln n(r) + 1 \} + a + q \psi(r) \tag{14}$$

$$
\psi(r) = \frac{q}{4\pi\varepsilon\_0\varepsilon} \left(\frac{1}{r} + \int\_V d\mathbf{r}' \frac{\Delta n(\mathbf{r}')}{|\mathbf{r} - \mathbf{r}'|}\right) \tag{15}
$$

where we recognise the mean-field electrostatic potential *ψ*(*r*). Setting this derivative to zero and writing *β* = 1/(*kBT*) yields the equation

$$n(r) = \exp\left(-1 - \beta \mathfrak{a} - \beta q \psi(r)\right) = n\_b \exp\left(-\beta q \psi(r)\right) \ . \tag{16}$$

This is the well known Poisson-Boltzmann equation: a mean-field approximation which can be solved by iteration. Linearizing this equation we obtain the Debye-Hückel equation

$$
\Delta n(r) = \frac{-\kappa^2 \exp(-\kappa r)}{4\pi r} \quad , \qquad \kappa^2 = \frac{n\_b q^2}{\varepsilon\_0 \varepsilon k\_B T} \; . \tag{17}
$$

Now, after defining the dimensionless coupling strength Γ by

$$
\Gamma = \frac{q^2}{4\pi\varepsilon\_0\varepsilon k\_B T a} \tag{18}
$$

where *a* = � 4*πnb*/3�−1/3 is the ion-sphere radius, we get the internal potential energy per particle

$$
\mu = -\frac{\sqrt{3}}{2} \, k\_B T \, \Gamma^{3/2} \, . \tag{19}
$$

The electrostatic contribution to the free energy is then

$$f\_{\rm el} = \int\_0^\Gamma \frac{d\Gamma'}{\Gamma'} \, u(\Gamma') = -\frac{1}{\sqrt{3}} \, k\_B T \, \Gamma^{3/2} \, . \tag{20}$$

The Debye-Hückel theory is both a mean-field and a linear response theory. The radial distribution function is

$$g\_{DH}(r) = 1 - \frac{\kappa^2}{4\pi n\_b} \frac{1}{r} \exp(-\kappa r) \ . \tag{21}$$

Since physically *g*(*r*) is greater or equal to zero the linear response cannot be applicable for small *r* where *g*(*r*) becomes negative. The Debye-Hückel hole (DHH) theory (Nordholm , 1984a) rephrases the Debye-Hückel theory to assign *g*(*r*) the value zero at these small *r* and then adds the exponential part continuously to give

$$g(r) = \begin{cases} 0 & r < h \\ 1 - \frac{h}{r} \exp\{-\kappa(r-h)\} & r \ge h \end{cases}.\tag{22}$$

Thus the *κ*-value is retained but the hole is reshaped to be physically realizable and consistent with the asymptotic DH screening of the central charge. At the density *nb* the "hole" should contain one ion. Thus we have

$$n n\_b \int\_0^\infty dr \, 4\pi r^2 \left(1 - g(r)\right) = 4\pi n\_b \left[\frac{h^3}{3} + \int\_h^\infty dr \, r h \exp\left(-\kappa (r - h)\right)\right] = 1 \quad , \tag{23}$$

so that

$$h = \kappa^{-1} \left[ \left( 1 + (3\Gamma)^{3/2} \right)^{1/3} - 1 \right] \,. \tag{24}$$

In the DHH theory the charge screening distribution can be integrated over the Coulomb energy (12) to yield

$$\mu = -\frac{k\_B T}{4} \left[ \left( 1 + (3\Gamma)^{3/2} \right)^{2/3} - 1 \right] \,\, . \tag{25}$$

As before, the electrostatic free energy can be obtained by a charging integration

$$f\_{\rm el} = \frac{k\_{\rm B}T}{4} \left( 1 + \frac{2\pi}{3\sqrt{3}} + \ln\left(\frac{\omega^2 + \omega + 1}{3}\right) - \omega^2 - \frac{2}{\sqrt{3}}\arctan\left(\frac{2\omega + 1}{\sqrt{3}}\right) \right) \tag{26}$$

where

$$
\omega = \left( 1 + (3\Gamma)^{3/2} \right)^{1/3} \,\, . \tag{27}
$$

The DHH theory (Nordholm , 1984a; Penfold et al. , 1991) produces thermodynamic properties in good agreement with simulation results for the full fluid range of coupling constants 1 *<* Γ *<* 160. We note that when the ions are imbedded in hard spheres, as in the RPM of electrolytes, the DHH theory above should be applicable when the hard sphere diameter *σ* is smaller than the electrostatic hole size *h* above. The theory has been applied in the study of polyelectrolyte screening (Penfold et al. , 1990) and we shall find use for the DHH model of coion correlations in our treatment of ILs below. While the DHH provides an analytic solution to the a charged fluid problem, the more general DFT's discussed later will require numerical iterative solutions. With this in mind, it is useful to consider the convergence properties of the DHH theory, where an analytical solution is available, in order to elucidate some of the issues that may confront the numerical approaches.

#### **4.1.1 An iterative scheme**

8 Will-be-set-by-IN-TECH

4*πnb*/3�−1/3 is the ion-sphere radius, we get the internal potential energy per

) = − 1 √3

1

0 , *r < h*

� ∞ *h*

*dr rh* exp�

− 1 �

> − 1 �

> > √3

arctan

− *ω*2 − 2

−*κ*(*r* − *h*)

1 + (3Γ)3/2�1/3

1 + (3Γ)3/2�2/3

�

1 + (3Γ)3/2�1/3

In the DHH theory the charge screening distribution can be integrated over the Coulomb

Thus the *κ*-value is retained but the hole is reshaped to be physically realizable and consistent with the asymptotic DH screening of the central charge. At the density *nb* the "hole" should

> � *h*3 3 +

The Debye-Hückel theory is both a mean-field and a linear response theory. The radial

4*πnb*

Since physically *g*(*r*) is greater or equal to zero the linear response cannot be applicable for small *r* where *g*(*r*) becomes negative. The Debye-Hückel hole (DHH) theory (Nordholm , 1984a) rephrases the Debye-Hückel theory to assign *g*(*r*) the value zero at these small *r* and

4*πε*0*εkBTa* , (18)

2 *kBT* Γ3/2 . (19)

*kBT* Γ3/2 . (20)

*r* exp(−*κr*) . (21)

� , *r* ≥ *h* . (22)

� �

. (24)

. (25)

��

, (26)

�2*ω* + 1 √3

. (27)

= 1 , (23)

−*κ*(*r* − *h*)

Γ = *q*2

√3

*u* = −

*d*Γ� Γ� *u*(Γ�

*gDH*(*r*) = 1 − *κ*2

Now, after defining the dimensionless coupling strength Γ by

The electrostatic contribution to the free energy is then

then adds the exponential part continuously to give

*g*(*r*) =

⎧ ⎨ ⎩

1 − *g*(*r*) �

*h* = *κ*−1

*u* = −*kBT* 4

1 − *h*

*r* exp�

= 4*πnb*

��

��

As before, the electrostatic free energy can be obtained by a charging integration

*ω* = �

� *ω*2 + *ω* + 1 3

*f*el =

� Γ 0

where *a* = �

distribution function is

contain one ion. Thus we have

*dr* 4*πr* 2�

> � 1 +

2*π* 3 √3 + ln

*nb* � ∞ 0

energy (12) to yield

*f*el = *kBT* 4

so that

where

particle

By considering a simple statistical mechanical theory of bulk point plasmas, we aim to elucidate some general convergence characteristics of successive approximation schemes that are often deployed in free energy density functional calculations for Coulomb fluids. Although only a very crude level of theory is discussed here, we suggest that the conclusions are relevant for more sophisticated analyses of ILs.

Although the explicit DHH solution for the OCP is very valuable, the corresponding generalisation for arbitrary plasmas presents an ill-posed inverse problem (Penfold et al. , 2005) since the number of parameters grows faster with mixture components than does the number of constraint equations. Nevertheless, a simple iterative scheme can be formulated that resembles the successive approximation algorithms often used for more sophisticated free energy density functional calculations of Coulomb fluids. To study the convergence properties of this scheme in detail, it is useful to re-cast the procedure in terms of the DHH theory for the OCP. Recall the essential idea here is to "renormalise" the Debye-Hückel radial distribution function

$$\mathcal{g}(r) = \begin{cases} 0 & r < h \\ 1 + \mathfrak{F}(h) \left( \mathcal{g}\_{DH}(r) - 1 \right) & r \ge h \end{cases} \tag{28}$$

where *gDH* is given by (21) and *ξ* = *ξ*(*h*) is an asymptotic charge scaling factor fixed to ensure that *g* remains continuous and non-negative everywhere. By virtue of the previous analysis the appropriate DHH value is *ξ*(∞) = *h* exp(*κh*)/(*a*Γ), but the following successive approximation scheme immediately suggests itself. Set the iteration counter *p* = 0 and specify an initial guess *ξ*(0). Now, invoke the continuity constraint to fix the current value for electrostatic exclusion hole size as the positive root of

$$M\exp(l) - \sqrt{3}\Gamma^{3/2}\,\xi^{(p)} = 0 \quad , \tag{29}$$

where *l* = *κh*. Next, appeal to the local electroneutrality condition (23) and obtain a new estimate for the renormalisation factor

$$
\xi^{(p+1)} = \frac{1 - (\Im \Gamma)^{-3/2} \, l^3}{(1+l) \exp(-l)} \, . \tag{30}
$$

Increment the counter *p* → *p* + 1 and repeat the nested cycle of equations (29) and (30) until convergence. This iterative procedure is related to the normal method of solving the nonlinear PB equation for screening charge densities. Indeed, the DHH analysis can be regarded as

Fig. 1. A representative family of iteration function *G* for generating the sequence of asymptotic charge renormalisation factors. Coupling constant Γ values corresponding to each curve are indicated and span the range consistent with an OCP in the fluid state (0 *<* Γ *<* 140). Lines of slope ±1 are indicated by the symbols.

a very simplified PB theory. It will be convenient to express the iteration sequence by the recurrence relation

$$\mathfrak{F}^{(p+1)} = \frac{\left( (\mathfrak{J}\Gamma)^{3/2} - l^3 \right)}{\mathfrak{J}l(1+l)} \mathfrak{F}^{(p)} \equiv \mathcal{G} \left( \mathfrak{F}^{(p)} \right) \; , \tag{31}$$

where *l* = *l ξ*(*p*) is given implicitly by (29) and the result (24) is clearly manifest as fixed points of the iteration function *G*.

Recall that, in general, the differentiable function *G* = *G*(*ξ*) is a contraction mapping on the compact interval [*a*, *b*] if:

i) *ξ* ∈ [*a*, *b*] ⇒ *G*(*ξ*) ∈ [*a*, *b*]

ii) there exists a (Lipschitz) constant *L <* 1 such that, |*G*� (*ξ*)| ≤ *L* for all *ξ* ∈ [*a*, *b*].

For the OCP, charge reversal phenomena can be ignored so that *ξ* must remain positive. Consequently, by virtue of (29) and (31), the renormalisation factor is bounded

$$0 \le \mathfrak{F}^{(p)} \le \mathfrak{f}\_{\text{max}} = \frac{\exp\left(\sqrt{3\Gamma}\right)}{\Gamma} \,. \tag{32}$$

From Figure 1, it is evident that *ξ* ∈ 0, *ξ*max implies *G*(*ξ*) ∈ 0, *ξ*max and we are encouraged to expect that *G* will be a contraction map on a suitable sub-interval that also satisfies the gradient criterion (ii). Differentiation yields

$$G'(\xi) = \frac{l}{3(1+l)^3} \left( (3\Gamma)^{3/2} - \left( (1+l)^3 - 1 \right) \right)$$

$$\geq \frac{l}{3(1+l)^3} \left( l^3 - \left( (1+l)^3 - 1 \right) \right) = -\left( \frac{l}{1+l} \right)^2 > -1 \,, \tag{33}$$

where *ξ >* 0 demands (3Γ)3/2 *> l* 3 and the inequality is confirmed by Figure 1. The condition *G*� (*ξ*) *<* 1 is satisfied if and only if