View Single Page Version

# Bulk Silica

30 minutes

Summary

Amorphous silica (SiO₂) is found in a variety of places, and is a relatively simple, network forming glass. Here we’ll set up a bulk silica simulation in a unit cell of our own design, and which can serve as the basis for other tutorials.

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: SiO2_NIMROD.mint01

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

Step 1 - Set up Atom Types

We’ll be setting up our silica system using atomic (ionic) species, rather than creating an explicitly bound network. So, we need two atomic species for this simulation - an oxygen, and a silicon. First, however, we’ll create some suitable atom types by hand - this is instead of applying a built-in forcefield or loading one in from elsewhere.

Select the Forcefield tab

In the Atom Types group, click the Add button located at the top-right

Select silicon in the periodic table dialog

Atom types created in this way are named after the symbol of the relevant element.

This creates an empty atom type in the table, and which we’ll now set some useful interaction parameters.

Change the “Charge” of the new atom type to 2

Change the “SR Form” to LJ to request a Lennard-Jones short-range potential for the atom type

In the “SR Parameters” box enter epsilon=0.175 sigma=1.03

We’ll now use the same process to create an oxygen atom type:

Click the Add button again, and choose oxygen from the periodic table

Change the “Charge” of the oxygen atom type to -1

Change the “SR Form” to LJ

Set the “SR Parameters” to epsilon=0.165 sigma=3.5

You might ask where these values for the short-range interaction parameters come from, and why the charges don’t reflect what you might expect to be “formal” oxidation state charges for the ions. Well, those are both good questions, and the real answer is that there are going to be a variety of combinations of Lennard-Jones epsilon, sigma, and atomic charges that give reasonable structures for silica, and those listed here are good enough for a starting point for us. There are no doubt other forcefields in the literature which suggest different parameters for an atomic/ionic representation of the system - you are free to try those if you wish!

Step 2 - Set up Atomic Species

Now let’s generate our atomic species, assign the atom types that we created in the last step, and create a basic site for each so we can do some analysis afterwards. First, the silicon:

Select silicon from the periodic table

In the Atoms section, change the AtomType (AT) for the atom to Si

With the atom selected in the table, right click the atom in the viewer and click Create site from ⇨ Atoms

And now, the oxygen:

Select oxygen from the periodic table

Go to the Forcefield tab

In the Atoms section, change the AtomType (AT) for the atom to O

With the atom selected in the table, right click the atom in the viewer and click Create site from ⇨ Atoms

Step 3 - Set up a Configuration

To create our silica configuration we’ll use the “relative random mix” generator as a basis.

Choose both species and press Next

Leave the configuration type as Mixture and press Next

Leave the box style as Fixed Geometry, Undefined Size but set the angle γ to 120 degrees (this will give us a monoclinic box) and press Next

Now we need to set up a relative mix of our two species:

In the table change the Population / Ratio of oxygen from 1 to 2

Set the Density of the box to 0.066 - note that this is in units of atoms Å^-3

Set the Multiplier to 1500

Press Finish to complete the wizard.

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

Note that We use units of atoms Å^-3 here as this allows us to add both atomic components using the same density value. If we had chosen units of g cm^-3 we would have to know the individual densities of silicon and oxygen in amorphous silica. Anyway, if you zoom the viewer out (mouse wheel) then you should then see a monoclinic box of silica!

Step 4 - Set up Processing

Let’s create two processing layers now - one to evolve our configuration, and one to calculate the current RDFs and structure factors.

We have only atomic species, so we’ll choose the basic atomic evolution layer for our system:

Our reference data is neutron-weighted so:

We just need to set the instrumental Q broadening in the SQ module and assign the reference data to our NeutronSQ module.

First, the Q-broadening:

Click on the SQ module to display its Options

Set the QBroadening to OmegaDependentGaussian with a FWHM of 0.03

Now for the reference data. Dissolve can perform some rudimentary operations on datafiles when reading them in - for our present experimental reference data (“SiO2_NIMROD.mint01”) it contains a point at Q = 0.0 Å^-1 which we will trim off by setting the minimum x (Q) value to read.

Click on the NeutronSQ module to display its options

For the Reference keyword in the Reference Data settings group, select the file SiO2_NIMROD.mint01 and set the format of the data to mint

Open the settings for the Reference datafile by clicking the Options button near to the file type selector, and in the Manipulations tab set XMin to 0.1

Step 5 - Equilibrate the Box

