Surface X-diffraction tutorial

At the present GenX implements two different facilities to handle sxrd data.

  • SXRD data loader

  • SXRD model to calculate sxrd data

Loading SXRD data

First the data loader has to be changed to load files for sxrd. Open the menu Settings‣Data Loaders (or press Shift-Ctrl-D) and an dialog pops up. Choose the sls_sxrd data loader in the drop-down box that appear. The data is loaded by marking the data set where the rods should be inserted and choose File-Import-Import Data (It is also possible to click on the small open icon above the data sets). A dialog appears and the data file to load is chosen. An example for an valid data file can be found in rods.dat or see below.

# Test file created: 2009-04-14
# Headers:
#h        k         l         I         Ie
0.000e+00 0.000e+00 1.000e-01 0.000e+00 0.00e+00
0.000e+00 0.000e+00 1.100e-01 0.000e+00 0.00e+00
0.000e+00 0.000e+00 1.200e-01 0.000e+00 0.00e+00
0.000e+00 0.000e+00 1.300e-01 0.000e+00 0.00e+00
0.000e+00 0.000e+00 1.400e-01 0.000e+00 0.00e+00
0.000e+00 0.000e+00 1.500e-01 0.000e+00 0.00e+00
0.000e+00 0.000e+00 1.600e-01 0.000e+00 0.00e+00

The data format that this data loader wants is a 5-column consisting of h, k, l, Intensity and Intensity error separated with a whitespace (comma, tabs..). The default token for a comment is a #. The settings for the data loaders are accessed through Setting-Import Settings (or Shift-Ctrl-I). If you would happen to need to create your own test files you can use the python script make_sxrd_data.py. .. image::_attachments/surfdiff/data_loaded.png

When the data is loaded new data sets are appended automatically. Each data set is given the name of the type ‘(h, k)’. If you started from the beginning you should also remove the data set ‘Data0’ since the data sets are appended to the previous sets.

Creating an SXRD model

Okay this is the tricky part first off we will create a simple model consisting of and 2 unit cell thick L aAlO3 film on top of a SrTiO3 substrate. The script will look something like this:

# 1
import models.sxrd as model

# 2 Defining the unit cell parameters
unitcell = model.UnitCell(3.945, 3.945, 3.945, 90, 90, 90)
# 3 Define the instrument
inst = model.Instrument(wavel = 1.0, alpha = 1.0)

# 4 Defining the bulk
bulk = model.Slab()
# 4.a Define the atoms
bulk.add_atom('Sr', 'sr', 0.0, 0.0, 0.0, 0.08, 1.0)
bulk.add_atom('Ti', 'ti', 0.5, 0.5, 0.5, 0.08, 1.0)
bulk.add_atom('O1', 'o', 0.5, 0.0, 0.5, 0.08, 1.0)
bulk.add_atom('O2', 'o', 0.0, 0.5, 0.5, 0.08, 1.0)
bulk.add_atom('O3', 'o', 0.5, 0.5, 0.0, 0.08, 1.0)

# 5 Creating an LaAlO3 unit cell - note asymmetric unit cell
laouc =  model.Slab(c = 1.05)
# 5.a Define the atoms
laouc.add_atom('La', 'la', 0.0, 0.0, 0.0, 0.08, 1.0, 1)
laouc.add_atom('Al', 'al', 0.5, 0.5, 0.5, 0.08, 1.0, 1)
laouc.add_atom('O1', 'o', 0.5, 0.5, 0.0, 0.08, 1.0, 1)
laouc.add_atom('O2', 'o', 0.5, 0.0, 0.5, 0.08, 1.0, 2)

# 6 Symmetry operations
p4 = [model.SymTrans([[1, 0],[0, 1]]), model.SymTrans([[-1, 0],[0, -1]]),
           model.SymTrans([[0, -1],[1, 0]]), model.SymTrans([[0, 1],[-1, 0]]) ]

# 7 Creating the sample with one unit cell
sample = model.Sample(inst, bulk, [laouc], unitcell)

#8 Put the right symmetry in the surface
sample.set_surface_sym(p4)

# 9 Define the Sim function
def Sim(data):
    I = []
    #9.a loop through the data sets
    for data_set in data:
        # 9.b create all the h,k,l values for the rod (data_set)
        h = data_set.extra_data['h']
        k = data_set.extra_data['k']
        l = data_set.x
        # 9.c. Calculate the struct factor
        f = sample.calc_f(h, k, l)
        # 9.d Calculate the intensity
        i = abs(f)**2
        # 9.e Append the calculated intensity to the list I
        I.append(i)
    return I

