View Single Page Version

A Simple Aromatic

Structure of liquid benzene from isotopic neutron data

80 minutes

Summary

You might always be asking the question “How much does quantity X differ between the equilibrated simulation at the reference forcefield parameters, and the refined simulation at the empirical parameters?”. This is a very important question, and this example will show how to checkpoint simulation data for comparison later on. It will also introduce the utility of size factor scaling when dealing with molecular systems containing rings.

You will set up and equilibrate a small liquid benzene simulation before calculating a few properties of interest, and then refine the potential against experimental data and recalculate the properties of interest to see how much they change. The data are three neutron scattering datasets measured on the NIMROD instrument at the ISIS Pulsed Neutron and Muon Source in 2014. We’ll be referencing the paper by Headen et al. along the way - “Structure of π-π Interactions in Aromatic Liquids”, T. F. Headen, C. A. Howard, N. T. Skipper, M. A. Wilkinson, D. T. Bowron and A. K. Soper, J. Am. Chem. Soc. 132, 5735 (2010).

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 Species

Draw the Molecule

The first thing we’ll do is draw a benzene molecule as best we can:

Species ⇨ Create ⇨ Draw…

A new editor window opens in which we can draw the rough geometry and connectivity for our molecule.

Select draw mode from the toolbar above the species viewer

Draw a roughly hexagonal ring of six carbon atoms (carbon is the default drawing element) by left-click-dragging in the viewer

Change the drawing element from C to H by clicking on the button next to the draw mode icon

Connect a single hydrogen atom to each carbon by left-click-dragging from each carbon site

You can create bonds by left-click-dragging between two existing atoms.

Your final species will look something like this:

A very badly drawn benzene molecule

Click OK to close the editor and create the new species

Double-click on the   NewSpecies   tab and change its name to   Benzene

Apply a Forcefield

Time to make it a little prettier! We’ll assign a standard forcefield to it, and optimise the geometry:

Species ⇨ Add Forcefield Terms…

From the available forcefields choose OPLSAA2005/Aromatics and click Next

You can filter forcefields by keywords in name and description by using the filter box at the top-right of the forcefield selection controls.

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

Click the Medic icon in the species viewer toolbar to optimise the geometry of the molecule using the forcefield we’ve just applied

We will also get ahead here and edit the master terms to reflect the geometry observed in the experimental data, since the forcefield we’ve applied here doesn’t get things quite right.

Go to the   Forcefield   tab, Master Terms section

In the Bonds table change the equilibrium bond length parameter (eq) of the ‘CA-HA’ bond term from 1.08 to 1.09 Å

Change the equilibrium bond length of the ‘CA-CA’ bond term from 1.40 to 1.38 Å

Create Deuterated Isotopologue

Since some of the experimental data was measured on deuterated benzene, we’ll need to create a suitable C6D6 isotopologue:

Go to the   Benzene   species tab, and click on the Isotopologues section

Click Add to create a new isotopologue for the species

Change the isotope for the HA atom type from Natural (bc = -3.739 fm) to 2 (bc = 6.671 fm) by clicking on the isotopologue entry and choosing from the drop-down menu

Change the name of the isotopologue to ‘C6D6’ (double-click on its name to do so)

The Natural isotopologue is always available and doesn’t need to be created, but doesn’t appear in the list of user-defined isotopologues. It is also not necessary to create “mixed” isotopologues (e.g. for 50:50 mixtures of H:D) as these are created by mixing individual isotopologues.

Add Analysis Sites

We’ll locate our analysis site at the centre of the benzene ring and give it some axes so that we may calculate orientational / spatial functions around it. The figure below shows the atoms we’ll select to define the origin (purple), x-axis (red) and y-axis (blue). Using these atoms as reference points for our coordinate system will set the XY plane to that of the ring, with the z axis perpendicular to the ring, pointing out of its center.

Origin (purple), x-axis (red) and y-axis (blue) atoms defining the oriented benzene site

Go to the   Benzene   species tab

Select all six carbon atoms by shift-clicking on them in the viewer (you may need to click reset view first, to see the whole molecule)

Right-click the selected atoms and click Create site from ⇨ Atoms This creates a new site with the Carbon atoms as the Origin Atoms

Now select the single carbon atom (as depicted in the image above) and right click it.

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

Finally, shift-click the pair of adjacent carbon atoms and right click on them

Click Modify current site ⇨ Set Y-Axis Atoms to define the direction of the y axis

