View Single Page Version

Ionic Liquid (EMIM-OAc)

The ionic liquid 1-ethyl-3-imidazolium acetate

60 minutes

Summary

This example focuses on setting up an ionic liquid (1-ethyl-3-imidazolium acetate) and reproducing some of the published results in “Structure and Dynamics of 1-Ethyl-3-methylimidazolium Acetate via Molecular Dynamics and Neutron Diffraction”, D. T. Bowron, C. D’Agostino, L. F. Gladden, C. Hardacre, J. D. Holbrey, M. C. Lagunas, J. McGregor, M. D. Mantle, C. L. Mullan, and T. G. A. Youngs, J. Phys. Chem. B 114, 7760 (2010). The authors used the SANDALS instrument to collect neutron data on seven isotopic samples at 323 K - the data as provided here have been reprocessed in full, and so may differ slightly from those presented in the original article.

Preparation

Download the example data as a zip or tar.gz and unpack it if you haven’t done so already. Alternatively, download the data files listed below and save them to a location of your choice.

File ⇨ New

File ⇨ Save As…

Save your input file under a sensible name in same directory as the example data

Step 1 - Create the Ionic Liquid Species

We need to create the cation and anion species for our simulation. Here we’ll just import them from some existing xyz coordinate files that we just happen to have laying around.

1-Ethyl-3-methylimidazolium Cation

First, the 1-ethyl-3-methylimidazolium cation, which we’ll call “EMIM” for short from here on in:

Species ⇨ Import ⇨ From XYZ…

Choose the emim.xyz file.

Rename the species to   EMIM   by double-clicking on the tab name

Once you’ve selected the xyz file Dissolve will allow you to perform any last minute edits to the structure should you need to. The principal reason might be to correct the bonding between atoms, which Dissolve automatically calculates when the xyz file is read in.

The structure in the file is fine as-is, so:

Exit the species editor by clicking OK

Once a species has been created you cannot add or remove atoms or bonds, as these now represent immutable data in your simulation. If you need to change something, you’ll have to create a new species.

The creation of the species is now complete, and should look like this:

The imported EMIM cation structure

For the uninitiated, the EMIM cation is a five-membered imidazole which has been alkylated on both nitrogen atoms (with an ethyl and a methyl group respectively), resulting in a +1 overall charge.

Acetate Anion

Second, the acetate anion, which we’ll refer to as “OAc”.

Species ⇨ Import ⇨ From XYZ…

Choose the oac.xyz file.

Again, the structure in the file is fine as-is, so click OK to accept the structure.

Rename the species to   Acetate   by double-clicking on the tab name

The imported OAc anion structure

And there you have it - acetic acid which has been deprotonated to give us an anion with -1 overall charge to balance out our cation.

Step 2 - Set up a Starting Forcefield

Of course we need a forcefield to describe our ionic liquid, and we’ll focus on the widely-used parameters of Padua and Canongia-Lopes. In Dissolve this forcefield is implemented in two distinct parts, one covering cations and one covering anions.

EMIM Cation

Go to the   EMIM   species tab

Species ⇨ Add Forcefield Terms…

From the available forcefields choose PCL2019/Cations and click Next

We will use the default Determine atom types for all atoms option to add atom types for every atom in the species, so click Next

There will be no conflicts with existing atom types as there are no atom types already defined, so click Next

For the intramolecular terms we want to assign them and reduce to master terms (the default settings) so click Next to proceed

There will be no conflicts with existing master terms, so click Finish to exit the wizard

Now we have a full description of our cation - if you look on the Forcefield section on the species tab you can see that Dissolve has calculated the total charge on the molecule from both atom types and species atoms. Both of these should be +1.

You’ll notice that there are warning icons next to the total species atom and atom type charges. These are there only to indicate that the total charge is non-zero - we know that our species are ions, so we ignore these warnings.

OAc Anion

Now for the anion:

Go to the   Acetate   species tab

Species ⇨ Add Forcefield Terms…

From the available forcefields choose PCL2019/Anions and click Next

We will use the default Determine atom types for all atoms option to add atom types for every atom in the species, so click Next

At this point you will see that there is a conflict between one of the atom types Dissolve wants to apply to the OAc species (HC) and those that we’ve just created on the cation. While from a forcefield perspective the atom type HC is transferable between the cation and the anion, when we’re thinking about refining pair potentials against experimental data it is probably wise to separate it into two types so that we can allow both to vary independently in the refinement process. So:

