We have been using atomistic-level simulations to study ionic liquids for over ten years. Since then we have helped advance our understanding of these fascinating fluids by making predictions of thermodynamic and transport properties and providing insight into the link between the structure and properties of these liquids.
Volumetric Properties and Gas Solubilities
In very early work, we showed that volumetric properties such as density and isobaric expansivity could be computed with a high level of accuracy using simple intermolecular potential functions of the kind used for biological simulations (see Morrow and Maginn, and Shah et al.). In later work, we computed the solubility of various gases in ionic liquids, and explained the origin of the high solubility of CO2 in these liquids, a feature that has subsequently been exploited for use in gas separations and CO2 capture.
Using advanced Monte Carlo methods, we have computed full isotherms of CO2 in an ionic liquid (see an example on the left). Agreement with experimental measurements (open symbols) is remarkably good. We have applied these and related techniques to many other gases and ionic liquids.
In developing new ionic liquids for CO2 capture, it was proposed to add reactive amine groups to the cations to enable CO2 to be removed from flue gas at low partial pressures. While this does remove CO2, the resulting reacted fluids become extremely viscous. By carrying out molecular dynamics simulations of the reacted and unreacted species, we discovered that the large increase in viscosity was caused by the formation of salt bridges (like the "B-C" salt bridge shown on the right) between the reacted carbamate and ammonium groups on the cation. Details of this study by Gutowski and Maginn can be found here.
Aprotic Heterocyclic Anions for CO2 Capture
We have also carried out simulations to study the properties of a new class of ionic liquids based on "aprotic heterocyclic anions" (AHAs). These anions react with CO2 but do not form salt bridges, thereby elimminating the problem of viscosity increase. A recent paper on this work was published by Wu et al.
Vapor Pressure, Boiling Point and Critical Point
In some of our early work in 2002, we computed solubility parameters and the internal energy of vaporization. It was long thought that ionic liquids were "non-volatile" liquids, but these early simulations showed that while the enthalpy of vaporization is high, it is not infinite. Working with Steve Leone's group at Berkeley, we showed that the most likely volatile species from an ionic liquid at low temperature and pressure is a single ion pair (see image below).
We also computed the enthalpy of vaporization for a series of ionic liquids, and showed that the trends agree with experimental data, although there is a lot of variation.
If ionic liquids are volatile, what is their vapor pressure? Boiling point? Critical point? Do these latter quantities even exist? We used advanced Gibbs ensemble Monte Carlo simulations with Cassandra to address these questions (the graphic below depicts this type of simulation). We find that the boiling points and critical points are "theoretical" in that the estimates from our simulations are above the decomposition temperature. Nevertheless, the simulations give insight into the critical behavior and phase equilibria of these liquids, and it is different from conventional liquids like alkanes. We find that as the size of the cation increases, the boiling point and critical temperature decrease (due to a competition between van der Waals and Coulombic forces). At the highest temperaratures and pressures, significant clustering is observed, while at very low pressures, single ion pairs dominate. This work was described in two publications: a Faraday Discussions article (for TF2N ionic liquids) and a J. Phys. Chem. Lett. paper (for BF4 ionic liquids).
Melting point is one of the most fundamental and important properties of an ionic liquid. in fact, the common "definition" of an ionic liquid (a pure salt that melts below 100 oC) involves the melting point. Predicting the melting point of an ionic liquid using molecular simulation is a challenging task. If you simply run a constant pressure-constant temperature simulation and wait for the crystal to melt, you will have to wait a long time! That's because the free energy barrier for homogeneous nucleation of a liquid is large and the time- and lengthscales of a molecular dynamics calculation are short. Click on the movie below to see a simulation of a 1-n-butyl-3-methylimidazolium chloride crystal at a temperature well above the melting point. You will note that it doesn't melt! We showed (J. Chem. Phys. 136, 144116 (2012)) that even though standard interfacial melting point calculation methods work well for simple systems like the Lennard-Jones fluid, they fail for systems with sluggish dynamics (like ionic liquids). The best approach for computing melting points is to use free energy-based methods, where one can be confident of locating the thermodynamic melting point.
We have developed a procedure which we refer to as the "pseudo-supercritcal path sampling" or PSCP method for computing melting points. In this procedure, we wish to find the temperature at which the Gibbs energy of the liquid and crystal are equal. To do this we follow a two-step procedure. First, we compute the free energy of the liquid and crystal phases within a constant using Gibbs-Helmholtz integration:
On the left is the -H/RT2 at various temperatures for a model of argon (see: J. Chem. Phys., 2005, 122, 014115; J. Chem. Phys., 2006, 124, 164503; J.Chem. Phys., 2007, 127, 214504; I&EC Res., 2010, 49, 559-571; J. Chem. Phys., 2012, 136, 144116). On the right are the integrals referenced to an arbitrary reference temerature of 65K.
Next, we need to find the absolute difference in free energy between the crystal and liquid at one temperature. Once we have this, we can shift the curves on the right by this absolute difference and locate the point at which Gsolid = Gliquid. We do this by constructing a reversible path between the solid and liquid phases. In reality, the path between solid and liquid has a first order phase transition and is not integrable by standard free energy methods. We use a fictitious path that is continuous, much like the continuous path between liquid and vapor that passes above the critical point (hence the name of the method). This procedure is shown schematically below.
We will describe the method by following the red arrows from the crystal phase to the liquid phase. The interaction energy of the atoms in the crystal are gradually weakened while a tethering potential is turned on to keep the atoms in their lattice sites. This results in a "weak crystal" and the free energy along this path can be computed using standard methods (we use thermodynamic integration). Next, the density of the weak crystal is adjusted to the target liquid density; the free energy of this step is easily computed. Then the tethering potentials are turned off, resulting in a "weak liquid". We also use TI along this path to compute the free energy change for this process. Finally, the full potential of the liquid is restored gradually and the free energy for this is computed. By summing up the free energy differences along each step, the total free energy difference between liquid and crystal is determined at one temperature. This gives the melting point.
The PSCP method has been applied to a number of ionic liquids and found to give reliable results. See, for example, Phys. Chem. Chem. Phys., 2012,14, 12157-12164. The table below shows that the method is capable of predicting how the effect subtle changes in cation structure (such as adding a method group to the C2 position of the imidazolium ring) have on melting point can be captured with simulations.
Results from Zhang and Maginn, PCCP, 2012, 14, 12157.