Introduction to Molecular Modelling:

Alumina surfaces and nanoparticles

1. Aims and Objectives

The aim of this task is to introduce the use of molecular mechanics for modelling solids and their surfaces. You will learn how to build a simple model of a periodic ionic solid (alpha alumina, Al₂O₃) using a package called METADISE, and use this to calculate the structures and energies of the bulk, surfaces and nanoparticles.

The code uses simple interatomic potentials to calculate the forces between the atoms and energy minimization to find the most stable structures.

By the end of this task you will be able to:

1. Generate the unit cell of a crystal structure given its cell dimensions, space group and atomic coordinates.

2. Appreciate what is required to define the potential interaction between ions.

3. Calculate surface energies, simulate morphology and construct some nanoparticles.

4. Evaluate the segregation of impurities to the different surfaces and investigate the effect on surfaces and morphologies

5. Optional: evaluate the effect of changing the extent of hydroxylation on surfaces and to set up files for molecular dynamics

2. Basic preliminary steps

This tutorial is on Alumina and you will have different tasks (different runs in each Step1, Step2, etc) so it’s recommended that you create a different folder for each calculation performed.

3. Molecular mechanics calculations using METADISE

3.1 Generating structural input files for METADISE

The structure of a three-dimensional solid requires the input of three main sets of information: the unit cell, the fractional coordinates and types of the atoms, and finally the space group symmetry.