Select the HC entry in the atom types list

Click the Suffix and type in A. Click Ok to change the name from HC to HCA to resolve the naming conflict.

Click Next to continue

For the intramolecular terms we want to assign them and reduce to master terms (the default settings) so click Next to proceed

Click Finish to exit the wizard

Step 3 - Apply Charge Scaling

At this point we’ll use some prior knowledge of the EMIM-OAc system. Many ionic liquid forcefields to make the quite reasonable assumption of integer (i.e. +1 and –1) charges on the ions, however the polarisation of (or charge-transfer between) ionic liquid ions has been the subject of significant debate in the literature, and various simulation studies show that scaling the ion charges to mimic a “static” charge transfer between ions can yield improvements to calculated quantities such as the diffusion coefficient and ion-ion structure.

We’ll apply a charge scaling on the ions following the work in “Structure and dynamics of 1-ethyl-3-methylimidazolium acetate via molecular dynamics and neutron diffraction”. Therein, the authors scaled charges on the ions by a factor of 0.85 and showed improved agreement with experimental data (neutron, of course!):

Go to the   EMIM   species tab

Species ⇨ Scale Charges…

Change the scale factor to 0.85

You’ll see the scaling reflected in the charges within the table on the species tab’s Forcefield, and the total Atomic Charge there should now be +0.85. We need to do the same for the anion, of course:

Go to the   Acetate   species tab

Species ⇨ Scale Charges…

Apply the scale factor to 0.85

The final important step is to explicitly request that Dissolve use charges from the species atoms, as the default automatic choice is to use atom type charges wherever possible.

Go to the   Forcefield   tab, Pair Potentials section

In the Charge Source controls untick the Choose Automatically option and choose the Species Atoms option below

If we don’t do this then Dissolve will elect to use the charges on the atom types (which are still valid).

Step 4 - Make the Ionic Liquid

Time to create a configuration to represent the ionic liquid.

Generate the Configuration

Configuration ⇨ Create…

Choose both species and press Next

Leave the configuration type as Mixture and press Next

Leave the box style as Fixed Geometry, Undefined Size and press Next

Set the Density to 0.09985 and the units to atoms/A3, and the Multiplier to 100. This gives us a box with side length of 29.6398 Å which is quite small for an ionic liquid, but in the interests of time we’ll leave it as that. Press Finish to complete the wizard.

Double-click on the configuration tab’s title and change its name to   Bulk

We need to adjust our pair potential range which is currently set to 15 Å and exceeds the half-cell width:

Go to the   Forcefield   tab, Pair Potentials section

In the Pair Potentials controls, change the Range to 12.0

The final thing to do is change the temperature (our experimental data was measured at 323 K) and set a size factor for our box - since we have rings in our cation we need to apply a size factor during the initial equilibration to help prevent occurrences of interlocking rings.

Go to the   Bulk   configuration tab

On the configuration viewer toolbar, click the Temperature button

Set the Temperature to 323.0

Click the Show Avaliable Nodes button

Drag the SizeFactor node into the Node list

Select the SizeFactor node

Set the SizeFactor to 10.0

Initial box of 1-ethyl-3-methylimidazolium acetate

Equilibrate

Let’s now add a suitable evolution layer for the system and set things running:

Layer ⇨ Create ⇨ Evolution ⇨ Standard Molecular (MC/MD)

Now we can equilibrate our system:

Ctrl-F

Set the number of iterations to 1000 and hit OK

You might want to track the change in energy of the system in the Output tab of the Energy module, in the   Evolve (Standard)   tab.

After one thousand iterations the size factor of the configuration should have returned to 1.0 (check this on the configuration tab!) and the intermolecular energy of the system should have reached a stable negative value of about –100 kJ mol-1 (for 100 ion pairs). The charged nature of the species in the present system means that reaching equilibrium can take a lot longer than for a system containing only small, neutral molecules - bear this in mind when running such systems. Also note that the total energy of the configuration remains large and positive due to the intramolecular energy contributions.

Equilibrated box of 1-ethyl-3-methylimidazolium acetate

Step 5 - Structure Factors

Time to create some isotopologues for our system and calculate the structure factors.

Isotopologues

We need a partly-deuteriated isotopologue for the cation where the alkyl groups are fully deuteriated but the imidazole ring hydrogens remain as protons.

Go to the   EMIM   species tab, Isotopologues section

Select the Isotopologues section to create a new isotopologue definition

Click Add to create a new isotopologue definition

Right-click on the name of the new isotopologue and choose Set hydrogens to ⇨ Deuteriated from the context menu to set all hydrogen types to deuterium

Change the isotopes for the HCW and HCR atom types back to to Natural (bc = -3.739 fm)

Double-click the name of the isotopologue and change it to H3D8

For the anion we need a fully deuteriated version:

Go to the   Acetate   species tab, Isotopologues section

Click Add to create a new isotopologue definition

Change the isotope for the HCA atom type to 2 (bc = 6.671 fm)

Double-click the name of the isotopologue and change it to D3

Structure Factors

Now we can create a new layer to calculate the RDF and structure factors:

Layer ⇨ Create ⇨ Correlations ⇨ RDF and Neutron S(Q)

Our reference datasets are listed in the following table - we’ll refer to them in the form X-Y where X is the isotopologue of the cation and anion respectively:

Cation Anion ID Datafile
1 Natural Natural H-H SLS39546-H_H.mint01
2 Natural D3 H-D3 SLS39548-H_D.mint01
3 H3D8 Natural H3D8-H SLS39549-H3D8_H.mint01
4 H3D8 D3 H3D8-D3 SLS39550-H3D8_D.mint01
5 Natural H:D3 50:50 H-50 SLS39567-H_50.mint01
6 H:H3D8 50:50 H:D3 50:50 50-50 SLS39568-50_50.mint01

Note that there are no exchangeable hydrogen atoms in this system.

The existing NeutronSQ represents sample 1 in the table since it will default to use the natural isotopologues for the cation and anion, so all that remains to do for that one is to set its identifying name and the reference datafile:

Click on the first NeutronSQ module to select it

Double-click on its name and change it to H-H

Find the Reference Data settings group, and for the Reference keyword select the file SLS39546-H_H.mint01 and set the format of the data to mint

For the remaining datasets you need to drag a new NeutronSQ into the layer from the module palette (click the Show Avaliable Modules button if you can’t already see it) and set the relevant isotopologues, short name, and datafile as listed above.

Finally, we need to set the instrumental broadening in the SQ module first:

Click on the SQ module to select it

Set the QBroadening to OmegaDependentGaussian with a FWHM of 0.02

Step 6 - Initial Comparison

Now that we have our structure factor calculations set up let’s run for a few steps and see what they look like.

Ctrl-5

You can go to the Output tab of any of the NeutronSQ modules and see the comparison of the F(Q) - we re-plot all seven together for convenience below (this is output from the EPSR module):

Calculated total structure factors (dashed lines) compared to the experimental data (solid lines) in the equilibrated simulation

Overall the calculate F(Q) agree rather well with the experimental curves, with the exception that the peak at Q = 1.5 Å-1 is too intense and well-defined, particularly for the H3D8-D and H3D8-D samples. This peak is related to ion-ion neighbour interactions - we can lightly abuse an equation related to Bragg scattering, $r = 2 \pi / Q$, to calculate the real-space correlation distance as 4.19 Å - and so suggests that the forcefield, even with the scaled charges, is giving a liquid that is overstructured from too-strong interactions between the cation and anion.

We’ll let the EPSR module sort that out for us!

Layer ⇨ Create ⇨ Refinement ⇨ Standard EPSR

Because of the stronger intermolecular interactions manifesting in a system of ions we will need to push the refinement a bit harder than usual. You can run the refinement with the default value of EReq of 3, but you will see very little improvement in the structure factors. Standard practice is to then sequentially increase this value in small steps (e.g. increments of 3 at a time) until no further improvement to the fits (as measured by the total R-factor) can be achieved. Simply putting in a huge number that is “definitely enough” is not a good idea - experimental data is imperfect to some degree, and at some point the generated potentials are highly likely to force the system into a worse state than it was before.

For the present system a value of about 17 is optimal, and we’ll put this in straight away:

Go to the   Refine (EPSR)   layer tab

Select the EPSR module

In the module Options change the EReq value to 17.0

Click the Output tab

Running the simulation for another 1000 iterations should give us a good refined structure:

Ctrl-F

Set the number of iterations to 1000 and hit OK

Step 7 - Refined Data

With the 1000 iterations finished the structure factors should look something like this:

Refined total structure factors (dashed lines) compared to the experimental data (solid lines).

