Simulation examples

In the following a few code-snippets are shown which should help you getting started with reflectivity and diffraction simulations using xrayutilities. All simulations in xrayutilities are for layers systems and currently there are no plans to extend this to other geometries. Note that not all of the codes shown in the following will be run-able as stand-alone scripts. For fully running scripts look in the examples directory in the download.

Building Layer stacks for simulations

The basis of all simulations in xrayutilities are stacks of layers. Therefore several functions exist to build up such layered systems. The basic building block of all of them is a Layer object which takes a material and its thickness in ångström as initializing parameter.

import xrayutilities as xu
lay = xu.simpack.Layer(xu.materials.Si, 200)

In the shown example a silicon layer with 20 nm thickness is created. The first argument is the material of the layer. For diffraction simulations this needs to be derived from the Crystal-class. This means all predefined materials in xrayutitities can be used for this purpose. For x-ray reflectivity simulations, however, also knowing the chemical composition and density of the material is sufficient.

A 5 nm thick metallic CoFe compound layer can therefore be defined sublayers

rho_cf = 0.5*8900 + 0.5*7874  # mass density in kg/m^3
mCoFe = xu.materials.Amorphous('CoFe', rho_cf)
lCoFe = xu.simpack.Layer(mCoFe, 50)

Note

The Layer object can have several more model dependent properties discussed in detail below.

When several layers are defined they can be combined to a LayerStack which is used for the simulations below.

1sub = xu.simpack.Layer(xu.materials.Si, float('inf'))
2lay1 = xu.simpack.Layer(xu.materials.Ge, 200)
3lay2 = xu.simpack.Layer(xu.materials.SiO2, 30)
4ls = xu.simpack.LayerStack('Si/Ge', sub, lay1, lay2)
5# or equivalently
6ls = xu.simpack.LayerStack('Si/Ge', sub + lay1 + lay2)

The last two lines show two different options of creating a stack of layers. As is shown in the last example the substrate thickness can be infinite (see below) and layers can be also stacked by summation. For creation of more complicated superlattice stacks one can further use multiplication

lay1 = xu.simpack.Layer(xu.materials.SiGe(0.3), 50)
lay2 = xu.simpack.Layer(xu.materials.SiGe(0.6), 40)
layerstack = xu.simpack.LayerStack('Si/SiGe SL', sub + 5*(lay1 + lay2))

Pseudomorphic Layers

All stacks of layers described above use the materials in the layer as they are supplied. However, epitaxial systems often adopt the inplane lattice parameter of the layers beneath. To mimic this behavior you can either supply the Layer objects which custom Crystal objects which have the appropriate lattice parameters or use the PseudomorphicStack* classes which to the adaption of the lattice parameters automatically. In this respect the ‘relaxation’ parameter of the Layer class is important since it allows to create partially/fully relaxed layers.

1sub = xu.simpack.Layer(xu.materials.Si, float('inf'))
2buf1 = xu.simpack.Layer(xu.materials.SiGe(0.5), 5000, relaxation=1.0)
3buf2 = xu.simpack.Layer(xu.materials.SiGe(0.8), 5000, relaxation=1.0)
4lay1 = xu.simpack.Layer(xu.materials.SiGe(0.6), 50, relaxation=0.0)
5lay2 = xu.simpack.Layer(xu.materials.SiGe(1.0), 50, relaxation=0.0)
6# create pseudomorphic superlattice stack
7pls = xu.simpack.PseudomorphicStack001('SL 5/5', sub+buf1+buf2+5*(lay1+lay2))

Note

As indicated by the function name the PseudomorphicStack currently only works for (001) surfaces and cubic materials. Implementations for other surface orientations are planned.

If you would like to check the resulting lattice objects of the different layers you could use:

for l in pls:
    print(l.material.lattice)

Special layer types

So far one special layer mimicking a layer with gradually changing chemical composition is implemented. It consists of several thin sublayers of constant composition. So in order to obtain a smooth grading one has to select enough sublayers. This however has a negativ impact on the performance of all simulation models. A tradeoff needs to found! Below a graded SiGe buffer is shown which consists of 100 sublayers and has total thickness of 1 µm.

1buf = xu.simpack.GradedLayerStack(xu.materials.SiGe,
2                                  0.2,  # xfrom Si0.8Ge0.2
3                                  0.7,  # xto Si0.3Ge0.7
4                                  100,  # number of sublayers
5                                  10000,  # total thickness
6                                  relaxation=1.0)

Setting up a model

