In [1]:
import numpy as np
import escape as esc

esc.require("0.9.8")

Loading material database from C:\dev\escape-core\python\src\escape\scattering\..\data\mdb\materials.db

Specular reflectivity. Parratt32¶

This notebook fits a neutron specular reflectivity data set for a simple Ni film on a Si substrate. Specular reflectivity measures how the reflected intensity changes with the normal component of momentum transfer, and it is especially sensitive to layer thickness, interface roughness, and scattering length density.

The data here were generated with Parratt32 and then perturbed with Poisson noise so the example resembles an experimental workflow. ESCAPE can fit either density-based parameters or SLD parameters directly, and the notebook shows how that choice affects the model setup.

Next we create the layers and assemble the sample.

In [2]:
sample = esc.multilayer( formula="Ni(10nm,1.5nm)//Si", bydensity=True)
src = esc.neutrons(0.45, units="nm")
# let's have a look on the profile
sample.show(source=src)
Out[2]:

Now we create the calculation kernel for specular reflectivity.

In [3]:
Qz = esc.var("qz")
R = esc.specrefl(Qz, sample, "matrix", source=src)

We now load the data file generated by Parratt32. The original curve is normalized and does not include realistic uncertainties. To mimic an experiment more closely, we scale the intensity by I0, draw Poisson noise, estimate the counting error for each point, and then normalize the data again. In a more advanced fit, I0 can also be treated as a free parameter.

In [4]:
qz, y, err = np.loadtxt("data/Ni-100A-15A_Si_refl.dat", unpack=True)
I0 = 1e6

y = np.random.poisson(y * I0) / I0
err = np.sqrt(y) / np.sqrt(I0)

dobj = esc.data(qz * 10, y, err, copy=True)

Below we create the model object. The residuals_scale="none" option means that no additional scaling is applied to the residuals. ESCAPE also supports logarithmic residual scaling when that is more appropriate.

The weight_type option controls whether the data uncertainties are used during fitting. Use weight_type="data" to include the error bars, or weight_type="none" to ignore them.

In [5]:
mobj = esc.model( R, dobj, residuals_scale="none", weight_type="data")
mobj.show().config(ylog=True, xlog=False, xlabel="Q[1/nm]", ylabel="Intensity")
Out[5]:

We first use the Levenberg-Marquardt optimizer. For this simple model and the chosen parameter ranges, it is usually stable even when the starting values change. The Shake button can be used to randomize the initial parameters before fitting.

In [6]:
opt = esc.levmar(
    mobj, maxiter=1000, nupdate=1
)
opt.shake()
opt()
opt.show().config_model(ylog=True, xlog=False, xlabel="Q[1/nm]", ylabel="R")
Out[6]:
In [7]:
# Let's print parameters and their errors after fit
opt
Out[7]:
Name: Levenberg-Marquardt 	Parameters number:          4
Parameter           	Value          	+-	Error     	Units     	Fixed
Ni Density          	         8.8981	+-	0.0085523 	          	    0
Ni Thickness        	         9.9952	+-	0.010273  	nm        	    0
Ni Roughness        	         1.4887	+-	0.016218  	nm        	    0
Si Density          	         2.3369	+-	0.0087571 	          	    0

Below we demonstrate the Differential Evolution optimizer. For this simple example, and with data that is already fairly clean, it does not usually offer an advantage over Levenberg-Marquardt.

In [8]:
opt = esc.diffevol(
    mobj,
    popsize=7,
    maxiter=20,
    ftol=1e-3,
    nupdate=1,
    strategy="sqgabin",
    polish_final_maxiter=150,
    polish_candidate_maxiter=0,
)
opt.shake()
opt()
opt.show().config_model(ylog=True)
Out[8]:
In [9]:
opt
Out[9]:
Name: Differential Evolution	Parameters number:          4
Parameter           	Value          	+-	Error     	Units     	Fixed
Ni Density          	         8.8981	+-	0.0085523 	          	    0
Ni Thickness        	         9.9952	+-	0.010273  	nm        	    0
Ni Roughness        	         1.4887	+-	0.016218  	nm        	    0
Si Density          	         2.3369	+-	0.0087571 	          	    0
In [ ]: