View Single Page Version

Liquid Water

Structure of liquid water at 298 K, refining against three isotopically-labelled neutron datasets

60 minutes

Summary

It should be no surprise that water is present in a significant fraction of disordered materials studies, owing to its importance to just about everything, and its presence just about everywhere! Here we’ll set up, run, refine, and analyse the structure of liquid water at 298 K using three neutron scattering datasets measured on the SANDALS diffractometer at the ISIS Pulsed Neutron and Muon Source in 2001. This will illustrate most of the core workflows in Dissolve when applying it to the analysis of neutron scattering data, and is an excellent place to start if you’re new to the analysis of disordered materials.

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 a Water Species

Our first task is to create our water molecule species and give it a suitable forcefield description.

There are several ways to get atomic coordinates in to Dissolve, but since water is a simple molecule in this example we’ll draw the molecule in the GUI itself.

We’ll then set up the isotopologues we need to enable us to calculate the same weighted structure factors as those measured experimentally on the instrument.

Step 1a - Draw the Molecule

We will draw our molecule using Dissolve’s built-in editor:

Species ⇨ Create ⇨ Draw…

A new editor window opens in which we can create our new species.

Choose the drawing mode from the toolbar above the species viewer

Change the drawing element to oxygen - the current drawing element is represented on the button next to the drawing mode tool, defaulting to C , so click the button to bring up a periodic table for you to select a new element from

Left-click somewhere in the viewer below to place an oxygen atom

Change the drawing element to hydrogen

Left-click-drag from the existing oxygen atom to draw a hydrogen bound to it

Repeat for the other hydrogen

Click OK to close the editor and create the new species

Finally, let’s rename it:

Double-click on the new species tab’s title and change the name from NewSpecies to Water

A very badly drawn water molecule