This sectiondescribes the parameters which are common for all diffraction models in xrayutilties-simpack. All models need a list of Layers for which the reflected/diffracted signal will be calculated. Further all models have some common parameters which allow scaling and background addition in the model output and contain general information about the calculation which are model-independent. These are

  • ‘experiment’: an Experiment/HXRD object which defines the surface geometry of the model. If none is given a default class with (001) surface is generated.

  • ‘resolution_width’: width of the Gaussian resolution function used to convolute with the data. The unit of this parameters depends on the model and can be either in degree or 1/AA.

  • ‘I0’: is the primary beam flux/intensity

  • ‘background’: is the background added to the simulation after it was scaled by I0

  • ‘energy’: energy in eV used to obtain the optical parameters for the simulation. The energy can alternatively also be supplied via the ‘experiment’ parameter, however, the ‘energy’ value overrules this setting. If no energy is given the default energy from the configuration is used.

The mentioned parameters can be supplied to the constructor method of all model classes derived from LayerModel, which applies to all examples mentioned below.

m = xu.simpack.SpecularReflectivityModel(layerstack, I0=1e6, background=1,
                                         resolution_width=0.001)

Reflectivity calculation and fitting

This section shows the calculation and fitting of specular x-ray reflectivity curves as well as the calculation of diffuse x-ray reflectivity curves/maps.

Specular x-ray reflectivity

For the specular reflectivity models currently only the Parrat formalism including non-correlated roughnesses is implemented. A minimal working example for a reflectivity calculation follows.

 1import numpy
 2# building a stack of layers
 3sub = xu.simpack.Layer(xu.materials.GaAs, float('inf'), roughness=2.0)
 4lay1 = xu.simpack.Layer(xu.materials.AlGaAs(0.25), 75, roughness=2.5)
 5lay2 = xu.simpack.Layer(xu.materials.AlGaAs(0.75), 25, roughness=3.0)
 6pls = xu.simpack.PseudomorphicStack001('pseudo', sub+5*(lay1+lay2))
 7
 8# reflectivity calculation
 9m = xu.simpack.SpecularReflectivityModel(pls, sample_width=5, beam_width=0.3)
10ai = numpy.linspace(0, 5, 10000)
11Ixrr = m.simulate(ai)

In addition to the layer thickness also the roughness and density (in kg/m^3) of a Layer can be set since they are important for the reflectivity calculation. This can be done upon definition of the Layer or also manipulated at any later stage. Such x-ray reflectivity calculations can also be fitted to experimental data using the FitModel class which is shown in detail in the example below (which is also included in the example directory). The fitting is performed using the lmfit Python package which needs to be installed when you want to use this fitting function. This package allows to build complicated models including bounds and correlations between parameters.

 1import lmfit
 2import numpy
 3
 4import xrayutilities as xu
 5
 6# load experimental data
 7ai, edata, eps = numpy.loadtxt('data/xrr_data.txt'), unpack=True)
 8ai /= 2.0
 9
10# define layers
11# SiO2 / Ru(5) / CoFe(3) / IrMn(3) / AlOx(10)
12lSiO2 = xu.simpack.Layer(xu.materials.SiO2, numpy.inf, roughness=2.5)
13lRu = xu.simpack.Layer(xu.materials.Ru, 47, roughness=2.8)
14rho_cf = 0.5*8900 + 0.5*7874
15mat_cf = xu.materials.Amorphous('CoFe', rho_cf)
16lCoFe = xu.simpack.Layer(mat_cf, 27, roughness=4.6)
17lIrMn = xu.simpack.Layer(xu.materials.Ir20Mn80, 21, roughness=3.0)
18lAl2O3 = xu.simpack.Layer(xu.materials.Al2O3, 100, roughness=5.5)
19
20# create model
21m = xu.simpack.SpecularReflectivityModel(lSiO2, lRu, lCoFe, lIrMn, lAl2O3,
22                                         energy='CuKa1', resolution_width=0.02,
23                                         sample_width=6, beam_width=0.25,
24                                         background=81, I0=6.35e9)
25
26# embed model in fit code
27fitm = xu.simpack.FitModel(m, plot=True, verbose=True)
28
29# set some parameter limitations
30fitm.set_param_hint('SiO2_density', vary=False)
31fitm.set_param_hint('Al2O3_density', min=0.8*xu.materials.Al2O3.density,
32                    max=1.2*xu.materials.Al2O3.density)
33p = fitm.make_params()
34fitm.set_fit_limits(xmin=0.05, xmax=8.0)
35
36# perform the fit
37res = fitm.fit(edata, p, ai, weights=1/eps)
38lmfit.report_fit(res, min_correl=0.5)
39# export the fit result for the full data range (Note that only data between
40# xmin and xmax were actually optimized)
41numpy.savetxt(
42    "xrrfit.dat",
43    numpy.vstack((ai, res.eval(res.params, x=ai))).T,
44    header="incidence angle (deg), fitted intensity (arb. u.)",
45)