So to create a model script the following actions has to be made:

  1. The first row imports the necessary model import models.sxrd as model.

  2. Create a UnitCell(a, b, c, alpha, beta, gamma) this object contains all the lattice parameters of the substrate. This also relates the (h, k, l) values in the measurement to the real world.

  3. Create an Instrument(wavel, alpha). This contains all the parameters related to the measurement. Wavelength in Angstroms and incidence angle in degrees.

  4. So now we define the bulk Slab, this takes the instrument as input.

    1. Add all the atoms to the bulk (here we use a P1 symmetry), the command has this syntax bulk.add_atom(id, element,x, y, z, u, occ). id denotes a unique id for the atom added. element is the element in question which has to defined in the f and rho table. x, y, z is the base position of the atom (movements of atoms is made through the variables dx, dy, dz}) . u is the Debye-Waller parameter (isotropic) and occ is the occupancy of the atom.

  5. Next up is the same procedure for the film atoms. Here we use the same class Slab but we create it with a twist. As an optional arguments we insert c = 1.05. This will make the c-axis constant of this Slab 5% longer than the one defined in the unitcell.

    1. Now we add all the symmetry in-equivalent atoms to this Slab, laouc, assuming that we have in plane p4 symmetry (more about the symmetry later on). Also there is then an extra argument to the add_atom method. The call is now on the form: laouv.add_atom(id, element, x, y, z, u, occ, mult). The optional parameter mult is then the multiplicity of the site, i.e. how many atoms there are of this atomic site in the unitcell (it is tabulated in the international tables of crystallography <http://it.iucr.org/>). For a non-special position (general) this will be 4. For the (0.5,0.5) position this will then be 1 as it will occur only one time.

  6. Now both the bulk unit cell and the surface is defined as slabs. What we have to do now is to define the symmetry operations for the top layer. This is done by creating SymTrans objects and placing them in a list so we build up all allowed symmetry operations of the plane group. The call is: SymTrans(P = [[1, 0],[0, 1], t = [0, 0]) where also the default values are shown. P is the rotation matrix of the transformation and t is the displacement vector. The default value is creating the identity matrix for P and the no displacement at all for t.

  7. A sample is creating by combining the instrument, bulk slab, surface slabs and the unitcell. Note that several surface slabs can be stacked on top of each other to create more complex structures. Call: Sample(inst, bulk_slab, surface_slabs, unit_cell) .

  8. To be added

  9. Now we need to define the Sim function this is a mandatory function that has to included in all model scripts. It takes a data structure as input. First we create an empty list which will contain the calculated data. This is the function that will be executed when you press on the simulate and when you fit your data.

    1. Then we loop through all data sets that we have loaded.

    2. For each data_set we extract the h,k and l values from the data file. Note that h,k is loaded as extra data but the l values which is the independent variable in the plots is retrieved by data_set.x. (As a side note data.extra_data is an dictionary of data which is created by the sls_sxrd data loader)

    3. Then we calculate the _complex_ structure factor for the sample. The syntax should be obvious.

    4. The intensity the square of the absolute of the structure factor is calculated and is appended to our result list I

    5. The result list (which has to have exactly the same length as the number of data sets), I, is returned from the Sim function.

Making a simulation

To simulate the data click on the yellow lightning in the main tool bar. The following screen should appear.

![Screenshot - simulation, linear scale](/apps/trac/genx/raw-attachment/wiki/IntroGuides/SurfDiff/simulation_lin.png)

In order to see anything the scaling has to change from linear to logarithmic. Right click on the plot and choose y-scale‣log.

This should make the plot a bit clearer. However, we still have a bunch of lines on top of each other.

Lets make som of the data sets non-active. Mark all data sets in the right data list except the first one, right click and choose Toggle show. This turn on and off the plotting of the selected data sets. Toggle Active sets a flag whether or not to the data should be simulated. Naturally this flag must be read and used in the simulation function to have an effect.

Adding more layers

Usually I am more interesting in fitting and simulating films that are a couple of unitcells thick. So here a couple of tricks to make that that easy. If you just want to add identical unitcells on top of each other you can just expand the list, to simulate lets say 3 UC:

# 7 Creating the sample with one unit cell
sample = model.Sample(inst, bulk, [laouc, laouc, laouc], unitcell)

Another possiblity is to use the smart list functionality that python offers, which makes the following code equivalent to the top one:

# 7 Creating the sample with one unit cell
sample = model.Sample(inst, bulk, 3*[laouc], unitcell)

If you have several non-identical layers lets say another Slab called stouc we can make an 2 UC over layer on top of the loauc.

# 7 Creating the sample with one unit cell
sample = model.Sample(inst, bulk, [laouc, laouc, laouc, stouc, stouc], unitcell)

which can also be written as:

# 7 Creating the sample with one unit cell
sample = model.Sample(inst, bulk, 3*[laouc] + 2*[stouc], unitcell)

An important consequence of reusing the defined unitcells is that all the parameters in the repeated unitcells is the same. For example if you would like to fit identical displacement parameters for several layers this is the way to go. Eventually you might want to decouple the different layers and fit them separately. Then you have to make a copy each slab into a separate object and place into a list. The following code copies one slab into a new one with another name.

sto_other = stouc.copy()

The object sto_other will then contain atoms at the same position as stouc but when a property of stouc is changed it will not affect the copy, sto_other. Note that if you would write

sto_other = stouc

a change in stouc or sto_other would affect the other! Naturally, in this object has also to be included in the sample in order to be simulated.

# 7 Creating the sample with one unit cell
sample = model.Sample(inst, bulk, 3*[laouc] + [sto_other, stouc], unitcell)

Defining parameters to fit

Before we proceed with the more advanced topics we will briefly review how to define which parameters to fit. As a user you should be aware that in GenX _everything_ can be fitted. However, it is no guarantee that the parameters that you have chosen to fit make any sense at all! So you should use your judgment and think about the physical model before you decide to fit (or freeze) a parameter. Fitting is also a little bit like craftsmanship, so in order to be good at it you have to practice in using the tools.

The definition of which parameters the program fits is done in the tab Grid (lower panel). There are 6 columns. The first, Parameter,is the name of the _function_ that sets the parameter, next, Value}}, the default/refined value to set the parameter to. Column number three, {{{Fit selects if the parameter in this row should be refined. Next the boundaries for the refinement is given by the Min and Max columns. The last column Errors display the calculated result of the error calculation that can be conducted after a fit.

No to find the name of the parameter right click on the a box in the parameter column (Note the box is not allowed to be in edit mode with a cursor blinking). This will cause a pop-up menu to appear with all possible set functions for the objects you have defined in the script. All function that appear in this menu and as parameters will appear on the form object_name.set_parameter for example to fit the slab-global c-axis of the stouc

stouc.set_c

would be chosen. Most parameter will also define its current value in the Value column and define a +/- 25% span between the Min and Max columns. So, care must be taken that the boundaries is physical and relevant to the parameter you have chosen.

The grid can also be used as a quick way to play around with the parameters in the model. Something to keep in made is the following:

  1. When you press simulate the script will be compiled.

  2. The Sim function will not be evaluated yet.

  3. The values in the Grid will be set to their values in the Value column.

  4. The Sim function will be evaluated.

When you build your model you think about that when fitting the following procedure will be followed:

  1. Set the parameter values.

  2. Evaluate the Sim function.

  3. Calculate the figure of merit.

So if you want to do something special to your parameters, more about that later on, you can not do it outside the Sim function as that code will only be evaluated as you simulate and recompile the model.

If you would like to try a quick fit just press the green arrow in the tool bar and of you go. Keep an eye on the status bar at the bottom of the window. This will display important information about how the it proceeds and what the program does. To stop the fit press the red stop sign in the toolbar and wait until a dialog appear which asks you if you want to keep the Values the program has fitted. If you would like to resume the fit, if you, for example, stopped it to early, press the round green arrow to the right of the stop button.

During fitting you should keep an eye on the FOM tab and the Pars tab which shows the progress of the algorithm. For a more detailed discussion about this see the x-ray tutorial :ref:_tutorials-xrr-fitting.

Scaling the simulation

When fitting and simulatiing one very important thing is to also fit/change the scale factor. One can fit the parameter inst.inten which is the “incoming intensity”. This parameter multiplies with the structure factor when using the sample.calc_f function. A more effective way of doing it when fitting data, to my opinion, is to scale each simulation separately. The sxrd model contains two function that accomplish this task: model.scale_sim(data, I) and model.scale_sqrt_sim(data, I). Both of these functions will return a new intensity array that are the least-squared fitted simulation to the data with respect to the scale factor. The minimization is done by analytically by solving \(\mathrm{min}_s \sum_i \left(D_i - sI_i\right)^2\) or \(\mathrm{min}_s \sum_i \left(\sqrt{D_i} - s\sqrt{I_i}\right)^2\). This function should be placed right before returning the value in the Sim function.

    ...
    # 9.e Append the calculated intensity to the list I
    I.append(i)
I = model.scale_sqrt_sim(data, I)
return I

Grouping atoms

A common sense approach can be that several atoms should move together and not independently of each other. For example, strain field will compress/expand an epitaxial film and consequently the lattice parameter will change. We have already briefly dealt with this as we talked about Slab`s and its parameter c which can scale the thickness of the slab. This should be the first attempt to use when fitting the data. Other, more elaborate, models might be to move atoms of the same element/site in the same manner. One obvious example of this if an alloyed site is simulated. This is where atoms group can come in handy. The class that handles this is called AtomGroup.

Basic usage

Naturally there are different ways of creating and working with these groups. There are two ways to create a group:

  1. The Slab.add_atom returns an AtomGroup object. Easy to use when only a small number of atoms has to chosen. For example:

    la_atom = laouc.add_atom('La', 'la', 0.0, 0.0, 0.0, 0.08, 1.0, 4)

  2. Fetch one atom from a Slab with its unique identity (name), id. This is done by typing the id between brackets (cmp. dictionary lookup):

    la_atom = laouc['La']

  3. Atoms can be searched by the Slab.find_atoms method. This lets you use a logical expression to locate the atoms. To locate an atom all at a site (0,0,0) one would write:

    la_atom = laouc.find_atoms('x == 0 and y == 0 and z == 0')

If we would have included two atoms at (0,0,0), perhaps interdiffusion between Sr and La, these two would now be moved together if the functions la_atom.setdx, la_atom.setdy,``la_atom.setdz`` would be called. The same is valid for the occupation oc and the Debye-Waller parameter, u.

There are more to this thing with groups as they stand now they are handy but groups can also be added together to form super groups. For example I want to have the same Debye-Waller parameter of all oxygen in both the stouc slab and the laouc slab. This could be done with:

all_ox = loauc.find_atoms('el == "o"') + stouc.find_atoms('el == "o"')

As a parameter in the grid we would choose the function all_ox.setu as a parameter.

Finally as a practical note. It is preferable to use version 2 with brackets([]) if you would like to group together atoms by adding them since you will not clutter up the name space with a lot of parameters you will never use by assigning variables to each atom.

Composition coupling

One common task for coupling atoms together is for fitting compositions instead of occupancies. This is especially true for solid-solid interfaces where the total occupancy usually can be considered to be 1.0. Lets assume we have a LaAlO/SrTiO interface and we want to fit the composition of the interface layer. First we define a interface Slab consisting of a mix the two elements.

slato =  model.Slab(c = 1.05)
# 5.a Define the atoms
slato.add_atom('La', 'La', 0.0, 0.0, 0.0, 0.08, 1.0, 1)
slato.add_atom('Sr', 'Sr', 0.0, 0.0, 0.0, 0.08, 1.0, 1)
slato.add_atom('Al', 'Al', 0.5, 0.5, 0.5, 0.08, 1.0, 1)
slato.add_atom('Ti', 'Ti', 0.5, 0.5, 0.5, 0.08, 1.0, 1)
slato.add_atom('O1', 'O', 0.5, 0.5, 0.0, 0.08, 1.0, 1)
slato.add_atom('O2', 'O', 0.5, 0.0, 0.5, 0.08, 1.0, 2)

So for creating a composition pair of the Sr and La we type:

sl = slato['La']|slato['Sr']

The AtomGroup will now have member functions sl.setcomp to set the La composition and sl.setoc to set the occupancy. Note that the occupancy is calculated as: oc_Sr_ = (1 - comp)*oc, oc_La_ = comp*oc. Per default the positions and debye waller parameters are not coupled (In this case they will refer to the La atom). To do this use the exclusive or operator instead

sl = slato['La']^slato['Sr']

, which couples all parameter through sl. The distinction becomes more important if one has more than two atoms and do nesting of the operators (assume we have an Y atom as well):

sl = slato['La']^slato['Sr']
ysl = slato['Y']^sl

will connect all positional and Debye-Waller parameters together but…

sl = slato['La']^slato['Sr']
ysl = slato['Y']|sl

this will couple the positions of La and Sr but not to the Y. These operators offer a large degree of flexibility but they also probably create hard to track error - so they should be used wisely.

General coupling

As you probably understand the concept with AtomGroup`s cant take you all the way to complete generality. This naturally more complicated but incredibly powerful as _anything_ can be implemented. The class which is used a container class for this is called UserVars. To use it you have to import it by inserting the following text (preferably in the top of the script)

from models.utils import UserVars

Next you create a new object of the class:

cp = UserVars()

You can as many of the objects as possible, a good way to structure your parameters, with as many parameters as possible inside. A new parameter is created by calling the method new_var:

cp.new_var('comp', 0.5)

Now there is a new variable with the name cp.comp and it is initialized to have a default value of 0.5. It can also be accessed by right clicking on the grid and choose it as a fitting parameter. However, now we have to connect the parameter to our sample. Lets assume out model has been extended to include an interdiffused slab between the laouc and stouc.

# 5 Creating an LaAlO3 unit cell - note asymmetric unit cell
laouc =  model.Slab(c = 1.05)
# 5.a Define the atoms
laouc.add_atom('La', 'la', 0.0, 0.0, 0.0, 0.08, 1.0, 1)
laouc.add_atom('Al', 'al', 0.5, 0.5, 0.5, 0.08, 1.0, 1)
laouc.add_atom('O1', 'o', 0.5, 0.5, 0.0, 0.08, 1.0, 1)
laouc.add_atom('O2', 'o', 0.5, 0.0, 0.5, 0.08, 1.0, 2)
# Create an SrTiO3 unit cell
stouc =  model.Slab(c = 1.05)
# 5.a Define the atoms
stouc.add_atom('Sr', 'sr', 0.0, 0.0, 0.0, 0.08, 1.0, 1)
stouc.add_atom('Al', 'al', 0.5, 0.5, 0.5, 0.08, 1.0, 1)
stouc.add_atom('O1', 'o', 0.5, 0.5, 0.0, 0.08, 1.0, 1)
stouc.add_atom('O2', 'o', 0.5, 0.0, 0.5, 0.08, 1.0, 2)
# Create an interface layer
lstouc =  model.Slab(c = 1.05)
# 5.a Define the atoms
lstouc.add_atom('Sr', 'sr', 0.0, 0.0, 0.0, 0.08, 1.0, 1)
lstouc.add_atom('La', 'la', 0.0, 0.0, 0.0, 0.08, 1.0, 1)
lstouc.add_atom('Al', 'al', 0.5, 0.5, 0.5, 0.08, 1.0, 1)
lstouc.add_atom('O1', 'o', 0.5, 0.5, 0.0, 0.08, 1.0, 1)
lstouc.add_atom('O2', 'o', 0.5, 0.0, 0.5, 0.08, 1.0, 2)

# 6 Symmetry operations
p4 = [model.SymTrans([[1, 0],[0, 1]]), model.SymTrans([[-1, 0],[0, -1]]),
           model.SymTrans([[0, -1],[1, 0]]), model.SymTrans([[0, 1],[-1, 0]]) ]

# 7 Creating the sample with one unit cell
sample = model.Sample(inst, bulk, [stouc, lstouc, laouc], unitcell)

# Couple the parameter!
# Create the UserVar
cp = UserVars()
cp.new_var('comp', 0.5)

So the explicit coupling has to be done within the Sim function otherwise the coupling will not work when you fit it at a later stage. So the Sim function would change to

# 9 Define the Sim function
def Sim(data):
    # Parameter coupling goes before simulation!
    lstouc.setSroc(1.0 - cp.comp)
    lstouc.setLaoc(cp.comp)
    I = []
    #9.a loop through the data sets
    for data_set in data:
        # 9.b create all the h,k,l values for the rod (data_set)
        h = data_set.extra_data['h']
        k = data_set.extra_data['k']
        l = data_set.x
        # 9.c. Calculate the struct factor
        f = sample.calc_f(h, k, l)
        # 9.d Calculate the intensity
        i = abs(f)**2
        # 9.e Append the calculated intensity to the list I
        I.append(i)
    return I

As seen above the different occupancies are set with the help of the cp.comp. Note that the parameter coupling is within the Sim function and before we calculate the rods. Also a good code of conduct is never assign a user variable a value from within the Sim function. This can cause a very erratic behavior when fitting (the evaluation depends on its history).

I hope you have understood that with UserVars you can do just about anything. Unfortunately, this will also involve a lot of typing and long scripts. If you want to know more about the general philosophy of the layout with UserVars and the simulation the tutorial about writing models might be good next step Writing a custom model.

Exporting Results

The Sample class allows export of the resulting atom position in the XYZ format. The most convenient way to perform this export is to use the UserFuncs plugin with the following export function in your script:

def export_xyz(xyz_file: TextIO):
    sample.export_xyz(xyz_file)

After simulating the “User funcs” menu should show the entry “export_xyz”. Activating this menu will open a file selection dialog and the selected file will be used to save the XYZ data to.