View Single Page Version
20 minutes
Summary
Liquid argon isn’t a particularly exciting system, but it has been measured experimentally by neutron scattering, and it serves as an excellent example to introduce some of the basic workflows of Dissolve. This example will show you how to set up a complete working simulation from scratch, with the purpose to simulate liquid argon at 85 K and compare the simulation to reference neutron data from “Structure Factor and Radial Distribution Function for Liquid Argon at 85 °K”, J. L. Yarnell, M. J. Katz, R. G. Wenzel, and S. H. Koenig, Phys. Rev. A 7, 2130.
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.
- Neutron data file: yarnell.sq
File ⇨ New
File ⇨ Save As…
Save your input file under a sensible name in same directory as the example data
Step 1 - Creating a Species
We need a suitable species to represent the argon atom in our simulation.
Species ⇨ Create ⇨ Atomic…
Select argon (Ar) from the periodic table
You will now see that a new tab called Ar has appeared containing the new species. In the Atoms section on the left you can see the single argon atom located at (0,0,0). and underneath that the summary of intramolecular terms holding your species together which, for a single atom, isn’t a lot! Our argon atom currently has no atom type assigned to it (the “AtomType” (AT) column in the table is empty), but we will remedy this in the next step. All atoms must be assigned an atom type before the simulation can proceed.
An atom type contains descriptive interaction parameters for a particular type of atom, and optionally an associated atomic charge. They describe the core interatomic interactions in your system, and define the pair potentials between atoms.
Note also the presence of the Isotopologues section which allows us to define isotopes for specific atom types, and which we’ll return to later.
Isotopologue definitions depend on atom types, so all atom types should be defined before setting up isotopologues.
Step 2 - Set up a Forcefield
To simulate any system we need to supply suitable parameters to describe the interactions between atoms, both internally within the same species and externally between atoms of the same (and different) species. The Forcefield tab is always available, and summarises the atom types used over all defined species and the resulting pair potentials, as well as master intramolecular definitions (which we aren’t using in this example). You can also set other important options, such as the range of the generated pair potentials and the truncation schemes.
We’ll now describe the interactions for your argon species by taking the terms from one of Dissolve’s built-in forcefields.
Go to the Ar species tab
By selecting the tab containing the Argon species you are making it ‘active’, and so it will be the target of any actions triggered from the Species menu.
Species ⇨ Add Forcefield Terms…
From the available forcefields choose
OPLSAA2005/NobleGases
and click
Here we choose how to assign atom types to the species - the default option of Determine atom types for all atoms uses connectivity descriptions in the forcefield to try and automatically choose which type to use.
Leave this option selected and click
There will be no conflicts between the proposed atom types and existing ones defined in the main simulation (since there aren’t any) so click
We will let Dissolve apply intramolecular terms to the whole species (there aren’t any anyway), so leave this section as it is and click
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
to exit the wizard
For more complicated (molecular) species the wizard also handles how intra-molecular terms are generated. You also may have to deal with potential naming conflicts with existing terms in the main simulation, which the Add Forcefield Terms wizard will help you resolve.
If you now return to the
Forcefield
tab you will see in the Atom Types section that we have an atom type for our argon atom which contains the necessary interaction parameters. The short range type (SR Form in the table) specifies the functional form of the van der Waals interaction parameters for each atom type, which in our case is LJGeometric
indicating that the parameters reflect a standard Lennard-Jones potential utilising geometric combination rules. Those parameters (epsilon
and sigma
) are shown in the SR Parameters column and describe how our argon atoms will interact with each other in the simulation.
Atom types are strictly associated with a chemical element, and can only be applied to atoms of the same element.
By default, Dissolve will generate all the necessary pair potentials for the current set of atom types automatically - these are listed in the Pair Potentials section where you can also select individual potentials and visualise them in the associated plot.
The current Pair Potentials range is set too high for our box size. We must alter this value to be less than the radius of largest possible inscribed sphere for box:
In the Control group change the Range to
10.0
.
Before we move on we need to set an option related to the charges to use. Along with van der Waals parameters, atomic charges are the second important component in describing interactions between atoms. In Dissolve these charges can either be taken from the atom types and included in the tabulated pair potentials, or located on species atoms and calculated analytically. Dissolve automatically detects the most appropriate scheme but is wary of systems or species where there are too many atoms with zero charge, as is the case here. As such, we must force a choice:
In the Charge Source group at the top-right of the Forcefield tab deselect the Choose Automatically option and enable the Force Choice option below it.
The actual charge source (Atom Types vs Species Atoms) is unimportant since charges are zero in both, but the default option is faster.
Step 3 - Add a Configuration
A configuration represents the target system using one or more of the available species. Our next step is to create a configuration that is consistent with the experimentally-measured system investigated in the original paper.
Configurations in Dissolve are created from a sequence of instructions that define box shape and size, the populations of species to add etc. This permits Dissolve to automatically regenerate configurations when it needs to, and is particularly useful when screening structures with variable parameters (density, composition, pore size etc.).
Dissolve has a wizard to guide you through the process of preparing standard configurations:
Configuration ⇨ Create…
Target Species
The first step of the wizard is to select all of the species which will go in to it - we only have one (our argon atom) so highlight it in the list and press
Configuration Type
The next choice is which type of configuration to generate. The available choices depend on the types of species selected - we can only make a random mix from the argon species, so press
to continue.
Box Geometry
With our basic type of configuration chosen we must now decide on the size and geometry of the periodic box (unit cell). The two options are Fixed Geometry, Undefined Size (where the final size of the box will depend on the how much and at what density we add in to it) and Fixed Geometry and Size where we set both explicitly now. We will use the former choice here.
When selecting Fixed Geometry, Undefined Size the axis lengths given are relative since they will be scaled as species are added to the box. Setting them all to
1.0
and leaving the angles at90
(the default) results in a cubic box which should be appropriate for most simulations. Press to continue.
Species Information
On the final page of the wizard we set the populations of the species and, if we chose Fixed Geometry, Undefined Size, the required density of the system. Change the Density to
0.0213
atoms/A3, then press to complete the wizard.
You will now see that the
section has been populated with several steps, and a new configuration of atomic coordinates has been created based on these steps. We’ll go through these one-by-one in the order they appear. To see the settings for any step, click on it to show its options in the panel immediately below the list of steps.Node | Purpose |
---|---|
Temperature | Defines the temperature of the configuration. |
Parameters | Defines one or more numerical parameters that can be referenced by other nodes in the generator. We have only one value called “rho” which defines the density of the configuration. Note that “rho” is just a number and has no associated units - these are specified in the Add node. |
Box | Defines side lengths and angles for the periodic box. Note that relative lengths can be given, and the box expanded automatically (when adding species, for instance) to achieve some specific density. |
Add | Adds a number of copies of the specified species to the configuration. The BoxAction option controls how/if the box volume should be scaled to give the density specified once the species has been added, and is most useful when supplying relative box lengths. If set to None the box would remain at its current size (defaulting to 1.0 Å3) but the full population of molecules will still be added. |
In the Add node note how we have referenced the “rho” parameter for the density.
Many numerical options can be given in the form of equations referencing variables such as those set in a Parameters node. A green tick indicates if the equation is valid.
Finally, we need to set the correct temperature for the configuration.
On the configuration viewer toolbar, click the
button
Set the Temperature to 85 K to match that of the experimental measurement
You can recreate a configuration at any time by clicking the
button, but bear in mind that other quantities calculated by modules may also need to be cleared. Remove all of this data to begin a ‘clean’ run with the Clear Module Data option in the Simulation menu.
Step 4 - Set up the Simulation
It’s time to tell Dissolve exactly what we want to do with this argon system. All of the useful work in a simulation achieved by sequences of modules, with each module performing a specific task, function, or calculation. Modules exist within ’layers’, and it usually makes sense to break up simulations into distinct parts (e.g. keeping evolution of the system separate from refinement or analysis) as layers can be turned on and off as the situation requires. Each layer also has a frequency at which it will run, relative to the main simulation loop counter - a frequency of 1 means the layer will be executed every iteration. Modules within the layer also have an associated run frequency, but this is relative to the number of times the layer has been executed, rather than the main loop counter.
We’ll split our simulation up into two layers:
- An ’evolution’ layer which moves the contents of our configuration around and calculates the total energy of the system
- A ‘calculation’ layer which periodically calculates the radial distribution function and neutron-weighted structure factor
You can have as many layers in a simulation as you need - whatever makes the most sense in what you’re trying to do.
Step 4a - Evolve the System
We almost always want to move the contents of our configurations around, and this is the job of the evolution layer.
Layer ⇨ Create ⇨ Evolution ⇨ Standard Atomic (MC/MD)
Our new layer contains the following modules:
Module | Purpose |
---|---|
AtomShake | Performs standard Monte Carlo moves on individual atoms |
MD | Performs a number of molecular dynamics steps, evolving the system according to Newton’s equations of motion |
Energy | Calculates the total energy, including contributions from intramolecular terms if present |
Selecting any module in the list on the left will show its full options in the panel on the right - select the AtomShake to see its list of options. AtomShake performs standard atomic Monte Carlo on one or more target configurations. It tries to move each atom in the configuration in turn, randomly displacing it by a certain amount, and checking the energy to see if that displacement constituted a ‘good’ move. As such it has options like TargetAcceptanceRate, which governs the number of moves that can / should be accepted in any one cycle, and StepSize which indicates the maximum distance by which an atom can be moved in any translation attempt (over time this will change automatically to give the acceptance ratio requested).
Most modules need configuration targets to work on, and all available configurations (one in our case) will have been automatically set as a target in each of the new modules in the layer. As mentioned, each module has its own defined frequency at which it will run within the layer - AtomShake and Energy modules will run every time the layer is run (frequency = 1) while the MD module will only run every fifth step. Each module can be individually enabled / disabled using the slider next to the frequency box, with a green tick indicating that the module is active.
Finally, modules may have entire control panels of other functionality, most commonly graphing output for various properties etc. This can be accessed by selecting the
button at the top of the controls. Go back to the configuration options for the module by selecting the button.All of the default settings for the modules within the evolution layer are sensible enough for our needs, so take a look around at what’s there if you want to, but there’s no need to change anything.
Step 4b - Calculate g(r) and S(Q)
The experimental data we’ll be fitting to is neutron scattering data, so we will need to calculate neutron-weighted structure factors.
Layer ⇨ Create ⇨ Correlations ⇨ RDF and Neutron S(Q)
The new layer contains the following modules:
Module | Purpose |
---|---|
GR | Calculates partial g(r) between every pair of atom types, and sums them into the total G(r) |
SQ | Calculates partial S(Q) between every pair of atom types by Fourier-transforming a source set of g(r), and sums them into the total F(Q) |
NeutronSQ | Takes the S(Q) calculated by an SQ module to generate the neutron-weighted partial and total structure factors |
Note that the ordering of these three modules in the layer is important - calculating the neutron-weighted structure factors requires the unweighted structure factors from the SQ module, which in turn requires the radial distribution functions from the GR module. Modules are run “top to bottom” in the list, fulfilling the requirement that the GR module runs before the SQ module etc.
We will now need to set a few parameters in the NeutronSQ module, in particular informing it of the isotopic composition of our system and loading in reference data.
A NeutronSQ module calculates S(Q) and F(Q) for a single isotopic composition. Where there are multiple isotopic samples, you need to add a NeutronSQfor each distinct isotopic composition.
Set up Isotopes
The NeutronSQ module will use isotopic natural abundances to calculate the neutron weights for all species unless we tell it otherwise. We’ll first define the correct isotopologue for our argon species, and then tell NeutronSQ to use it. The experimental measurement was made using Ar36 since its coherent scattering cross-section (24.9 fm) is considerably higher than that of the naturally-occurring mix (1.91 fm).
Go to the Ar species tab
Click on the Isotopologues section
Click
to create a new isotopologue definition assigning the default (natural) isotope to each atom type present in the species
Change the entry for the Ar atom type from
Natural (bc = 1.909 fm)
to36 (bc = 24.9 fm)
For sanity’s sake, you may also want to double-click on the name of the isotopologue in order to change it to something more meaningful (‘Ar36’ perhaps)
Now we’ll go to our calculation layer and set the isotopologue for our NeutronSQ module:
Go to the G(r) / Neutron S(Q) layer tab
Select the NeutronSQ module and open its main
The Isotopologue keyword currently shows that all species will “Default to Natural” isotopologues
Click the button for the Isotopologue keyword to open its full options
Press the
button to populate the list with the default isotopic selection for each species
Change the isotopologue for the argon species from
Natural
toAr36
(assuming that you changed the name earlier!)
The ‘Natural’ isotopologue for each species is defined internally by Dissolve, and is always available. It does not appear in the list of defined isotopologues on the species tab.
Import Reference Data
The NeutronSQ module itself looks after any related experimental reference data that we might wish to compare our calculated data to, and which we’ll now set up.
Go to the G(r) / Neutron S(Q) layer tab
Select the NeutronSQ module and look for the Reference Data settings group
For the Reference keyword open the file
yarnell.sq
, and set the format of the data toxy
The data, along with its Fourier transform, can now be seen in the module’s
page. You may notice that the data have been normalised to the average squared value of the atomic scattering and oscillate around 1.0 - we will need to tell Dissolve to convert the data back to absolute units and make it oscillate around zero.
Select the NeutronSQ module and find the Reference Data settings group (go back to the page if you need to)
Set the ReferenceNormalisedTo style to
SquareOfAverage
This tells Dissolve that the data have been normalised, and allows Dissolve to remove that normalisation in order to get the data in the correct units
Open the options for the file import at the extreme right of the ReferenceData keyword
Find the Manipulations group and set the RemoveAverage value to
9.0
This will instruct Dissolve to work out the average value of the data from a Q of 9.0 Å-1 onwards, and subtract it from the data
Step 5 - Equilibrate the System
With our two processing layers set up we can start to run the simulation and monitor how the system is evolving. We’ll do this from the NeutronSQ module.
page of the
Go to the RDF and Neutron S(Q) layer tab
Select the NeutronSQ module and open the page
We’ll visually compare the simulated total structure factor to the experimental one
Simulation ⇨ Run
You can also use
Ctrl-R
to start a simulation running
Note the counter towards the right-hand side of the status bar at the bottom of the main window which tracks the current iteration, and the status indicator to the far left of the status bar telling you what Dissolve is doing (and also note that pressing
Esc
stops the current simulation).
While the simulation is running you cannot edit any input values, keywords etc., but you can investigate the simulation’s progress and output as it happens. For example, you could go to the Standard Atomic (MC/MD) evolution layer and look at the page of the Energy module to see what the total energy of the configuration is doing.
After the simulation has been running for a little while (perhaps 100 iterations), you’ll see that the calculated data compare quite favourably with the reference data, with the G(r) and F(Q) looking something like this:
Equilibrated total G(r) for liquid argon
Equilibrated total F(Q) for liquid argon
When there is no further dramatic change in the calculated data you can stop the simulation:
Simulation ⇨ Stop
Or…
Esc
Keep in mind that the simulation will not actually stop until the current iteration is completed - most parts of the GUI will remain grayed out until then.
Step 6 - Next Steps
Now that our simulation is equilibrated, we will consider the quality of our calculations against the reference experimental data. The reference potential which we took from an off-the-shelf forcefield actually does an excellent job of reproducing the data, so there is very little else to do in terms of refining the potential in order to get better agreement. We also have limited options for analysing further the structure of the system - for a single-component atomic fluid the radial distribution function pretty much tells you everything you need to know.
You’re now prepared for more useful / interesting / exciting (delete as appropriate) examples - a good next step would be to try the liquid water example.