import numpy as np
import escape as esc
esc.require('0.9.7')
from escape.utils.widgets import show
Specular reflectivity. Parratt32¶
In this notebook using a simple model of Ni layer on Si substrate, we demonstrate how one can fit a Neutron Specular Reflectivity Data. The experimental data was generated using famous Parratt32 software. A poisson noise has ben added to the data. The data and errors have been rescaled and normalized. In Parratt32 program the user provides SLD parameters, real and imaginary parts and fit them separately. In industial packages for the Specular reflectivity fitting, however, one can meet another approach using single mass density parameter. In this case the SLDRe and SLDIm are calulated from the density at each iteration during the fit. This approach reduces the number of fitting parameters, but requires a certainty in sample composition, i.e. atoms and their concentrations, which cannot be always garanteed. The ESCAPE reflectivity module can work with both approaches, moreover in ESCAPE user can choose which approach to use for each material. It there is no confidence about the material quality and its composition, it is adviced to fit SLD parameter instead of mass density.
Below we give a description of our model Ni layer with thickness of 10 nm and roughness 1.5 nm.
Next we create layers and sample.
sample = esc.multilayer("Sample", formula="Ni(10nm,1.5nm)//Si", bydensity=True)
src=esc.neutrons(0.45, units="nm")
#let's have a look on the profile
show(sample, source=src)
Now we create calculation kernel for specular reflectivity
Qz=esc.var("qz")
R = esc.specrefl("Specrefl", Qz, sample, "matrix", source=src)
We read the data from the file, generated using Parratt32 software. The generated data is normalized and has undefined errors. In order to make the data to be closer to real experimental data and to obtain realistic fit errors for our fit parameters, we multiply the data by I0, apply poisson noise and find errors for every point folllowed by normalizing the data again. In general case, the I0 can be also a parameter.
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("Ni/Si", qz*10, y, err, copy=True)
Below the model is created. The residuals scale is set to none, which means that there is no scaling. There is also a support for logarithmic log scaling of residuals. The weight_type parameter can have data value or none, where for the former case the data errors are taken into account or no weights for the latter case.
mobj = esc.model("Model", R, dobj, residuals_scale="none", weight_type="data")
show(mobj, ylog=True, xlog=False, xlabel="Q[1/nm]", ylabel="Intensity")
Let's first use Levenberg-Marquardt optimizer. For the given parameter ranges it behaves quite stable independently from starting parameter values. Shake button allows to choose initial parameters randomly.
opt = esc.levmar("LM", mobj, maxiter=1000, xtol=1e-15, ftol=1e-15, gtol=1e-15, nupdate=1)
opt.shake()
opt()
show(opt, ylog=True, xlog=False, xlabel="Q[1/nm]", ylabel="R")
#Let's print parameters and their errors after fit
opt
Name: LM Parameters number: 4 Parameter Value +- Error Units Fixed Ni density 8.8979 +- 0.0088536 g/cm^3 0 Layer-Ni: Thickness 10.006 +- 0.01058 nm 0 Layer-Ni: Roughness Sigma 1.5014 +- 0.016831 nm 0 Si density 2.3347 +- 0.0088124 g/cm^3 0
Below we demonstrate Differential Evolutionary optimizer. For such a simple model and data with nearly ideal quality DE doesn't give any advantage comparing to LM.
opt = esc.diffevol("DiffEvol", mobj, popsize=10, maxiter=1000,
mutation=0.5, crossover=0.8, minconv=1e-3, mincost=1e-8, nupdate=1,
polish_final_maxiter=150, polish_candidate_maxiter=0)
opt.shake()
opt()
show(opt, ylog=True)
opt
Name: DiffEvol Parameters number: 4 Parameter Value +- Error Units Fixed Ni density 8.8979 +- 0.0088536 g/cm^3 0 Layer-Ni: Thickness 10.006 +- 0.01058 nm 0 Layer-Ni: Roughness Sigma 1.5014 +- 0.016831 nm 0 Si density 2.3347 +- 0.0088124 g/cm^3 0