The total r-factor should be around a value of 4.8×10-4. This is a significant improvement from where we began with our reference forcefield, and the agreement for all datasets is now pretty good. We can of course be critical and say, for example, that the H3D8-H sample (green data in the above plot) has the worst agreement with the “doublet” ion-ion peak too pronounced in the simulation. However, what we can say is that the current simulation is overall a lot more consistent with the experimental data than the reference forcefield was, and so whatever properties we calculate from it are a better reflection of the measured reality.

Step 8 - Analysis Sites

Time to calculate some properties, but first we need to create some analysis sites on our ions.

Cation

To be consistent with the original study our focal site on the cation will be at the midpoint of the two nitrogens of the imidazolium:

Go to the   EMIM   species tab, and click on the Sites section

Select both nitrogen atoms by holding down the Shift key and clicking on them in the viewer (you can also draw individual boxes around them)

Right click the selection and click Create site from ⇨ Atoms to create a new site with the origin at the midpoint of the these two atoms

Now select one of the nitrogens on its own and right click it

Click Modify current site ⇨ Set X-Axis Atoms to define the direction of the x axis

Finally, select the carbon atom in between the nitrogens and right click it, selecting Modify current site ⇨ Set Y-Axis Atoms in the toolbar to define the direction of the y axis

Double-click the site in the list and rename it to RingNN

Anion

For the anion we’ll use a simple centre-of-mass site:

Go to the   Acetate   species tab, and click on the Sites section

Select all atoms in the acetate (draw a box around everything with the mouse)

Right click the selection and click Create site from ⇨ Atoms to create a new site

In the site Definition the Weight by mass option

This locates the site more-or-less on the carbon atom of the carboxylate group

Double-click the site in the list and rename it to COM

Step 9 - Ion-Ion RDFs

Here we’ll set up calculation of some basic ionic liquid properties in the form of the three unique ion-ion radial distribution functions - the cation-anion, cation-cation, and anion-anion flavours - and the cation-anion coordination number.

So, let’s create and prepare suitable layer:

Layer ⇨ Create ⇨ Analysis ⇨ RDF & Coordination Number

Go to the   Analyse RDF/CN   layer

Double-click the “SiteRDF” module in the list and rename it to RDF-CA

Change the DistanceRange maximum to 14.0 (slightly less than our half-cell length)

Set SiteA to RingNN on the cation

Set SiteB to COM on the anion

The SiteRDF module can also calculate simple coordination numbers over specific distance ranges, however typically we don’t know what this range is until we’ve seen the actual RDF data, so we’ll leave that calculation turned off for now.

Now to add on two more SiteRDF modules for the other correlations of interest:

Drag a SiteRDF module from the palette up to the module list

Double-click the new SiteRDF module and rename it to RDF-CC

Change the DistanceRange maximum to 14.0

Set SiteA and SiteB to RingNN on the cation

Drag another SiteRDF module from the palette up to the module list

Double-click the second new SiteRDF module and rename it to RDF-AA

Change the DistanceRange maximum to 14.0

Set SiteA and SiteB to COM on the anion

Now start up the simulation again:

Ctrl-F

Set the number of iterations to 250 and hit OK

Because our system is quite small (100 ion pairs) the RDFs will take a few thousand iterations to become properly averaged, but 250 steps should be enough to get a feel for the results. If you look at the output from the RDF-CA module you should see that the first minimum occurs at around 7 Å, so we’ll set this value as the upper limit for our coordination number calculation:

Find the Coordination Numbers options in the SiteRDF module

Set RangeA to span 0.0 to 7.0

Click RangeAEnabled to request the calculation

Start the simulation up once more:

Ctrl-R

The coordination number should average out at about 6.2 anions (within 7.0 Å of a central cation), compared to 6.8 (EPSR) or 6.9 (pure molecular dynamics) in the 2010 paper.

Exercises

Follow-on exercises:

    1. Partial g(r) - Set up and calculate partial g(r) between the cation ring hydrogens and anion oxygens

Exercise 1 - Partial g(r)

While the ion-ion radial distribution functions give us a general picture of the liquid structure, there are some specific interaction sites which drive a lot of the finer details. In the present ionic liquid those interactions are between the ring hydrogens of the imidazolium (which have been assigned the HCR and HCW atom types) and the acetate oxygen atoms. The asymmetry of the cation and the different chemistry of the ring hydrogens means that the partial g(r) for all three sites are going to be slightly different.

Tips:

  1. You’ll need to add individual analysis sites for the oxygen atoms of the acetate anion, and reference both in the analysis layers.