This script can interactively show the fitting progress and after the fitting shows the final plot including the x-ray reflectivity trace of the initial and final parameters.

XRR fitting output

The picture shows the final plot of the fitting example shown in one of the example scripts.

After building a SpecularReflectivityModel is built or fitted the density profile resulting from the thickness and roughness of layers can be plotted easily by

m.densityprofile(500, plot=True)  # 500 number of points
XRR density profile resulting from the XRR fit shown above

Diffuse reflectivity calculations

For the calculation of diffuse x-ray reflectivity the LayerStack is built equally as shown before. The only difference is that an additional parameter for the lateral correlation length of the roughness can be included: lat_correl. The DiffuseReflectivityModel also takes special parameters which change the vertical correlection length and the way how the diffuse reflectivity is calculated (to be document in more detail). For a Si/Ge superlattice with 5 periods the calculation of the diffuse reflectivity signal at the specular rod is calculated using the simulate() method. A map of the diffuse reflectivity which can be obtained in the coplanar reflection plane can be calculated with the simulate_map() method.

 1from matplotlib.pylab import *
 2import xrayutilities as xu
 3sub = xu.simpack.Layer(xu.materials.Si, inf, roughness=1, lat_correl=100)
 4lay1 = xu.simpack.Layer(xu.materials.Si, 200, roughness=1, lat_correl=200)
 5lay2 = xu.simpack.Layer(xu.materials.Ge, 70, roughness=3, lat_correl=50)
 6ls = xu.simpack.LayerStack('SL 5', sub+5*(lay2+lay1))
 7
 8alphai = arange(0.17, 2, 0.001)  # for the calculation on the specular rod
 9qz = arange(0, 0.5, 0.0005)  # for the map calculation
10qL = arange(-0.02, 0.02, 0.0003)
11
12m = xu.simpack.DiffuseReflectivityModel(ls, sample_width=10, beam_width=1,
13                                        energy='CuKa1', vert_correl=1000,
14                                        vert_nu=0, H=1, method=2, vert_int=0)
15d = m.simulate(alphai)
16imap = m.simulate_map(qL, qz)
17
18figure()
19subplot(121)
20semilogy(alphai, d, label='diffuse XRR')
21xlabel('incidence angle (deg)')
22ylabel('intensity (arb. u.)')
23ylim(1e-6, 1e-4)
24
25subplot(122)
26pcolor(qL, qz, imap.T, norm=mpl.colors.LogNorm())
27xlabel(r'Q$_\parallel$ ($\mathrm{\AA^{-1}}$)')
28ylabel(r'Q$_\perp$ ($\mathrm{\AA^{-1}}$)')
29colorbar()
30tight_layout()

The resulting figure shows the simulation result. Currently you have to refer to the docstrings and implementation for further details.

Diffuse x-ray reflectivity of a Si/Ge multilayer

Diffraction calculation

From the very same models as used for XRR calculation one can also perform crystal truncation rod simulations around certain Bragg peaks using various different diffraction models. Depending on the system to model you will have to choose the most apropriate model. Below a short description of the implemented models is given followed by two examples.

Kinematical diffraction models

The most basic models consider only the kinematic diffraction of layers and substrate. Especially the semiinfinite substrate is not well described using the kinematical approximation which results in considerable deviations in close vicinity to substrate Bragg peak with respect to the more acurate dynamical diffraction models.

Such a basic model is employed by

en = 9000  # eV
mk = xu.simpack.KinematicalModel(pls, energy=en, resolution_width=0.0001)
Ikin = mk.simulate(qz, hkl=(0, 0, 4))

A more appealing kinematical model is represented by the KinematicalMultiBeamModel class which implements a true multibeam theory is, however, restricted to the use of (001) surfaces and layer thicknesses will be changed to be a multiple of the out of plane lattice spacing. This is necessary since otherwise the structure factor of the unit cell can not be used for the calculation.

It can be employed by

mk = xu.simpack.KinematicalMultiBeamModel(pls, energy=en,
                                          surface_hkl=(0, 0, 1),
                                          resolution_width=0.0001)