Rename the site to COG by double-clicking its name

Step 2 - Set up a Configuration

Time to create a suitable liquid configuration for our simulation using the predefined “simple random mix” generator:

Configuration ⇨ Create…

Choose the benzene 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.876 and the units to g/cm3, and the Multiplier to 100. You’ll see that this gives us a box with side length of just over 24.5 Å which is fairly small, but in the interests of speed we’ll go with it for now so hit Finish .

Double-click on the   NewConfiguration   tab’s title and change its name to   Liquid

Initial, randomised benzene box containing 100 molecules

Since the molecules were added to the box completely at random, there will most likely be some that overlap quite a lot. Benzene has a ring which makes it possible for molecules to become permanently “interlocked”, but we can do our best to remove such contacts by using what Dissolve (and EPSR) call “size factor scaling”. Basically, this allows the configuration box and the positions of the molecules to be scaled (but keeping the internal geometry intact), removing bad overlaps and letting molecules move safely apart. The factor by which the box is scaled is reduced gradually until the initial requested system size is regained (i.e. corresponding to a size factor of 1.0).

In anticipation of some bad contacts:

Go to the   Liquid   configuration tab

Click the Show Available Nodes on the bottom-left of the window

Drag the SizeFactor node into the Procedure section to add it

Set the value for the Size Factor to 10.0

Finally, before we move on, we will have to reduce our pair potential range from the default of 15 Å since this is larger than the half the maximum width of our simulation box.

Go to the   Forcefield   tab, Pair Potentials section

In Control, reduce the pair potential Range from 15 to 12 Å

For a cubic box, you can’t have a cutoff value greater than half the box length, as this means that a given atom will ‘see’ other atoms twice because of the periodic boundary conditions.

Step 3 - Set up Equilibration

To equilibrate our system we’ll now add a standard Monte Carlo / Molecular Dynamics processing layer:

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

We can now start running our simulation and start to equilibrate the configuration box.

Ctrl-R

If you’re watching the box in the configuration tab you’ll notice that it intially “explodes” because of our defined size factor of 10.0, but the box will quickly reduce in size back to normal and result in a condensed phase with no interlocked benzene rings.

When the size factor has reduced the inter-molecular energy of the simulation will fluctuate a little, but should oscillate around some equilibrium value, at which point you can proceed. You can check this from the graphs of the Energy module:

Go to the   Evolve (Standard)   layer tab

Select the Energy module and go to its Output tab

Once this has been achieved, we can make a snapshot of the current simulation and store this equilibrated point as a backup. The restart file written by Dissolve as it runs stores all the information necessary to continue a stopped simulation, but there is no going back to a previous point in the simulation (e.g. to undo a bad analysis choice, or to reference earlier data) unless we have a suitable restart point to go back to. At any time we can write a new restart file from the GUI, independent of the standard restart file, and which we can keep to load in at a later date.

Esc

File ⇨ Save Restart Point…

Change the filename to equilibrated.restart and click Save

Now we can move on to set up the analysis of the structural properties we’re interested in.

Step 4 - Set up Analysis

For our “properties of interest” we’ll calculate both the centre-of-geometry RDF and the out-of-plane angle between molecules as a function of distance. We can get both from the AxisAngle module:

Module Purpose
AxisAngle Calculates the distance vs. axis angle map between two sites

Since our site (COG) is set up for the Z axis to be pointing out of the plane of the benzene ring, it is the angle formed between these axes on the different molecules that we’ll focus on.

Let’s add an empty processing layer to to the simulation, add on our module, and set it up:

Layer ⇨ Create ⇨ Empty

Double-click the layer tab and rename it to Analysis

Show the module palette for the current layer by clicking the Show Available Modules button at the bottom of the module list on the left

Drag a AxisAngle module up to the module list

Select the new AxisAngle module to show its Options

For the DistanceRange change the Bin Width to 10

For the AngleRange change the Max to 90

Enable the Symmetric option

This will account for the planar symmetry of the benzene ring, and map any calculated angles in the range 90 - 180° back in to 0 - 90°

Now we must tell the module which sites to use for its calculation:

Set both SiteA and SiteB to COG

Set the AxisA and AxisB to ZAxis

Enable the ExcludeSameMolecule to prevent the unwanted self-correlation spike at $r$ = 0 in the resulting RDF

