View Single Page Version
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.
- Cation Coordinates: emim.xyz
- Anion Coordinates: oac.xyz
- Fully Protiated IL Data:
SLS39546-H_H.mint01
- Cation-H / Anion-D3 Data:
SLS39548-H_D.mint01
- Cation-D8 / Anion-H Data:
SLS39549-H3D8_H.mint01
- Cation-D8 / Anion-D3 Data:
SLS39550-H3D8_D.mint01
- Cation-H / Anion-H:D3 Data:
SLS39567-H_50.mint01
- Cation-H:D8 / Anion-H:D3 Data:
SLS39568-50_50.mint01
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
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
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
We will use the default Determine atom types for all atoms option to add atom types for every atom in the species, so click
There will be no conflicts with existing atom types as there are no atom types already defined, so click
For the intramolecular terms we want to assign them and reduce to master terms (the default settings) so click
to proceed
There will be no conflicts with existing master terms, so click
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
We will use the default Determine atom types for all atoms option to add atom types for every atom in the species, so click
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
and type inA
. Click to change the name fromHC
toHCA
to resolve the naming conflict.
Click
to continue
For the intramolecular terms we want to assign them and reduce to master terms (the default settings) so click
to proceed
Click
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, 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
Leave the configuration type as Mixture and press
Leave the box style as Fixed Geometry, Undefined Size and press
Set the Density to
0.09985
and the units toatoms/A3
, and the Multiplier to100
. 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 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, 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
button
Set the Temperature to
323.0
Click the
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
You might want to track the change in energy of the system in the Energy module, in the Evolve (Standard) tab.
tab of theAfter 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
section to create a new isotopologue definition
Click
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
andHCR
atom types back to toNatural (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, section
Click
to create a new isotopologue definition
Change the isotope for the
HCA
atom type to2 (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 tomint
For the remaining datasets you need to drag a new NeutronSQ into the layer from the module palette (click the 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 of0.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 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):
tab of any of the
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
change the EReq value to17.0
Click the
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
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 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
⇨ 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
⇨ to define the direction of the x axis
Finally, select the carbon atom in between the nitrogens and right click it, selecting
⇨ 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 section
Select all atoms in the acetate (draw a box around everything with the mouse)
Right click the selection and click
⇨ 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
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
to7.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:
-
- 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:
- You’ll need to add individual analysis sites for the oxygen atoms of the acetate anion, and reference both in the analysis layers.