import escape as esc
esc.require('0.9.6')
import numpy as np
from escape.utils import show
import warnings
warnings.filterwarnings("ignore")
In this notebook we demonstrate modelling of small angle scattering data from bimodal silica particles. More details as well as experimental data you can find in the original paper from Applied Crystallography Journal
http://journals.iucr.org/j/issues/2015/05/00/vg5026/
This notebook is a demonstration of the core functionality of ESCAPE package.
The investigated population of silica nanoparticles consists of particles of two sizes: small particles of type "A" and large particles of type "B", each with its own Gaussian distribution of sizes with characteristic mean radii $R_A$, $R_B$ and their distribution widths $\sigma_A$ and $\sigma_B$.
The resulted intensity is calculated as a sum of the scattering intensity contributed by each scatterer, i.e. each assembly of particles as following
$I_{mod}(q)=\sum_c^N\left[\int_a^b P(q, x)f(x)dx\right]S(q)$
Each contribution $c$ consists of a form-factor $P(q, x)$ determining the shape of a scatterer, $f(x)$ - distribution of certain parameter of the form-factor in the case of disperse aspects. $S(q)$ is a structure factor for each contribution, which is not relevant for our case and is equal unity.
We start implementation of our model with definition of model parameters. In the comments below we give a short description of each parameter. For the details please have a look on the original publication.
#Mean Radiuses of silica particles
R0A = esc.par("Mean_Radius0A", 8.52224, units="nm")
R0B = esc.par("Mean_Radius0B", 37.6503, units="nm")
#We still need a normal radius parameter to perform integration of intensity
RA = esc.par("RadiusA", 8.52224, units="nm")
RB = esc.par("RadiusB", 37.6503, units="nm")
#widths of radii distributions
s_RA = esc.par("Sigma_RA", 2.00668)
s_RB = esc.par("Sigma_RB", 8.29554)
#contrast with water
eta = esc.par("Contrast", 1.017, scale=1e11,
fixed=True, units = "cm^-2")
#Represent number concentration of particles
NA = esc.par("NA", 1.02979, scale=1e-27)
NB = esc.par("NB", 6.5107, scale=1e-31)
#Background
B = esc.par("Background", 0, userlim=[0, 1], fixed=True)
Now we can implement form-factors for our model. These are standard form-factors for homegeneous particles with spherical shape.
#definition of functors for intensities
q = esc.var("q")
PA = esc.pow(4.0/3.0 * np.pi * esc.pow(RA, 3) * eta *
(3 * (esc.sin(q*RA)-q*RA*esc.cos(q*RA))/esc.pow(q*RA, 3.0)), 2)
PB = esc.pow(4.0/3.0 * np.pi * esc.pow(RB, 3) * eta *
(3 * (esc.sin(q*RB)-q*RB*esc.cos(q*RB))/esc.pow(q*RB, 3.0)), 2)
$P_A$ and $P_B$ are the form-factors for the corresponding scatterers without taking into account the size distribution. Since ESCAPE operates mostly with simple functor-type objects, the user has the ability to look at and analyze the behavior of each object, before the final simulation. For example, to visually estimate the contribution of any part of the model or locate a calculation error.
qrange = np.linspace(0.0001, 5, 1000)
show(PA, coordinates=qrange, ylog=True, xlog=True, xlabel="Q, 1/nm", ylabel="PA")
show(PB, coordinates=qrange, ylog=True, xlog=True, xlabel="Q, 1/nm", ylabel="PB")
After implementing of form-factors and checking their curves we are ready to finalize our model and to implement integration over size distribution. Here we had to reformulate the expressions, because in ESCAPE there is no Gaussian number-weighted size distribution. For that purposes we introduced expressions for $C_A$ and $C_B$ normalization constants which we will use as multipliers for the form-factor expressions.
#ESCAPE Gaussian distribution requires FWHM
fwhm_RA = s_RA*2.355
fwhm_RB = s_RB*2.355
#Normalization constant for number size Gaussian distribution
CA = 2.0*NA/(1.0+esc.erf(R0A/(np.sqrt(2.0)*s_RA)))
CB = 2.0*NB/(1.0+esc.erf(R0B/(np.sqrt(2.0)*s_RB)))
#convolution over size distribution
PAr = esc.average_normal(CA*PA, fwhm_RA, RA, R0A, epsabs=1e-8, epsrel=1e-8, maxiter=30, numpoints=61, numstd=5)
PBr = esc.average_normal(CB*PB, fwhm_RB, RB, R0B, epsabs=1e-8, epsrel=1e-8, maxiter=30, numpoints=61, numstd=5)
#Intensity plus background
I = PAr + PBr + B
show(I, coordinates=qrange, ylog=True, xlog=False, xlabel="Q, 1/nm", ylabel="I")
We are ready now to perform fitting of experimental data. To do that we need to open the data file, which was taken from the supplementary material to the original article. To open the file we will use the numpy library and the loadtxt function. The data will be loaded into the arrays x, y, err. Then we create a data object Escape dobj. The next step is to create a model. Here everything is simple, we will use the function model. The arguments residuals_scale and weight_type are relevant to the calculation of the cost function. Weight_type set to "data" means that the data errors will be used for the weights.
#data
x, y, err=np.loadtxt("data/SilicaS2929.pdh", unpack=True)
dobj = esc.data("SilicaS2929.pdh", x, y, err, copy=True)
#model
mobj = esc.model("Model", I, dobj, residuals_scale="none", weight_type="data")
It's a good time to create and run an optimizer. We will start by testing the well-known Levenberg-Marquardt optimizer by shaking out the initial parameter values beforehand.
opt = esc.levmar("LM", mobj, maxiter=100, xtol=1e-15, ftol=1e-15, gtol=1e-15, nupdate=1)
opt.shake()
opt()
show(opt, ylog=True, xlog=True, xlabel="Q, 1/nm", ylabel="Intensity")
opt
The values of the parameters and their errors are very close to the values presented in the article, although there is some expected discrepancy. We try now Differential Evolution algorithm.
#Differential evolution fit with final polishing of parameters with LMA optimizer
opt = esc.diffevol("DiffEvol", mobj, maxiter=50, polish_final_maxiter=50,
polish_candidate_maxiter=0, minconv=1e-5, nupdate=1)
opt.shake()
opt()
show(opt, ylog=True, xlabel="Q, 1/nm", ylabel="Intensity")
opt
The results of LM and DE optimizers for the presented data are in a good agreement with results presented in the original paper.