It would also be nice here to calculate the spatial density function (the three-dimensional distribution of molecules around a reference point) and compare it, but this takes more iterations than a (sane) example will allow.

Step 5 - Calculate Baseline Properties

Now that the analysis processing layer is prepared, we need to run it for a while in order to get some good statistics on our quantities. A good 500 iterations (preferably 1000) is needed to get the statistics for the 2D map, so let’s ask Dissolve to run for a set number of steps:

Simulation ⇨ Run For…

Enter 1000 and click OK

Dissolve will now iterate for 1000 iterations and stop - you can see the estimated time remaining to complete these 1000 steps in the status bar at the bottom of the main window. Now might be a good time to have a drink, or write that email you’ve been putting off. If you think 1000 iterations is going to take too long, press Escape to stop the simulation and go for 500 iterations instead.

To see what’s being calculated, go back to the AxisAngle module and look at its output section:

Go to the   Analysis   layer

Click the AxisAngle module and select the Output section

We have three graphs - top-left is the radial distribution function between the centres-of-geometry of the benzene, while the bottom-left is the angle histogram of the z-axis angles (averaged over all distances). The main graph to the right shows the distance-angle map of the two quantities. You might want to explore the latter in 3D space, rather than a top-down view - you can change the view style for the plot in the toolbar at the top (change the view type from XY to 3D).

Once the iterations have completed, we need to do a bit of housekeeping.

Snapshot the Current Data

Let’s save the current data now, as it represents our reference point against which we’ll compare things later on:

File ⇨ Save Restart Point…

Change the filename to reference.restart and click Save

You might also find it useful to copy images of the graphs here to help with comparison later on - right-click on an empty part of the graph and select Copy to clipboard - or you can right-click on a particular curve and save the data directly in order to plot it in some other external software.

Clear the Current Data

Since we’re about to refine our simulation against experimental data and re-run the analysis, we don’t want to retain any of the currently-calculated quantities. We can request that Dissolve clear all of the data generated by modules and reset the simulation back to an “empty” state. Critically, the current contents of any configurations will remain untouched (i.e. our equilibrated box will remain in its present state, so we can proceed directly on with the simulation). It will also reset the iteration counter to zero.

Time to get rid of that data…

Simulation ⇨ Clear Module Data

Accept the consequences, and click OK

Turn off the Analysis Layer

We need to refine the simulation before we calculate our properties again, so we need to temporarily turn off our analysis layer:

Go to the   Analysis   layer

Disable the layer by clicking the Layer Enabled? button at the top-left of the tab so it displays

Step 6 - Set up Refinement

Time to refine. We’ll need two processing layers - one to calculate the g(r) and S(Q) from the simulation, and one to refine it against experimental data.

As noted in the introduction we have neutron-weighted experimental data, so we need a layer to calculate g(r) and neutron-weighted S(Q) for our three experimental datasets:

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

Open the module palette for the layer by clicking the Show Avaliable Modules button at the bottom left of the tab

Drag two additional NeutronSQ modules from the Correlation Functions category to the module list, placing them after the existing NeutronSQ module

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

Now let’s set up the three NeutronSQ modules to calculate what we want:

C6H6

Click on the first NeutronSQ module to select it

Double-click on its name and change it to C6H6

In the Reference Data settings group find the Reference keyword, select the file C6H6.mint01 and set the format of the data to mint

C6D6

Click on the second NeutronSQ module to select it

Change its name to C6D6

In the Control section click the button for the Isotopologue option

Press the Species button to add the natural isotopologue for the first species present

Change the isotopologue from Natural to C6D6

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

50:50 Mix

Click on the third NeutronSQ module to display its options

Change its name to 5050

In the Control section click the button for the Isotopologue option

Press the Species button to add the natural isotopologue for the first species present

Select either row in the table and click the Isotopologue button to insert the next unused isotopologue (i.e. the deuterated analogue)

The isotopologue mixture you’ve just defined has equal parts (weight = 1.0) of the Natural and C6D6 benzene isotopologues. You can use whatever ratios make sense - Dissolve will always normalise the total weight back to 1.0 internally.

Finally, find the Reference Data settings group and for the Reference keyword select the file 5050.mint01 and set the format of the data to mint

All that remains is to add our EPSR refinement layer:

Layer ⇨ Create ⇨ Refinement ⇨ Standard EPSR

All of the default values set there will be sufficient.

Step 7 - Refine the Simulation

Restart the simulation and monitor the progress of the r-factors and eReq value in the EPSR module.