You’ll now have a water molecule that might look like its gone ten rounds with Mike Tyson, so let’s clean it up. First, we’ll need to assign suitable forcefield terms. We’ll use a water model known as ‘SPC/Fw’ (Yujie Wu, Harald L. Tepper and Gregory A. Voth, ‘Flexible simple point-charge water model with improved liquid-state properties’, Journal of Chemical Physics 124 024503 (2006), http://dx.doi.org/10.1063/1.2136877) which basically means “simple point charge, flexible water”.

Species ⇨ Add Forcefield Terms…

From the available forcefields choose SPC/Fw and click Next

We need to assign new atom types for all of our water molecule’s atoms - the default option, Determine atom types for all atoms, does just that 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, which again are the default settings, so hit Next again

Dissolve will check to see if there are naming conflicts between new and existing master terms, just as it did for the atom types. There will be none as, again, we had no master terms to start with, so click Finish to exit the wizard

Take a look at the   Forcefield   section for the species and you’ll see that we now have atom types assigned to our atoms, and suitable intramolecular terms assigned to the bonds and angle. Note that the functional forms of the interactions are actually the names of master terms, and are preceded with @ to distinguish them as such (e.g. @HW-OW-HW). Master terms are global and can be referenced by one or more species, and are particularly useful when molecules possess high symmetry as there is no need to repeat the same parameter definitions. Furthermore, as we shall see later, adjusting species geometry by modifying the master terms is much easier than modifying all the individual values within a species.

We can now clean up the geometry of our molecule by doing a quick geometry optimisation.

Click the medic icon in the species viewer toolbar to optimise the geometry of the molecule

Step 1b - Add Isotopologues

We have three experimental datasets in this example - H2O, D2O, and a 50:50 mix of the two. The natural isotopologue (H2O) is defined automatically by Dissolve, so we don’t need to add it by hand. We will need to add a new one for D2O, but that is all. Mixtures (e.g. the H2O:D2O 50:50 sample) are created by blending isopologues, rather than defining specific isotopologues for the desired isotope ratios, and will be done when we set up the NeutronSQ modules.

Go to the   Water   species tab, and click the   Isotopologues   section

Click Add

This will create a new isotopologue containing every atom type in the species, defaulting to the natural isotope mix

Change the isotope for the HW atom type from Natural (bc = -3.739 fm) to 2 (bc = 6.671 fm)"

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

Users familiar with EPSR might expect to also set any exchangeable atoms here in the isotopologue definitions, but in Dissolve this is handled by specific settings in individual modules.

Step 2 - Set up a Configuration

Time to create a suitable bulk water configuration for our simulation.

Configuration ⇨ Create…

Choose the water 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

The Density we need is 0.1 atoms/A3 which just happens to be the default, and the standard Multiplier of 1000 is good for us, so press Finish to complete the wizard.

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

One thousand water molecules should have a cubic box size of just over 31 Å side length. You will see in the Add node that its Species option is not set to anything, which seems counter-intuitive. Specifying a species here instructs the Add node that $N$ identical copies of the species should be added to the box. In most cases it is preferable to introduce some variety to the geometry of the molecules, and so the Add node can instead reference a CoordinateSets node as the source, and which contains sets of coordinates representing randomised, evolved, or read in (from a file) coordinates for the species.

Before we proceed we’ll tweak the temperature of the configuration.

On the configuration viewer toolbar, click the Temperature button

Set the Temperature to 298 K

Initial, randomised water box containing 1000 molecules

Step 3 - Set up Basic Processing

We’ll now create two processing layers to handle the evolution of the configuration and the calculation of radial distribution functions and structure factors.

The Evolution Layer

We’ll use the standard molecule evolution layer for our system:

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

The new layer contains the following modules:

Module Purpose
MolShake Performs standard Monte Carlo moves on individual molecules
MD Performs a number of molecular dynamics steps, evolving the system according to Newton’s equations of motion
Energy Calculates the total energy of configurations

Monte Carlo simulation via the MolShake module will be the principal way our system moves and evolves, translating and rotating whole molecules but keeping the internal geometry of those molecules constant. The molecular dynamics run is critical for water (and, indeed, any molecular system) because it is this that will relax and evolve the intramolecular degrees of freedom in our molecules and give us realistic sampling of our bonds and angles. We don’t need this to happen every iteration of the simulation, so the frequency of the MD module is set to 5 by default.

g(r) and Structure Factors

We have neutron-weighted experimental data, so we need a layer to calculate RDFs and neutron S(Q):

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

First we need to set the instrumental broadening in the SQ module:

Click on the SQ module to display its options

Set the QBroadening to OmegaDependentGaussian with a FWHM of 0.02

This broadening is a known parameter of the SANDALS instrument on which the experimental data were collected.

For generality, some of Dissolve’s broadening functions refer to ‘omega’, which should be taken to mean the reciprocal space axis (in this case, ‘Q’).

Next, we will set up our calculation of the weighted structure factors. Since a NeutronSQ module calculates the partial and total structure factors for a single isotopic composition, we will need to add two more since we have three reference datasets.

Show the module palette for the layer by clicking Show Available Modules at the very bottom of the module list

Drag a NeutronSQ module from the list of available modules, placing it after the existing NeutronSQ module. Alternatively, double-click the NeutronSQ entry in the list of available modules to append one to the current list of modules

Add another NeutronSQ module, ensuring all three are after the SQ module

Note that each of the new NeutronSQ modules has a unique name (NeutronSQ, NeutronSQ01, and NeutronSQ02) - it is a requirement that modules within Dissolve can be uniquely identified by their name. We’ll now give the modules sensible names that describe our three datasets, and set the isotopologues and reference data files.

H2O

Click on the first NeutronSQ module (“NeutronSQ”) to display its options

Change its name from “NeutronSQ” to “H2O”

In the Reference Data settings group, for the Reference keyword select the file “SLS18498-H2O.mint01” and set the format of the data to mint

D2O

Click on the second NeutronSQ module (“NeutronSQ01”) to display its options

Change its name from “NeutronSQ01” to “D2O”

In the Isotopes & Normalisation section click the button for the Isotopologue option - it will currently say <Default to Natural>

Press the Species button to add a new isotopologue for each species present

Change the isotopologue from Natural to Deuterated

In the Reference Data settings group, for the Reference keyword select the file “SLS18502-D2O.mint01” and set the format of the data to mint

HDO

The HDO sample is a little different in respect of the isotopologue specification. In order to get the 50:50 mix we will basically add two isotopologues for the water species - one H2O and one D2O. Each will have the same relative weighting of 1.0. Importantly, we must also tell the module that the HW atom type is exchangeable - otherwise, the weighting of the intramolecular terms will be incorrect as, in effect, we will simulate a mix of distinct H2O and D2O molecules, which in reality is not what was measured since the hydroxyl hydrogens undergo fast exchange and appear as a mixture on the timescale of the neutron measurements.

As a general rule, any alcoholic or amine hydrogen is exchangeable.

Click on the third NeutronSQ module (“NeutronSQ02”) to display its options

Change its name from “NeutronSQ02” to “HDO”

In the Isotopes & Normalisation section click the button for the Isotopologue option - it will currently say <Default to Natural>

Press the Species button to add the natural isotopologue for the water species

Select the Water species

Click the Isotopologue button to insert another isotopologue (the next unused isotopologue will be added - in this case, the Deuterated one)

In the Isotopes & Normalisation section click the button for the Exchangeable option - it will currently say <None>

Tick the HW atom type in the list

In the Reference Data settings group, for the Reference keyword select the file “SLS18502-HDO5050.mint01” and set the format of the data to mint

Step 4 - Equilibrate the System

We’re now in a position to run our initial equilibration of the system, so let’s set the simulation running.

Simulation ⇨ Run

…or…

Ctrl-R

It’s useful to observe the energy of the configuration as we’re equilibrating our system, since once this reaches a steady value we can begin to think about what to do next. The Energy module provides plots of the various energy components:

Go to the   Evolve (Standard)   tab

Select the Energy module and click the Output button to see the module’s output graphs

You’ll see several curves on the graph, including the total energy (black line) which is the sum of the internal (intramolecular) bond and angle and the intermolecular (i.e. Lennard-Jones plus coulomb) terms. The MD module is set to only run if the total energy of the configuration is stable, but you should be able to notice when it kicks in as the bond (green line), angle (purple line), and total (blue) intramolecular energies which are constant during the Monte Carlo steps will suddenly drop to lower values. There will be some variation of the components up and down, but for a system of 1000 water molecules the total energy should approach a value of around -35×103 kJ mol-1 with the intermolecular energy around -44×103 kJ mol-1.

The energy will take a few hundred iterations to stabilise, so while we’re waiting we’ll take an initial look at the comparison between simulation and experiment in the next section.

Step 5 - Assessing the Reference Potential

With our equilibrated (or equilibrating…) system we’ll now make a basic comparison between our simulated total structure factors and the reference datasets.

Go to the   G(r) / Neutron S(Q)   tab

Click on the “H2O” NeutronSQ module and select the Output button

We’ll first check the agreement between the experimental total structure factor (the F(Q)) and the simulated data. You’ll see a row of buttons at the top of the graph giving you access to various different plotted quantities, with the default being the Total F(Q) which should look a little like this:

Calculated (black) and experimental (red) F(Q) for the equilibrated water (H2O) simulation

There are some small mismatches between the peak positions at $Q$ values of 3.1 and 5.1 Å-1 but overall the agreement is pretty good. The mismatch at $Q$ < 2 Å-1 is related to processing of the experimental data, and which we’ll not try to do anything to fix in this example.

Let’s take a look now at the total radial distribution functions:

Click the Total G(r) button on the “H2O” NeutronSQ module’s output page

Simulated (black), Fourier-transformed from simulated F(Q) (green), and experimental (red) G(r) for the equilibrated water (H2O) simulation

On this graph you will have three curves which we’ll discuss individually as the distinctions are important.

Calculated G(r) (black)

The black line represents the G(r) calculated directly from the atomic coordinates in the configuration, with broadening applied to the intramolecular terms (if requested).

Experimental G(r) (red)

The experimental G(r) is obtained by Fourier-transforming the reference data and as such is a subject to all the potential artefacts introduced by doing so. A window function is typically applied, and the operational range of the Fourier transform is often truncated in order to reduce ripple effects or to avoid known “bad” data at the lower Q end of the spectrum.

Simulated G(r) Calculated by Fourier Transform (green)

Since the experimental data is subjected to Fourier transform in order to generate a reference G(r), the most reliable comparison to make is to simulated data treated in the same way. Thus, the green line (“Via FT” in the plot) is obtained by double Fourier of the original calculated G(r) (black line), taking it to Q-space and back again.

So, any critical comparison between simulated and experimental G(r) should be made by comparing the red and green lines.

Step 6 - Fix the Water Molecule Geometry

There is a hint in the structure factors for the H2O sample (particularly the G(r)), that the intramolecular geometry of our off-the-shelf forcefield is not quite consistent with the real world. This is made clearly obvious if you look at the G(r) for the D2O sample:

Go to the   G(r) / Neutron S(Q)   tab

Click on the “D2O” NeutronSQ module and select the Output button

Click the Total G(r) button on the “H2O” NeutronSQ module’s output page

Simulated (black), Fourier-transformed from simulated F(Q) (green), and experimental (red) G(r) for the equilibrated water (D2O) simulation

Clearly we have a mismatch between the experimental (red) and simulated (green) peak positions at around 1 Å (related to the O-H bond) and 1.7 Å (related to the H-O-H angle). It is usually necessary to adjust the geometry of your species a little in order to be consistent with the experimentally-measured data, and in the present case we are lucky that we only have two parameters to adjust!

Here we are modifying the intramolecular terms based on comparison of the D2O data, but bear in mind that liquid water is amongst the systems most sensitive to isotopic substitution since all hydrogens are hydroxyl hydrogens, and subject to exchange as well as strong hydrogen bonding. As such, the differences in intramolecular geometry between H2O and D2O are measurable.[1]

Since we set up our simulation to use intramolecular master terms (via the Add Forcefield Terms… wizard) we can modify those to directly affect our water molecule’s geometry. For the O–H bond it is quite straightforward to read of the true distance (0.976 Å) from the reference g(r). The angle distance requires a touch more trigonometry but, given knowledge of the correct O–H distance, we can work out that the corresponding equilibrium angle we require is 107.134 °.

First of all, let’s stop the simulation from running:

Simulation ⇨ Stop

…or…

Esc

And now let’s make the changes to our intramolecular terms:

Go to the   Forcefield   tab and find the Master Terms tab

In the Bonds section, change the equilibrium bond length (’eq’) of the HW-OW bond term from 1.0 to 0.976 Å

In the Angles section the equilibrium angle (’eq’) of the HW-OW-HW angle term from 113.24 to 107.134 °

Now run the simulation for a little longer to let the species adjust to their new geometry:

Ctrl-R

You should see a marked improvement in the comparison of the D2O total G(r) and structure factor.

Simulated (black), Fourier-transformed from simulated F(Q) (green), and experimental (red) G(r) for the equilibrated water (D2O) simulation with adjusted intramolecular parameters

The change in the G(r) will not be instant as the majority of the evolution of the system is from the MolShake module which does not change the intramolecular geometry. Only the MD module will affect the intramolecular geometry. Also, the g(r) calculated by the GR are averaged over five calculations by default.

It’s also worth checking the other two samples, where the same kind of improvement should be noticeable (if a little less prominent).

References

  1. Quantum Differences between Heavy and Light Water, A. K. Soper and C. J. Benmore, Phys. Rev. Lett. 101, 065502 (2008).

Step 7 - Set up Potential Refinement

Let’s briefly recap what we’ve done so far:

  1. Set up a liquid water system based on a literature forcefield (SPC/Fw, see http://dx.doi.org/10.1063/1.2136877)
  2. Equilibrated the system and made an initial structural comparison with experimental data
  3. Adjusted the intramolecular geometry of the water molecule in order to better match the experimental data

Our agreement with experiment is relatively good, but it is possible to make it even better by modifying the inter-atomic interaction parameters contained in the atom types. However, generally this is not to be attempted by hand in all but the simplest of cases, as the effects of adjusting the interatomic are seldom as obvious as those for intra-molecular parameters. Also, for even a modestly-complicated system the number of adjustable parameters is simply too large to tackle with manual fitting.

Here we’ll employ the EPSR module in order to adjust the interatomic potentials automatically to give better agreement with the experimental reference data.

Esc

Layer ⇨ Create ⇨ Refinement ⇨ Standard EPSR

Our new layer contains only the EPSR module, and which Dissolve has set up with sensible targets and defaults. We’ll explore the various graphs as we proceed, but for now let’s check the set-up of the module. Brief descriptions of the important parameters are given below - for more in-depth explanations see the EPSR module page.

An initial value for EReq has been set (3.0) - this determines the magnitude or “strength” of the generated interatomic potentials

The Feedback factor is 0.9 - this states that we are 90% confident in the experimental data, and that the calculated partial structure factors should make up 10% of the estimated partials

The range of data over which to generate the potential in Q-space is determined by QMax (30 Å-1) and QMin (0.5 Å-1)

The experimental data to use in the refinement are set in the Target option, which lists all available modules by name that have suitable data for the EPSR module to use. You’ll see that Dissolve has added all of the available NeutronSQ modules have been selected automatically.

All of these default values are fine for our purpose, and there’s very little that you should have to change in the first instance. So, start the simulation running again to begin the refinement process:

Ctrl + R

While the simulation is running we’ll go through the different Output tabs in the EPSR module one by one to see what information we have available, and which help to illustrate the basic workflow of the EPSR methodology.

1. F(Q)

Basic comparison between reference experimental data (solid lines) defined as Targets and those simulated by Dissolve (dashed lines), including the difference functions (dotted lines).

Experimental reference F(Q), calculated F(Q), and difference functions for data targeted by the EPSR module

2. Delta F(Q)

Delta functions between experiment and simulation are replicated here along with the fits to the functions achieved by the Poisson / Gaussian sums. Note that the difference functions here have been multiplied by -1.

Delta functions (negated) and the Poisson/Gaussian fits to those functions

3. Estimated S(Q)

Estimated partial structure factors (solid lines) derived from combining experimental F(Q) and the calculated partial S(Q) (shown as dashed lines).

Estimated partial structure factors

4. Estimated g(r)

Partial radial distribution functions for each atom type pair in the system, showing estimated partials (solid lines, calculate by Fourier transform of the estimated S(Q)) and those from the simulation (dashed lines, calculated by Fourier transform of the simulated partial g(r)).

Note that the estimated g(r) from the combination of the experimental and simulated data are intermolecular partials - i.e. they contain no contributions from intramolecular terms related to the intramolecular bonding. The simulated partial g(r) (dashed lines), however, are the those calculated directly from the configuration(s) in the simulation, and include both intermolecular and intramolecular terms. The reason for this is to aid in detection of cases where the intramolecular structure of the molecules in the simulation is significantly different from that measured experimentally, which would manifest as small features in the estimated g(r) at or around the r-values associated to bond or angle (1,3) distances.

Estimated partial radial distribution functions

5. G(r)

Neutron-weighted total G(r) from Fourier transform of the reference F(Q) (solid line) and the Fourier transform of the simulated F(Q) (dashed line). As mentioned earlier, the reason for displaying the latter rather than the directly-calculated G(r) is to provide for a consistent comparison between simulation and experiment given the necessity of performing Fourier transforms (truncation errors are worse in Q-space data, so transforming to Q to r demands more judicious use of window functions in the transform).

Total radial distribution functions from Fourier transform of reference and simulated F(Q) data

6. Potentials

φ(r) represent the empirical potentials generated from the difference data, one per pair potential. These are applied on top of the reference potential.

Empirical potentials for each atom type pair (i.e. each distinct pair potential)

7. R-Factors

For each experimental dataset the “goodness of agreement” with simulation is represented by the r-factor, along with the total (average) r-factor over all datasets.

R-factors for experimental (reference) datasets

8. EReq

Here you can see the evolution of the current potential magnitude (whose limit is guided by the value of EReq).

Step 8 - Managing the Refinement

Now that our simulation is running and happily refining interatomic potentials, it’s fair to ask “How long should I refine for?”, or “What limit of EReq should I use?”, or even “How do I know when I’ve achieved the best agreement with experiment that I can get?”. The somewhat vague answer to all these questions is “When you can’t make the simulation any better.” Simply giving the EPSR module the freedom to do anything with a large EReq value is not a sensible approach, since it is highly likely that increasing the magnitude of the empirical potentials too much will harm the simulation rather than do it good. Unfortunately there is no direct way to determine what value of EReq is optimal a priori, and so some degree of trial and error is involved.

The EPSR algorithm implemented within Dissolve doesn’t contain any attempts at automating the adjustment of EReq, so a general procedure to follow is:

  1. Set a modest (i.e. low) value of EReq to begin with (the default of 3.0 is enough for neutral molecules, but should be higher for ionic systems).
  2. Run the simulation and observe the progression of the r-factors on the R-Factor tab of the module (hopefully they should decrease from their initial starting values).
  3. Once the r-factors have stabilised, increase EReq (increasing by an amount equivalent to the initial value each time is sensible).
  4. If the r-factors decrease again, repeat (3). If not, then you have probably reached the limit of what is possible within the context of the current simulation and its forcefield.

If you now take a look at the r-factors for your water simulation you should see that they have steadily decreased while you’ve been reading. In fact, an EReq value of around 3 represents the “limit of no improvement” for the present case - R-factors should be of the order of 1.4×10-4, and the simulated F(Q) should now agree pretty well with the experimental data:

Experimental (solid lines) and simulated (dashed lines) total F(Q) after refinement

We can observe some discrepancies at the low Q end of the data, but there is little that we can do to remedy these errors since they are likely due to imperfect subtraction of inelasticity effects in the processing of the original data. Simulated F(Q) (even un-refined ones) often give a good indication of what the experimental data “should look like”, especially in the region where inelasticity effects dominate, and can be useful in guiding data processing (e.g. with Gudrun).

At this point, we can start to think about what kind of structural properties we want to calculate from our refined simulation.

Step 9 - Setting up Analysis

Having now generated and refined a simulation so that it is consistent with experimental data, we can finally get round to the main goal of any total scattering experiment - understanding the structure of the system. Exactly what questions you want to ask depends on the what you really want to know about your system, and what analysis you will want to run on your simulation.

The process of refining a simulation against experimental structure factors means that we already have access to the atomic partial and total radial distribution functions. However, usually there are more interesting properties to know - for liquid water, we may wish to know the centre-of-mass radial distribution function (i.e. a g(r) between molecules) instead. We may want to quantify the geometry of the hydrogen bond interactions in the system, or we might be interested in the three-dimensional distribution of water molecules about a central, reference water molecule. The final “production” stage involves calculating these quantities and drawing some meaningful conclusions from them.

To achieve this we need to define one or more ‘sites’ on our water molecule. Sites are used in all of Dissolve’s analysis routines, and represent specific atoms, or averages of collections of atoms to be used as points of reference in the calculation of quantities. There are built-in analysis modules to calculate common quantities of interest, but Dissolve also allows custom analysis routines to be defined. Here, we shall only concern ourselves with Dissolve’s built-in functionality.

Defining Sites

Sites are defined within species and basically represent instructions for calculating a reference point (and possibly axes) for any molecule of that type in the simulation. Sites can be found on the relevant species tab:

Pause the simulation with Esc

Open the   Water   species tab and go to the Sites section

You’ll see in this section a list of defined sites (currently none) for your water species, along with a viewer of the species on the right which can be used to interactively select atoms for the site origin, or to define axes. It also illustrates the position and, if defined, the axes for the site.

Next Steps

The following subsections will describe how to set up several flavours of analysis in separate layers. Feel free to add all of them, or just a selection, to your simulation before moving on to Step 10 and starting it up again to see the results.

Step 9a - Centre-of-Mass RDF

Examining the centre-of-mass RDF between molecules is often useful since it gives a general picture of the arrangement between molecule. A related quantity is the coordination number, i.e. the number of molecules within solvation shells that typically form around molecules in the liquid. Here we’ll set up calculation of the centre-of-mass RDF and the associated coordination number for the first solvation shell. Dissolve provides a predefined layer to calculate both at once, however first we need to define a suitable site - the centre of mass of the water molecule.

Define the Centre-of-Mass Site

If you don’t already have a site COM representing the centre-of-mass of the molecule, perform the following steps:

Open the   Water   species tab and go to the Sites section

In the species viewer to the right, click-drag with the left mouse button and draw a box around all three atoms of the water molecule to select them.

Right-click the selcted atoms and click Create site from ⇨ Atoms

In the Sites section, under Definition select the Weight origin by mass option

Double-click on the new site in the list and change its name to COM

You’ll see the site represented in the species viewer as a 3D crosshair indicating the coordinates of the origin, surrounded by a small cube.

Create an Analysis Layer

We’ll now create one of Dissolve’s predefined layers to give us the processing modules we need:

Layer ⇨ Create ⇨ Analysis ⇨ RDF & Coordination Number

The new layer contains the following modules:

Module Purpose
SiteRDF Calculates the radial distribution function between two specified sites

We’ll need to set up both of these modules to calculate exactly what we need. First, SiteRDF:

Select the SiteRDF module and display its Options

Click the SiteA option and choose the COM site

Select the COM site for SiteB as well

Ensure that the ExcludeSameMolecule option is enabled

Since we are calculating the RDF of a site around itself, the ExcludeSameMolecule option prevents consideration of the same site with itself (which would give a distance of zero and result in a large undesirable spike at the beginning of the RDF). The distance range over which to calculate the RDF can be set in the Control settings group, but the defaults (0 - 10 Å with a step size of 0.05 Å) are fine for what we need.

For the coordination number calculation:

Tick the RangeAEnabled option

Set the maximum for RangeA to 3.4 Å

We can specify up to three separate ranges over which to calculate coordination numbers but we will focus on the first coordination shell in this example.

Step 9b - Hydrogen Bond Geometry

Water is the archetypal hydrogen bonding fluid, and so it makes sense to analyse the hydrogen bonding contacts within the fluid. We will be principally interested in the geometry of O–H···O contacts occurring in the liquid in terms of the H···O distance and the O–H···O angle.

Define Atomic Sites

We’ll need simple atomic sites for the oxygen and each of the hydrogen atoms in our water molecule, so let’s create them now:

Open the   Water   species tab and go to the Sites section

Click the Oxygen atom in the species viewer to select it

Right-click the selcted atom and click Create site from ⇨ Atoms . This will create a site called O

Repeat this for each hydrogen atom individually, naming those sites H1 and H2

Create the Analysis Layer

Let’s now create a new layer and add the analysis module that we want:

Layer ⇨ Create ⇨ Empty

Rename the layer by double-clicking on the tab name, and change it to ‘Analyse HBond’

Show the module palette for the layer by clicking Show Available Modules at the very bottom of the module list

Drag a DAngle module over to the Current Modules list

As its name suggests, the DAngle module calculates distance and angle histograms, and their 2D map, from three target sites. Since the module is general-purpose, we’ll need to set up the module to give us information on our specific geometry of interest. DAngle assumes in the target interaction A–B···C that A–B occur on the same molecule, and gives us a choice as to whether we exclude C sites that are also on the same molecule.

Select the DAngle module to display its Options

Change the maximum value of the DistanceRange to 5.0 Å

Press the button for SiteA and choose the O site

For SiteB select both the H1 and H2 sites

For SiteC choose the O site

Enable the ExcludeSameMolecule to ignore interactions where site C is on the same molecule as A and B.

Step 9c - Centre-of-Mass Spatial Density

Radial distribution functions provide a spherically-averaged picture of structural correlations between sites, so when there are strongly directional interactions (such as hydrogen bonds, for instance) or simply where the central molecule of interest lacks symmetry, understanding the distribution of other sites in 3D space around a given molecule type is much more informative. A spatial distribution function or spatial density function (SDF) tells us, in a manner analogous to the radial distribution function, the densities of some object in 3D space relative to the bulk density. In other words, given a central molecule, the SDF tells us where other molecules or sites prefer to locate themselves.

The SDF is an “orientation-dependent” function - rather than averaging spherically around a molecule as in the RDF, the SDF preserves spatial information relative to some reference set of axes related to the local geometry of the central molecule. So, our first task is to define a set of reference axes for our water molecule.

Set Up Site Axes

We will define axes on our centre-of-mass site. In general, a system of reference axes can be constructed quite easily from the following steps:

  1. Define a coordinate that will represent the origin of the axes (i.e. the local coordinate {0,0,0}). In the present example this will be the centre-of-mass of the molecule.
  2. Select a group of one or more atoms whose coordinates, when averaged, will define absolutely the direction of the x axis.
  3. Select a group of one or more atoms whose coordinates, when averaged, will define the rough direction of the y axis (and more importantly the XY plane), and which will be orthogonalised with respect to x (i.e. made to form a 90° angle).
  4. Define the z axis from the cross product of the x and y axes.

When defining a site within Dissolve it is your responsibility to provide atom indices for the origin and the representative x and y directions - Dissolve will do the rest. Thus:

Open the   Water   species tab and go to the Sites section

Make sure the COM site we created in Step 9a is selected in the list on the left

Select the oxygen atom

Right-click the oxygen atom and click Modify current site ⇨ Set X-Axis Atoms

Select either one of the hydrogen atoms

Right-click the selected hydrogen atom and click Modify current site ⇨ Set Y-Axis Atoms

You’ll note that, as soon as you pressed the Y Axis button a set of axes appeared in the site viewer, letting you know that your definition is complete and showing you how your system of axes are oriented.

Create the Analysis Layer

We will add another of Dissolve’s predefined analysis layers to get the modules that we want:

Layer ⇨ Create ⇨ Analysis ⇨ Average Molecule & SDF

Along with the calculation of the SDF, it is useful to have the actual average geometry of the species around the specified site, and so the layer provides the following modules:

Module Purpose
AvgMol Calculates the average geometry of the parent species around a specified reference site
SDF Calculates the spatial density function of a particular site around a central reference site

As the name implies, the AvgMol module will provide us with the average species geometry which we can use as a reference when we plot the spatial density function calculated from the SDF module. The SDF module takes two sites as input - the central site A about which the distribution of site B will be calculated. For both modules, the principal (central) site must have a system of axes defined.

Let’s proceed and set up the two modules in the layer. First, AvgMol:

Select the AvgMol module to display its Options

Find the Control settings group, and click the Site option.

Set the Species to Water and the Site to COM for the average molecule calculation (it will be the only one available, since it is the only one which has a set of axes defined)

And now SDF:

Select the SDF module to display its Options

In the Control settings group set the central SiteA to COM

Set the surrounding SiteB to O

Step 10 - Running the Analysis

With our analysis modules set up we are now ready to restart the simulation and begin calculating the quantities we requested - we’ll call this our “production run”. There are a couple of things to note at this point:

  1. The EPSR module / layer should remain enabled during the production run in order to maintain the simulation in its present state.
  2. All quantities calculated from analysis modules are averaged automatically over time.
  3. Simple quantities (e.g. radial distribution functions) tend to accumulate statistics quickly (in a small number of iterations) while more complicated quantities (e.g. spatial distribution functions) will take much longer to give a useful average.

When you’re ready, restart the simulation and we can look at the quantities start to form as the data accumulates. Brief overviews of what to expect are given below.

Ctrl + R

Centre-of-Mass RDF

The SiteRDF module output shows a strong peak at around 2.8 Å for the COM–COM distance in the liquid. Integrating the curve up to the first minimum (approximately 3.4 Å) gives us a first shell coordination number of around 4.7 molecules (displayed in the panel below the graph).

Is the coordination number of 4.7 realistic?

Hydrogen Bond Geometry

Go to the   Analyse HBond   layer

Click the DAngle module and go to the Output page

There are three graphs on display here - the B···C g(r) (top left), the A–B···C angle histogram (bottom left), and the 2D correlation map of the two (right). The B···C g(r) corresponds to the intermolecular H···O and you should see the principal hydrogen bonding peak at 1.8 Å followed by a second, broader interaction at 3.3 Å. So, the primary hydrogen bonding distance in the liquid is 1.8 Å, but how directional is this interaction? If we look at the distance-angle map and read off the x axis to 1.8 Å we can see that the A–B···C angle is always between 150 and 180°. In other words, a very directional hydrogen bond! Angles associated to the broader feature at 3.3 Å are around 60°.

What causes the broad feature at feature at 3.3 Å / 60°?

3D Liquid Structure

To make sense of the spatial distribution function output by the SDF module we will need to set the central reference molecule to the average molecule we also requested, and adjust the cutoff to get a sensible surface.

Esc

Go to the   Analyse AvgMol/SDF   layer

Click the SDF module and go to the Output page

Set the Reference Molecule to COM@Water (AvgMol)

Change the Lower cutoff to 0.065

Cutoffs for the surface display are given in terms of the number density per cubic Å of the surrounding molecule / site in the simulation box. In the present example we have 1000 water molecules in a cubic box of side length 31.0723 Å, so the number density of water molecules is 1000 ÷ 31.07233 = 0.033. The cutoff we will set below therefore gives a surface that represents twice the bulk density of water molecules, i.e. a higher than normal probability of finding water molecules.

Spatial density functions take a relatively long to accumulate to a point where the surfaces are smooth, but you should be able to make out the salient features after only a handful of iterations. Along each O–H bond vector there exists a small lobe of density, and represents the positions of molecules acting as hydrogen bond acceptors to the central one (this visually represents the distance and angle we found in the DAngle results above). A halo of density also exists around the oxygen atom, representing the positions of molecules acting as hydrogen bond donors to the central one.