Imult = mk.simulate(qz, hkl=(0, 0, 4))

This model is expected to provide good results especially far away from the substrate peak where the influence of other Bragg peaks on the truncation rod and the variation of the structure factor can not be neglected.

Both kinematical model’s simulate() method offers two keyword arguments with which basic absorption and refraction correction can be added to the basic models.

Note

The kinematical models can also handle a semi-infinitely thick substrate which results in a diverging intensity at the Bragg peak but provides a basic description of the substrates truncation rod.

Dynamical diffraction models

Acurate description of the diffraction from thin films in close vicinity to the diffraction signal from a bulk substrate is only possible using the dynamical diffraction theory. In xrayutilities the dynamical two-beam theory with 4 tiepoints for the calculation of the dispersion surface is implemented. To use this theory you have to supply the simulate() method with the incidence angle in degree. Accordingly the ‘resolution_width’ parameter is also in degree for this model.

resol = 0.001  # resolution in incidence angle
md = xu.simpack.DynamicalModel(pls, energy=en, resolution_width=resol)
Idyn = md.simulate(ai, hkl=(0, 0, 4))

A second simplified dynamical model (SimpleDynamicalCoplanarModel) is also implemented should, however, not be used since its approximations cause mistakes in almost all relevant cases.

The DynamicalModel supports the calculation of diffracted signal for ‘S’ and ‘P’ polarization geometry. To simulate diffraction data of laboratory sources with Ge(220) monochromator crystal one should use:

.. testcode::

import math qGe220 = linalg.norm(xu.materials.Ge.Q(2, 2, 0)) thMono = math.asin(qGe220 * xu.config.WAVELENGTH / (4*math.pi)) md = xu.simpack.DynamicalModel(pls,

Cmono=math.cos(2 * thMono), polarization=’both’)

Idyn = md.simulate(ai, hkl=(0, 0, 4))

Comparison of diffraction models

Below we show the different implemented models for the case of epitaxial GaAs/AlGaAs and Si/SiGe bilayers. These two cases have very different separation of the layer Bragg peak from the substrate and therefore provide good model system for our models.

We will compare the (004) Bragg peak calculated with different models and but otherwise equal parameters. For scripts used to perform the shown calculation you are referred to the examples directory.

(004) of AlGaAs(100nm) on GaAs

XRD simulations of the (004) Bragg peak of ~100 nm AlGaAs on GaAs(001) using various diffraction models

(004) of SiGe(15nm) on Si

XRD simulations of the (004) Bragg peak of 15 nm Si0.4 Ge0.6 on Si(001) using various diffraction models

As can be seen in the images we find that for the AlGaAs system all models except the very basic kinematical model yield an very similar diffraction signal. The second kinematic diffraction model considering the contribution of multiple Bragg peaks on the same truncation rod fails to describe only the ratio of substrate and layer signal, but otherwise results in a very similar line shape as the traces obtained by the dynamic theory.

For the SiGe/Si bilayer system bigger differences between the kinematic and dynamic models are found. Further also the difference between the simpler and more sophisticated dynamic model gets obvious further away from the reference position. Interestingly also the multibeam kinematic theory differs considerable from the best dynamic model. As is evident from this second comparison the correct choice of model for the particular system under condideration is crucial for comparison with experimental data.

Fitting of diffraction data

All diffraction models can be embedded into the FitModel class, which is suitable to refine the model parameters. Below (and in the examples directory) a runnable script is shown which shows the fitting for a pseudomorphic InMnAs epilayer on InAs(001). The fitting is performed using the lmfit Python package which needs to be installed when you want to use this fitting function. As one can see below the set_param_hint() function can be used to set up the respective fit parameters including their boundaries and possible correlation with other parameters of the model. It should be equally possible to fit more complex layer structures, however, I expect that one needs to adjust manually the starting parameters to yield something very reasonable. Since this capabilities are rather new please report back any success/problems you have with this via the mailing list.

 1from copy import deepcopy
 2import xrayutilities as xu
 3from matplotlib.pylab import *
 4
 5# global parameters
 6wavelength = xu.wavelength('CuKa1')
 7offset = -0.035  # angular offset of the zero position of the data
 8
 9# set up LayerStack for simulation: InAs(001)/(In,Mn)As(~25 nm)