Go to the   Refine (EPSR)   layer tab

Select the EPSR module and show its Output tab

Click the R Factors tab to display the errors

You’ll see here a graph of the individual R-factors for the three datasets, as well as the total (summed) R-factor. The current “strength” of applied potential is shown at the bottom-left of the plot (Current EReq). Once you’ve been running for a few hundred steps the ereq value should have increased to 3.0, and the fits might look something like this (on the F(Q) tab of the EPSR module):

Comparison of reference vs simulated structure factors at ereq=3.0

The fits should be really quite good, but you might observe some “ripples” in the simulated structure factors at lower $Q$ values, as is the case in the image above. These do not arise from structure in the system, but rather are high-frequency Fourier ripples caused by truncation of the partial radial distribution functions, made worse in the present case because our system isn’t that large (100 molecules). We can remove some of that effect by introducing a window function into the Fourier transform when generating our structure factors:

First, stop the simulation if it’s running:

Esc

Go to the   RDF / Neutron S(Q)   layer tab

Select the SQ module

In the Control settings group change the WindowFunction to Lorch0

Start the simulation running again and after a few steps you should see much smoother structure factors:

Comparison of reference vs simulated structure factors at ereq=3.0, and with the Lorch0 window function applied during Fourier transform of g(r) to S(Q)

Applying a window function in this way can have unwanted effects as well, with structural features “smoothed out”. Look closely at the C6D6 data around 1.5 < $Q$ < 2.0 and you can see that the pair of peaks has been somewhat blended together with the Lorch0 function applied.

With our refinement progressed to an acceptable level (in the data shown above the total r-factor is of the order of 5.5×10-4), we’re now in a position to save another restart point:

File ⇨ Save Restart Point…

Change the filename to ereq3.restart and click Save

Step 8 - Calculate Refined Properties

We now have our refined simulation that is in better agreement with the experimental data, so its time to recalculate our properties so that we can compare them to those we saved earlier. All the processing layers should now be enabled and remain on for the production run, including the refinement layer.

Go to the   Analysis   layer tab

Enable the layer by clicking the button in the Layer Control section so it displays

Simulation ⇨ Run For…

Enter 1000 and click OK

If you’re ready for another drink, now would be another opportune moment, or perhaps there’s another pressing email that needs hammering out?

Once those iterations have completed, we can move on to do the comparison of the data. Also, we can snapshot our production run data:

File ⇨ Save Restart Point…

Change the filename to production.restart and click Save

Step 9 - Compare Data

In the next two sections we’ll compare both the centre-of-geometry RDF and the axis angle map between the reference and refined simulations. How you do this is really up to you - you can load any of your restart points into Dissolve to see the data contained within, and then either:

  1. Export specific datasets to files and re-plot them side-by-side in your graphing software of choice.
  2. Take image snapshots / copy graphs to the clipboard and paste them in a document elsewhere.
  3. Run two copies of Dissolve, and load different restart points into each, and compare “interactively”.

In the following sections we’ve exported the basic 1D curves and plotted them in external software, and we compare the 2D data side-by-side. Your results won’t look exactly the same, but should show the same trends.

Centre-of-mass Radial Distribution Functions

We’ll begin with the centre-of-mass (ring centre) RDFs.

Equilibrated and refined centre-of-mass RDFs

As you can see the differences are quite subtle, but they are there. The refined (red) simulation shows a decrease in the main peak intensity, and an increase in probability of finding molecules at the shortest contact distances (around 4 Å). Looking at the probabilities between z-axis angles we observe that “flatter” angles are more likely in the refined simulation:

Equilibrated and refined Z-axis angle distributions

These functions only give us a hint at what the structural changes in the liquid are, however, and we need to combine them with the distance-angle map.

Distance/Angle Correlations

Equilibrated and refined distance/angle maps

The 2D function calculated by AxisAngle shows us variations in the centre-of-mass RDF with z-axis angle - the key feature is at the lowest angles (tending towards zero, at the bottom of the picture). The refined data on the right shows the emergence of a clear shoulder at 4 Å - as discussed in the paper by Headen et al. this corresponds to the presence of offset, parallel molecule stacking (see Figure 9 in that paper) which is minimal in the reference simulation. While it might be argued that this is a rather subtle change, it’s important to bear in mind that a typical off-the-shelf forcefield does not show this feature, and it only appears in a meaningful way when we refine against real-world, experimental data.