To ensure that our silicon and oxygen atoms are well mixed and distributed we’re going to “cook” the box at a high temperature to start with.

Go to the Bulk configuration tab

On the configuration viewer toolbar, click the Temperature button

Change the Temperature to 2000.0

Now we’ll run for 250 steps:

Set the number of steps to 250 and click OK

You can open the NeutronSQ module in the RDF / Neutron S(Q) layer to monitor the structure, or the Energy module in the Evolve (Basic Atomic) to check the energy as the simulation progresses. The total energy will be large and negative - of the order of –4×10⁶. When finished, the F(Q) and G(r) should look a little like this:

Calculated structure factor (black line) of amorphous silica at high temperature (2000 K) compared to reference data at 298 K (red line)

Calculated total G(r) (black line) of amorphous silica at high temperature (2000 K) compared to Fourier transform of reference data at 298 K (red line)

We’ll now reset the temperature of the box to 298 K and run a second equilibration.

Go to the Bulk configuration tab

On the configuration viewer toolbar, click the Temperature button

Change the Temperature to 298.0

Set the number of steps to 250 and click OK

Once finished, the structure should look a lot better:

Equilibrated structure factor (black line) of amorphous silica compared to reference data at 298 K (red line)

Equilibrated total G(r) (black line) of amorphous silica compared to Fourier transform of reference data at 298 K (red line)

Add the standard EPSR-style refinement layer to your simulation:

Select the EPSR module if it isn’t already. We’ll need to set a relatively high starting value for EReq as we are basically dealing with point charges in a box, and so whatever empirical potential is necessary to drive the structure of the system needs to grow large enough to have some effect against the strong Coulombic forces present.

Go to the Refine (EPSR) tab

In Options go to the Control settings group and set EReq to 50.0

If we run the simulation a bit longer:

On the R-Factor tab of the EPSR module you will see that the magnitude of the potential (top graph) increases rather quickly to 50 (within 100 steps in the plot below), but the quality of fit as measure by the r-factor (bottom graph) takes rather longer to reach a steady state.

Evolution of empirical potential magnitude (max = 50.0)

Evolution of r-factor (EReq max = 50.0)

The structure factor of the system at this point looks something like this:

Stable structure factor at EReq = 50.0

The question is, can it be better? The answer is, of course, yes - the intensity of the first strong peak in the F(Q) at 1.47 Å is a little low. The empirical potential will eventually sort this out (perhaps with an increase of EReq to 75 or so) but it takes time to get it exactly right.

What does the peak at Q = 1.47 Å correspond to, and why might it be difficult to reproduce?

For the purposes of the example, let’s move on and calculate something from the (semi-) refined structure.

First, stop the simulation:

Esc

Step 7 - Analyse the Structure

We’ll introduce another analysis module, Angle, to finish off this tutorial. Angle is similar to DAngle, but is more flexible, assumes less about the connectivity between the sites, and gives more comprehensive output. We’ll use it to calculate the Si–O–Si angle, and the correlations of that angle with the individual Si–O bonds.

Create an empty layer, and add a Angle module to it.

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

Drag a Angle module from the Analysis group in the Module Palette into the list above to add it to the layer

Now let’s set it up.

Select the Angle module to display its Options

Change the RangeAB Min and Max values to 1.2 and 2.1 respectively, and corresponding to the range of the Si–O bonds

Change the RangeBC Min and Max values to 1.2 and 2.1 as well

Change the AngleRange BinWidth to 5.0

Select the Si site for SiteA and SiteC

Select the O site for SiteB

Enable the ExcludeSameSiteAC option, as otherwise we’ll get a spike at zero in the resulting angle data

Now we can start Dissolve running again and look at the calculated data:

Ctrl-R

Select the Angle module and open the Output page

The three graphs on the first page of Angles tabs shows the two radial distribution functions A–B and B–C, and the A–B–C angle. The RDFs should mirror those calculated by the NeutronSQ module (peak at around 1.6 Å) while the normalised angle graph tells us that all the Si–O–Si angles are predominantly greater than 120°.

The second and third tabs display us the distance-angle correlation graphs, which are both identical in the present case.

View Single Page Version

Summary

Preparation

Step 1 - Set up Atom Types

Step 2 - Set up Atomic Species

Step 3 - Set up a Configuration

Step 4 - Set up Processing

Step 5 - Equilibrate the Box

Step 6 - Set up Potential Refinement

Step 7 - Analyse the Structure