10InAs = deepcopy(xu.materials.InAs)  # do not modify internal database
11InAs.a = 6.057
12lInAs = xu.simpack.Layer(InAs, inf)
13InMnAs = xu.materials.Crystal('InMnAs', xu.materials.SGLattice(
14    216, 6.050, atoms=('In', 'Mn', 'As'), pos=('4a', '4a', '4c'),
15    occ=(0.99, 0.01, 1)), cij=InAs.cij)
16lInMnAs = xu.simpack.Layer(InMnAs, 254, relaxation=0)
17pstack = xu.simpack.PseudomorphicStack001('list', lInAs, lInMnAs)
18
19# set up simulation object
20thetaMono = arcsin(wavelength/(2 * xu.materials.Ge.planeDistance(2, 2, 0)))
21Cmono = cos(2 * thetaMono)
22dyn = xu.simpack.DynamicalModel(pstack, I0=1.5e9, background=0,
23                                resolution_width=2e-3, polarization='both',
24                                Cmono=Cmono)
25fitmdyn = xu.simpack.FitModel(dyn)
26fitmdyn.set_param_hint('InMnAs_c', vary=True, min=6.02, max= 6.06)
27fitmdyn.set_param_hint('InAs_a', vary=True)
28fitmdyn.set_param_hint('InMnAs_a', expr='InAs_a')
29fitmdyn.set_param_hint('resolution_width', vary=True)
30params = fitmdyn.make_params()
31
32# plot experimental data
33f = figure(figsize=(7,5))
34d = xu.io.RASFile('inas_layer_radial_002_004.ras.bz2', path='examples/data')
35scan = d.scans[-1]
36tt = scan.data[scan.scan_axis] - offset
37semilogy(tt, scan.data['int'], 'o-', ms=3, label='data')
38
39# perform fit and plot the result
40fitmdyn.lmodel.set_hkl((0, 0, 4))
41ai = (d.scans[-1].data[d.scan.scan_axis] - offset)/2
42fitr = fitmdyn.fit(d.scans[-1].data['int'], params, ai)
43# print(fitr.fit_report())  # to get a summary of the fitted parameters

The resulting figure shows reasonable agreement between the dynamic diffraction simulation and the experimental data.

(004) diffraction data and fitted model of InMnAs epilayer on InAs

Powder diffraction simulations

Powder diffraction patterns can be calculated using PowderModel. A specialized class for the definition of powdered materials named Powder exists. The class constructor takes the materials volume and several material parameters specific for the powder material. Among them are crystallite_size_gauss and strain_gauss which can be used to include the effect of finite crystallite size and microstrain. Texture modelled by the March-Dollase pole density offers the preferred_orientation direction parameter as well as a preferred_orientation_factor variable.

The PowderModel internally uses PowderDiffraction for its calculations which is based on the fundamental parameters approach. The fundamental parameters approach code used here can was originally implemented and documented here: http://dx.doi.org/10.6028/jres. 120.014.c and http://dx.doi.org/10.6028/jres.120.014.

Several setup specific parameters should be adjusted by a user-specific configuration file or by supplying the appropriate parameters using the fpsettings argument of PowderModel.

If the correct settings are included in the config file the powder diffraction signal of a mixed sample of Co and Fe can be calculated with

 1import numpy
 2import xrayutilities as xu
 3
 4tt = numpy.arange(5, 120, 0.01)
 5Fe_powder = xu.simpack.Powder(xu.materials.Fe, 1,
 6                              crystallite_size_gauss=100e-9)
 7Co_powder = xu.simpack.Powder(xu.materials.Co, 5,  # 5 times more Co
 8                              crystallite_size_gauss=200e-9)
 9pm = xu.simpack.PowderModel(Fe_powder, Co_powder, I0=100)
10inte = pm.simulate(tt)
11pm.close()  # after end-of-use

Note that in MS windows and macOS you need to encapsulate this code into a dummy function to allow for the multiprocessing module to work correctly. See the Python documentation for details. The code then must look like:

 1import numpy
 2import xrayutilities as xu
 3from multiprocessing import freeze_support
 4
 5def main():
 6    tt = numpy.arange(5, 120, 0.01)
 7    Fe_powder = xu.simpack.Powder(xu.materials.Fe, 1,
 8                                  crystallite_size_gauss=100e-9)
 9    Co_powder = xu.simpack.Powder(xu.materials.Co, 5,  # 5 times more Co
10                                  crystallite_size_gauss=200e-9)
11    pm = xu.simpack.PowderModel(Fe_powder, Co_powder, I0=100)
12    inte = pm.simulate(tt)
13    pm.close()
14
15if __name__ == '__main__':
16    freeze_support()  # only required on MS Windows
17    main()