Start by logging into one of the Windows PCs, and open the Network Drive labeled H: (this is the fileserver and saving any files you create here means they will be available in future on any of the University PCs). Create a new folder called ‘Alumina’. Double-click to enter the folder. Go in the folder Alumina, then create the folder `Step1´:

In Alumina/Step1, create a new file called al2o3.txt by first creating a new text file called al2o3.txt (Right-click > New > Text Document)

We want to add an Al ion at fractional coordinate (0, 0, 0.35203) and an O ion at fractional coordinate (0.3064, 0.0, 0.25) by typing the following

fractional

Al core 0.0 0.0 0.35203

O core 0.3064 0.0 0.25

Ends

Next add the information about the unit cell. Al₂O₃ is a hexagonal cell and has a cell parameter of 4.7589 Å on the a/b direction and a height of 12.9910 Å along c. This information goes above the fractional key word. Finally we need to add information about the space group, R3CH in this case (you may already know about space groups, but if not then please accept this for now). The corresponding number of the space group is 167. There is a table on the numbers of space groups at http://img.chem.ucl.ac.uk/sgp/mainmenu.htm or http://www.cryst.ehu.es/cgi-bin/cryst/programs/nph-table. At this stage the input file should look like this:

cuto 1.0 15.0 0.7 60 # convergence to angstrom ; short range cutoff

# core shell separation (if core-shel<0.7,

# considered as the same atom) ;

# maximum deviation of the angle ;

cell 4.7589 4.7589 12.9910 90.0 90.0 120.0 #Values of a, b, c, a, b, g for Al₂ O₃ unit cell

Space full R3CH 167 1 #R3CH 167 is space group and number

#1 for origin, can replace full for redu for reduced cell

Fractional

Al core 0.0 0.0 0.35203

O core 0.3064 0.0 0.25

ends

Finally, we need to add some keywords to tell the program what type of calculation we want to do and to define the force field we want to use. At this stage we only want to generate the crystal structure, so all you need to add to the bottom of your file is:

Nopotential #No force field (or potential parameters) will be defined

Check #Just read in the input and write it back out again, inc xyz and cif files

Stop #Tells program not to go further (to calc surfaces etc.)

Start #We have finished giving input to the program so it can do any calculation requested

Stop #Gives us control of the computer back when the calculation has finished

Save the file, and rename it al2o3.met (hint: remember to disable ‘Hide extensions’ for known files)

On the University computers you cannot create your own batch file system (whereby you can simply double-click on the file and make it run) – if you want to do this on your own computer – do ask.
 
To run on the University computers: 
(a) Copy the file al2o3.met and call it input.txt
(b) Copy the executable metadise.exe into the same folder, and run by double-clicking on it.
(c) Once the simulation has finished and the files that you want are created (see below), you can then delete metadise.exe (to save space).
 
(or whichever directory you’re performing the calculations in)

The program will only take a second or so to run and will produce several output files. The code*.out file is the main output file from METADISE and sum*.out is a summary of this file. Both can be opened in text editors WordPad or Notepad. The remaining files display the crystal structure in formats that can be read by various visualization packages and by METADISE itself (in the case of the .res files). The .xyz and .cif are in a format that can be read by ACCELRYS_DISCOVERY_STUDIO or VESTA.

Experiment with using the viewer to manipulate the cell.

For example you might try to change the display type or grow the crystal.

3.2 Choosing a potential model

Being able to generate a structural model of our crystal is only the first step of the calculation. We also need to choose a suitable potential model to describe the interactions in our system and calculate key properties such as energies, elastic constants and electrostatic constants. We are now going to select the potential parameters required to model Al₂O₃. At this point you might want to refer back to the briefing notes on force fields and potential models if you are unclear about any terms used below.

Create a new directory Step2.

In Alumina/Step2, copy the input file you used to create the alpha-alumina structure in step 1 and change the keyword

nopotential to @include poten.txt. This tells the program to look in an external file (poten.txt) for the interatomic potentials; the next task will be creating this file in your working directory. You want to optimize your model with this potential and need to change a few keywords in your input file. Change the check keyword so it reads conp to request that METADISE performs an energy minimization calculation to constant pressure. Also ensure that the first stop keyword has been removed so that the end of your file reads:

@include poten.txt

conp

start

stop

Save your file as al2o3_pot.met

You will then create a file poten.txt. First add the long-range Coulombic interactions. This simply requires details of the atomic species in the cell and their charge. The keyword Buckingham (buck) is used to indicate how the short-range interactions will be defined. In our model we will assume the Al—Al interaction is totally Coulombic (a common assumption in these models) so we only need to add details of the Al—O and O—O interactions. The potential file should look like this:

potential # keyword to indicate the potential model will follow

species # keyword to indicate a list of ionic species and charges will follow

Al 3.0 # Ion name and charge: Al³⁺

O -2.0 # Ion name and charge: O^2-

ends # The end of the list of ionic species

BUCK AL CORE O CORE # buck keyword to indicate a list of potential parameters will follow

1460.300 0.2991 0.0 # Potential parameters for Al-O interaction

ends # The end of this buck parameter

BUCK O CORE O CORE

22764.3 0.1490 27.88 # Potential parameters for O-O interaction

ends # The end of this buck parameter

ends # The end of the potential model definition

Now, once you have added these parameters to your potential model, run METADISE by same method as in section 3.1. Note that you may need to edit RunMetadise.bat to output in your current directory.

Once METADISE has completed open the code*.out file and scroll to the bottom of the file. If the calculation was successful you will find the phrase valid minimization to constant pressure a few screens up from the bottom. Try searching for text string “valid” to find it.

Below this you will find the lattice vectors of the relaxed system, the lattice energy and various bulk properties. These serve as a useful check on how well our chosen model reproduces the true behavior of the system we are studying. Make a note of the lattice energy listed below that phrase: you will need this presently. The fin*.res file contains the coordinates of the relaxed system in METADISE input file format, useful as the starting point for additional calculations. Of course, you can also view the relaxed structure in ACCELRYS_DISCOVERY_STUDIO. The af*.car file contains the structure after minimization and the bf*.car the initial structure we defined in the input file. Compare the two to see if there has been an appreciable change in the structure.

3.3 Create surfaces

We will generate (hkl) surfaces of a-alumina from the optimized bulk obtained in section 3.2 with the potential. We want to calculate the attachment and surface energies which can be used to predict the morphology of the crystal in using a set of different surfaces. {Note: as we are using the hexagonal representation, the Miller indices have 4 digits (abcd), but as the 3^rd c=-a-b is unnecessary we do not include it, but in the text replace it by “.”}. Let´s consider the case on the (01.2) surface as follows:

Step3_hkl: For the case on the (01.2) surface, copy the fin_xxx.res in a new input file in Step3_012 directory, surf_012.met

After the inclusion of the poten.txt, you will transform the optimization keywords for the bulk as follows:

Check

surf # to generate a surface - you can replace by slab

miller 0 1 2 # for the case of the (01.2) surface

start

stackgen systematic # all code generated

start

stop

stop

Run METADISE. This step will prepare the file to generate the surface to relax.

In the pldro*.out you will see a schematic of the surface unit cell for the 4 most likely terminations. From the summ*.out file you will see the energies, and that code 13 (this might change depending on the order of coordinates) has both the lowest attachment energy -0.54130 (eV/ion) and the lowest unrelaxed surface energy 3.61701 (J/m^2).

Make a sub-directory, 13, copy surf_012.met (and poten.txt) into it and replace

stackgen systematic

by

stackgen 13

Then rerun METADISE. Then copy the stac_xxx.out as a new input file, code13.met and run METADISE again, this will optimise the surface geometry.

Note: METADISE had changed the axis on the bulk. The surfaces are perpendicular to the x direction (the stacking vector). The depth of the surface is specified by the region command in the stac file and it specifies the number of unit cells. If using this to generate coordinates for DFT, the first two numbers are normally reduced.

The details of the minimization are in the code*.out file, but again the summary can be found in summ*.out.

Counter CODE unrel_S_en/Jm^-2 Suface_en/Jm^-2

0 13 3.62645 2.82056

The unrelaxed and relaxed surface energies are given at the bottom of the file. These are calculated from the bulk (which has twice the number of ions as for the surface), the surface energy and area.

g = 16.021 ( Es – Eb/2)/A

Where 16.021 is the conversion from eV/A^2 to J/m^2, A=area of surface unit cell, Es is the optimised surface energy (on newt line) and Eb is the bulk (on grad line)

Es=-5786.507450 eV, Eb=-11598.341280 eV and A= 71.927821 A^2

This can now be duplicated for other terminations of this surface and the other low index surfaces.

fin_o#####.res contains the relaxed structure, and will be used in 3.6. Note final structure, and surface coordinations, obtained from viewing the structure

Another important surface is the (00.1) surface. We are now going to repeat section 3.3, but changing miller 0 1 2 to miller 0 0 1. However, the crystallographic unit cell is too small and only one surface termination is exposed. If you try it you will get:

height charge atom arrangement

0.0 3.0 AL--

0.2 0.0

0.4 0.0

0.6 0.0

0.8 -6.0 O--O--O--

1.0 0.0

1.2 0.0

1.4 3.0 AL--

Therefore, the unit cell is grown so that each oxygen layer will have 6 oxygen atoms, which can be conveniently cut into two.

Step3_hkl: For the case on the (00.1) surface, copy the fin_xxx.res (from 3.2) to make a new input file (naming it surf_001.met), and create a Step3_001 directory to put it in.

After the inclusion of the poten.txt, transform the optimization keywords for the bulk as follows:

Check

grow 1 2 1 # we want to grow one of the direction

surf # to generate a surface

miller 0 0 1 # for the case of the (0001) surface

start

stackgen systematic # all code generated

start

stop

stop

Run METADISE. We want to optimize the surface with the oxygen termination which is the code 6. If you look in the pldro*.out, we will see the type of termination in each code. The code 1 is terminated by Al atoms and the code 6 by Oxygen atoms. The attachment energies given in the sum*.out indicate that there are 2 types of different cut (code 1 and 6).

Note 1: the Al termination is in fact the most stable.

Note 2: If the bulk is without the growth, the O termination is not generated. The growth allows more configurations.

Make a sub-directory, 6, copy surf_001.met (and poten.txt) into it and replace

stackgen systematic

by

stackgen 6

Then rerun METADISE and copy the stac_xxx.out as a new input file, naming it code6.met, and use this to run METADISE again to optimise the surface geometry.

The details of the minimization are, as before, found in the code*.out file and summarized in the summ*.out. Verify that the calculation optimized correctly (code is valid in the sum*.out). Verify the structure of the first code (i.e. aft_o0001.cif) with DISCOVERY STUDIO.

Note: in the code13.met file, METADISE has chosen 3 of the oxygen atoms to remain at the surface and the other 3 transferred below to remove the dipole. If you wish to try different permutations, add trans 1 at the end of the line containing the first oxygen and trans -1 at the end of the line of the final oxygen. This then swaps them and maintains a zero dipole.

You should now be able to generate a table of the format shown below using the information found in sum_xxx.out:

surface	code	Surface area	Attachment Energy	Surf Energy Before relax (J/m^2)	Surf Energy After relax (J/m^2)
0001	1	xxx	xxx	xxx	xxx
x	x	xxx	xxx	xxx	xxx

Try other possible surfaces such as (10.1), (11.0), (01.2), (11.3) and (22.3).

Note1: These other surfaces should have a bigger surface area and therefore you do not need to grow the slab (i.e. delete the grow 1 2 1 command at the beginning of this section)

Note2: The higher the Miller indices the longer the calculations will take.

3.4 Generate morphologies

Step4: we will use the Gibbs-Curie-Wulff construction to simulate the morphology of the crystal under thermodynamic equilibrium. Basically, the morphology depends on the relative surface energies of each surface and the symmetry of the bulk.

After collecting all the surface energies calculated, you will use the bulk structure and complete the file with the calculated values. If some surfaces have not been calculated, just remove them or add a comment symbol in front of them (e.g. for the surface -2 1 6)

print pov 1 # to generate povray files

cell 4.7589 4.7589 12.9910 90 90 120

space full R3CH 167 1

fractional

Al core 0.0 0.0 0.35203

O core 0.3064 0.0 0.25

ends

nopoten

check

nano

start

gem

space full R3CH 167 1

index 0 0 1 xxx # replace xxx by the surface energy that you have

#obtain in section 3.3 step3_001

index 1 0 1 xxx

index 0 1 2 xxx

index 1 1 3 xxx

index 2 2 3 xxx

index 1 0 2 xxx

# index -2 1 6 xxx

ends

start

stop

Run METADISE.

Note: You do not need to have a potential, or to reoptimize the bulk.

Generate files to see the morphology such as .wrl to open with Flux Studio.

The codexx.out gives the relative surface area of each surface.

3.5 Construct nanoparticles

Constructing nanoparticles uses the same idea as constructing the morphology of particles (in the previous section 3.4). You will choose the surfaces which have been calculated and that you want to expose (stable surfaces as a start). So let´s start with 2-3 surfaces.

Two steps are needed. The first is to localize the particle without dipole moment and then focus on the center satisfying that property. Such as:

1/ Step5: Run METADISE with this example (nanoparticle of 25 A of radius with xxx surfaces):

PRIN pdb 0 car 0 xtl 0 xyz 0 rel 0

#print dlpoly 1 reg2 2

DIME 50000000

THERMAL

CUTO 1.000 15.0 0.5 60.0 60.0 1.0

cell 4.7589 4.7589 12.9910 90 90 120

space full R3CH 167 1

fractional

Al core 0.0 0.0 0.35203

O core 0.3064 0.0 0.25

O shel 0.3064 0.0 0.25

ends

@include poten.txt

#conp

print dlpoly 1 reg2 2 # to print a dl_poly input file

check

nano

start

gem

space full R3CH 167 1

index 0 0 1 2.42

#index 0 1 2 xxx

#index 1 1 3 xxx

index 1 1 0 2.67

index 1 0 0 2.86

#index 1 0 1 xxx

#index 2 2 3 xxx

nano

rad1 25. # radius of the particle

fixcharge # finds a nanoparticle with zero charge

#cent 0.0 0.0 0.0 # the default origin is the crystallographic

# but can be changed

ends

start

stop

The nanoparticles can be minimized by simply adding, after start

minimise

conj 500

start

However, it is far more efficient to anneal in an MD calculation. To do this follow the instructions in Optional Material and adapt the command for the nanoparticle (for example, ensure that the simulation cell size in CONFIG, is sufficiently large, e.g. a=b=c=100 Angstrom).

3.6 Adding Impurities

The surfaces can be a sink for impurities and hence we can calculate the energy of an impurity as a function of coverage and of depth. Select a surface and replace the lattice ion with a cation impurity (use viewerlite to identify surface and bulk ions), this may also in turn have to be charge compensated. This is most easily done by removing such units in the fin######.res file generated in 3.3 or in 3.5 (the final one generated from the optimized surface geometry) and renaming this defect.met and running it with METADISE, ensuring the keywords from “bulk” to “meta” have been removed. The potential models for a number of cation impurities are listed below in the appendices. The energies for impurities are found on the “newt valid” line of the sum_xxx.out file.

a. Calculate the energy of an isolated impurity in the bulk. The defect must not be too close to region2 (this is the fixed region) or indeed too close to the surface (the first part of region 1 in the input file). If the energy of the simulation cell without the defect is E_bp and the energy of the cell with the defect is E_bd then the defect energy, E_def=E_bd-E_bp

b. Calculate the energy of the same operation, but at the surface to give E_defs. The segregation energy E_seg=E_defs-E_def.

c. Calculate this quantity as a function of depth. Is there a simple relationship with depth?

d. Consider increasing concentration of impurity by scaling the surface area and hence the number of Al ions at the surface. You need to scale the simulation cell to have a larger surface area. The segregation energy becomes: E_seg=(E_defs-n.E_def)/n, where n is the number of defects and E_defs is the energy difference due to a surface containing n defects.

e. The new surface energy g_def = g + 16.021 (E_seg )n/A Where n/A is the number of defects per unit area (and the area may be scaled).

The new lowest surface energy for each Miller Index can be used to generate a new morphology.

Although perhaps beyond the scope of this exercise, the energies of segregation can be calculated to the nanoparticles themselves to investigate the adhesion to edges and kinks as well as flat surfaces.

Step6_hkl: Copy the fin*.res from 3.3 and put in new folder (as before)

Add at the end:

@include poten.txt

surface

xscale 2 2 # to scale the surflat vectors so that area is 4 times bigger

start

STOP

Optional Material

Add hydroxyl groups to surfaces

Many oxide surfaces are terminated by hydroxide groups when exposed to the atmosphere, which may affect the surface stability. There is a function in METADISE which enables you to add hydroxyl groups on the surface. The steps are described as follows:

Step7_hkl: Create a new directory Step7_001 for the adsorption of OH on (00.1) surface

Copy the input file from section 3.3 code6.met to a new file called 0001_OH_code6.met

Change the “surf” keyword to “slab”

Add hydroxyl after the “slab” keyword:

adsorb

hydrox Al core o core 2.2

ends

See Appendix 3 for the new poten.txt file containing the OH groups.

Run METADISE where you should have the poten.txt in your working directory.

Verify the structure before relaxation (bef_o0001.car or cif). Then you should optimise this structure as before.

The METADISE energies can be converted into surface energies, but unlike an ab initio calculation, the oxygen ions are different and a Born Haber cycle has to be applied, see for example De Leeuw J. Am. Ceram. Soc., 82 [11] 3209–16 (1999). In this case the energy is estimated to be 10.58 eV per water molecule. So that if 3 water molecules are added to surface with area 39.604163 A^2, the surface energy is corrected by 3*10.58*16.021/39.604163

Run METADISE.

At this stage, you have been able to generate a surface with hydroxyl groups. This is a structure that you will use to run a MD calculation with DL_POLY in adding water on the top of the surface.

Molecular Dynamics

DL_Poly files can be generated by adding

Print dlpoly 1

Once the program has run you should get the same output files as we generated in section 2, but in addition you should have three files with the suffix *.dlp

Open the file CONFIG*.dlp. This file contains information about the system’s configuration and is equivalent to the LATTICE and BASIS section of a METADISE input file. Close the file and rename it CONFIG (all capitals with no suffix – ignore any warnings about changing suffix type).

Now open the file FIELD*.dlp. This contains information about the force field and is analogous to the potential section of a METADISE file. Hence, the order in which the atoms appear is important, and must be exactly the same order as in the CONFIG file. Because we are only using a simple pair-wise potential model the FIELD file is quite simple, however this is not always the case. Close the file and rename it FIELD (again, all capitals with no suffix – ignore any warnings about changing suffix type).

Open the remaining file, CONTROL*.dlp. This file, as its name suggests, contains all the control variables for the simulation, i.e. it tells the program what to do. METADISE generates a template file with some standard values for a typical simulation; however for the simulation we are going to perform we will need to change a few of these values.

1.         Check that the time step is set at 0.001 ps (1 fs)

2.         Change the number of ‘steps’ to 20000

3.         Select an npt constant pressure ensemble by placing # in front of the line ensemble nve and removing the # from in front of the line below and changing the keyword nst to npt. The # is the comment character and means that the program will not read the remainder of the line. Running the NPT ensemble enables the size of the simulation cell to change during the MD run.

4.         Change the values nstraj= 1 istraj= 250 to nstraj= 0 istraj= 100. This changes how often the program writes out to the HISTORY file (more on this later)

Once you have made these changes save the file as CONTROL. (again, all capitals with no suffix – ignore any warnings about changing suffix type). See CONTROL in NVT in appendix 4.

NOTE: The reliability of the result will depend on the number of steps as this improves the statistics. Thus, if the computer is fast enough, or you are leaving it running etc, try increasing the number of steps, but be careful or you may spend too much time waiting.

Run MD with DL_POLY

All DL_POLY simulations should be run in separate folders. You need CONFIG, CONTROL and FIELD.

Step7: Run an NVE calculation with DL_POLY

Run an NVE calculation until the stabilization of the energy (see STATIS, third column).

Once the job has reached the given time, you should find a number of new output files in the working directory: HISTORY, OUTPUT, RDFDAT, REVCON, REVIVE and STATIS.

Step8: Run an NVT calculation for you optimized calculation (copy the REVCON in a new CONFIG file) for temperature at T=300K then choose another temperature between 400-1000K. See CONTROL in NVT in Appendix 4.

A proper optimization should take hours and hours. During this initiation, we will limit the running calculation to 15mn, the objective being to show how to run a calculation and how to analysis the output. If you find that your DL_POLY job is terminating immediately, please inform a demonstrator or/and you can look in DL_POLY manual to identify the type of error (in looking at the corresponding number of errors at the end of the manual).

In VMD:

In File menu, open New molecules, load HISTORY with the type file to load DLPOLY v2 HISTORY

In the Display menu, click on orthographic.

You can run a movie at the bottom of the VMD main menu.

Acknowledgement:

Thanks to Dr Corinne Arrouvel (Rio)

Appendix 1: Potential Model for Al₂O₃.

Potential

spec

Al 3.

O -2.

ends

buck Al core o core

1460.300 0.2991 0.0

ends

buck o core o core

22764.3 0.1490 27.88

ends

ends

Appendix 2a: Potential Model for cation doping of Al₂O₃

Potential

SPEC

O CORE -2.0

Y CORE 3.00

GD CORE 3.00

EU CORE 3.00

HO CORE 3.00

LA CORE 3.00

LU CORE 3.00

ND CORE 3.00

TB CORE 3.00

MG CORE 2.0

CR CORE 3.0

FE CORE 3.0

AL CORE 3.0

ENDS

BUCK

AL CORE O CORE 1460.3 0.2991 0.0

FE CORE O CORE 1102.4 0.3299 0.0

CR CORE O CORE 1734.1 0.3010 0.0

Y CORE O CORE 1345.1 0.3491 0.0

EU CORE O CORE 1358.0 0.3556 0.0

GD CORE O CORE 1336.8 0.3551 0.0

HO CORE O CORE 1350.2 0.3487 0.0

LA CORE O CORE 1439.7 0.3651 0.0

LU CORE O CORE 1347.1 0.3430 0.0

ND CORE O CORE 1379.9 0.3601 0.0

TB CORE O CORE 1369.7 0.3521 0.0

O CORE O CORE 22764.3 0.149 27.88

MG CORE O CORE 1428.5 0.29453 0.00

ENDS

ENDS

Appendix 2b: Potential Model for cation doping of Al₂O₃

poten

#

# LEWIS library - collection of potentials based

# around the Catlow oxygen-oxygen potential#

# Reference:

# (1) G.V. Lewis and C.R.A. Catlow, J. Phys. C: Solid

# State Phys., 18, 1149-1161 (1985)

#

species

# Table 1 species

Ca core 2.00000

Cr2 core 2.00000

Mn2 core 2.00000

Fe2 core 2.00000

Co core 2.00000

Ni core 2.00000

Zn core 2.00000

Zr core 4.00000

Cd core 2.00000

Hf core 4.00000

Ce core 4.00000

Eu2 core 2.00000

Tb core 4.00000

Th core 4.00000

U core 4.00000

Sc core 3.00000

Mn core 3.00000

Y core 3.00000

La core 3.00000

Nd core 3.00000

Eu core 3.00000

Gd core 3.00000

Ho core 3.00000

Yb core 3.00000

Lu core 3.00000

Pu core 3.00000

# Oxygen

O core -2.0

ENDS

BUCK

Ca core O core 1227.7 0.33720 0.00000 0.0 10.0

Cr2 core O core 619.8 0.33720 0.00000 0.0 10.0

Mn2 core O core 832.7 0.33720 0.00000 0.0 10.0

Fe2 core O core 725.7 0.33720 0.00000 0.0 10.0

Co core O core 684.9 0.33720 0.00000 0.0 10.0

Ni core O core 641.2 0.33720 0.00000 0.0 10.0

Zn core O core 700.3 0.33720 0.00000 0.0 10.0

Zr core O core 1453.8 0.35000 0.00000 0.0 10.0

Cd core O core 868.3 0.35000 0.00000 0.0 10.0

Hf core O core 1454.6 0.35000 0.00000 0.0 10.0

Ce core O core 1017.4 0.39490 0.00000 0.0 10.0

Eu2 core O core 665.2 0.39490 0.00000 0.0 10.0

Tb core O core 905.3 0.39490 0.00000 0.0 10.0

Th core O core 1144.6 0.39490 0.00000 0.0 10.0

U core O core 1055.0 0.39490 0.00000 0.0 10.0

# Table 1 potentials

Sc core O core 1299.4 0.33120 0.00000 0.0 10.0

Mn core O core 1257.9 0.32140 0.00000 0.0 10.0

Y core O core 1345.1 0.34910 0.00000 0.0 10.0

La core O core 1439.7 0.36510 0.00000 0.0 10.0

Nd core O core 1379.9 0.36010 0.00000 0.0 10.0

Eu core O core 1358.0 0.35560 0.00000 0.0 10.0

Gd core O core 1336.8 0.35510 0.00000 0.0 10.0

Ho core O core 1350.2 0.34870 0.00000 0.0 10.0

Yb core O core 1309.6 0.34620 0.00000 0.0 10.0

Lu core O core 1347.1 0.34300 0.00000 0.0 10.0

Pu core O core 1376.2 0.35930 0.00000 0.0 10.0

# Oxygen-oxygen potential

O core O core 22764.0 0.14900 27.87900 0.0 12.0

ends

ends

Appendix 3: Potential Model (2a) type for hydroxylated Al₂O₃ with water.

Potential

spec

Al 3.

O -2.

OH -1.4 16.0

OW -0.8 16.0

H 0.4 1.008

HW 0.4 1.008

ends

buck

Al core o core 1460.300 0.2991 0.0

AL CORE OH core 1022.300 0.2991 0.0

AL CORE OW core 584.300 0.2991 0.0

o core o core 22764.3 0.1490 27.88

O core OH core 22764.0 0.149 13.94

O core OW core 22764.0 0.149 28.92

OH core OH core 22764.0 0.149 6.97

OH core OW core 22764.0 0.149 17.14

H CORE O core 396.27 0.25 0.0 20.0

H CORE OH core 311.97 0.25 0.0 1.3 20.0

HW CORE O core 396.27 0.25 0.0 20.0

HW CORE OH core 311.97 0.25 0.0 20.0

ENDS

spring

OH core H core 13.0 0.9 1.0 0.0 1.3

ends

#TIP3P

SPRING HW CORE OW CORE

23.9907488 0.9572 1.0 0.0 1.2

ENDS

LENN OW CORE OW CORE 12 6

25246.059 25.80524056 0.0 20.0

ENDS

BOHA 1 4.3382909 104.52

ends

thbo

ANGA 1 OW CORE HW CORE HW CORE

1.2 1.2 1.65

ends

Appendix 4: input with adsorb hydroxyl group (from inclusion potential).

@include poten.txt

#bulk

#start

#minimise

#conj 1

#bfgs 0

#newt 0

#maxu 10

#nokill

#start

#meta

slab

#

# add hydrox

adsorb

hydrox al core o core 1.9

ends

start

minimise

conj 0

bfgs 20

newt 600

maxu 20

nokill

start

anal

start

stop

Appendix 5: DLPOLY. CONTROL file for NVT.

# for Al2O3

surfa___

temperature 300.00000

pressure 0.0

steps 1000000

equilibration 50000

timestep 0.000200000 ps

#restart

scale 1

cutoff 8.0 angstrom

#primary cutoff 1.73 angstrom

delr 0.8 angstrom

print 250

rdf 250

print rdf

job time 6000 seconds

close time 100 seconds

ewald precision 1d-5

cap forces 1000.0

#ensemble nve

ensemble nvt hoover 0.5

stats 250

#trajectory nstraj= 1 istraj= 250 keytrj=0

trajectory 1 250 0

stack 100

shake tolerance 1.d-8

finish

Appendix 6: DLPOLY. Example of FIELD file with water.

units eV

molecules 5

Alfix

nummols 1

atoms 1

AL 26.000 3.000 0 1

finish

Aluminium

nummols 719

atoms 1

AL 26.000 3.000

finish

OXYGEN

nummols 990

atoms 1

O 15.499 -2.0

finish

OH CORE

NUMMOL 0

atoms 2

OH 16.0000 -1.000 0 0 0

H 1.0080 0.400 0 0 0

bonds 1

harm 1 2 23.9907488 0.9000

finish

TIP3P WATER

nummols 0

atoms 3

OW 16.0000 -0.834

HW 1.0080 0.417

HW 1.0080 0.417

bonds 2

harm 1 2 23.9907488 0.9572

harm 1 3 23.9907488 0.9572

angles 1

harm 2 1 3 4.3382909 104.52

finish

vdw 14

AL O buck 1460.3000 0.2991 0.0000

AL OH buck 1022.3000 0.2991 0.0000

AL OW buck 584.3000 0.2991 0.0000

O O buck 22764.0000 0.1490 27.8800

O OH buck 22764.0000 0.1490 13.9400

O OW buck 22764.0000 0.1490 28.9200

O H buck 396.2700 0.2500 0.0000

O HW buck 396.2700 0.2500 0.0000

OH OH buck 22764.0000 0.1490 6.9700

OH OW buck 22764.0000 0.1490 17.1400

OH H buck 311.9700 0.2500 0.0000

OH HW buck 311.9700 0.2500 0.0000

OW OW 12-6 25246.0590 25.8052

OW H buck 396.2700 0.2500 0.0000

CLOSE

Appendix 7: Other Structures.

#Chromia

space full R3c 167 1

CELL 4.950700 4.950700 13.565600 90.000000 90.000000 120.000000

fractional

Cr core 0.000000 0.000000 0.347660

O core 0.305100 0.000000 0.250000

ends

# hematite

CELL 5.03594 5.03594 13.74439 90.000000 90.000000 120.000000

space r-3c 167 1

frac

Fe core 0.000000 0.000000 0.35528

O core 0.307100 0.000000 0.250000

ends