Fitting QENS data

Fitting QENS data#

Previously, some quasi-elastic neutron scattering (QENS) data has been simulated and reduced, and can now be analysed with easydynamics. If you get stuck, the documentation can hopefully help.

In the reduction step we produced two data sets:

an elastic sample, which scatters neutrons without changing their energy. We will use it to measure the instrument resolution.
a quasi-elastic sample, in which (some of) the scatterers diffuse. This is the data we actually want to interpret.

The workflow is as follows: We first load the data into EasyDynamics and inspect it. We then determine the resolution function using the elastic data. Next, we fit the QENS data to an empirical model and inspect the resulting fit. Finally, we extract physically meaningful information from these fits.

import os

import numpy as np
import scipp as sc

import easydynamics as edyn

# Quiz questions for this notebook
from jupyterquiz import display_quiz
from quizlib import qens as quiz

# Make the plots interactive
%matplotlib widget

Load and prepare the data#

First we load the elastic sample data. We will use it to determine the instrument resolution, in the same way that a vanadium measurement is used in a real experiment.

ℹ️ If you did not complete the QENS data reduction yesterday, don’t worry: we have prepared data for you that you can download.

Our reduced data are histograms of neutron counts, and they need a small adjustment before fitting. Future versions of EasyDynamics will not need this step.

Some of the bins have zero counts, and consequently zero variance. When fitting, each data point is weighted by 1/variance, which therefore diverges, giving those bins infinite weight and destabilising the fit. The solution that we will employ is to give all these bins a variance of 1 instead.

def prepare_data(experiment):
    """Give bins with too small variance a variance of 1, so that their weight in the fit stays finite."""
    data = experiment.data
    indices = data.variances <= 0.5
    data.variances[indices] = 1.0
    experiment.data = data

We load the data into an Experiment object and use its plot function to inspect the data. We can either plot a 2d map or use a slicer to see the intensity as function of energy for each Q.

filename = "../4-reduction/energy_QENS_elastic.h5"

# If you have not completed the reduction exercise, uncomment the following lines to fetch the pre-prepared data
# import utils
# filename = utils.fetch_data('4-reduction/energy_QENS_elastic.h5'

elastic_experiment = edyn.Experiment(display_name="Elastic")
elastic_experiment.load_hdf5(filename=filename)
prepare_data(elastic_experiment)

elastic_experiment.plot_data(slicer=True)

Look closely at the position of the elastic peak in the plot above (you may need to zoom in).

display_quiz(quiz.q1)

Next we load the quasi-elastic sample data. This is the data whose dynamics we want to understand. There is clearly a lot more going on here.

filename2 = "../4-reduction/energy_QENS_sample.h5"

# If you have not completed the reduction exercise, uncomment the following lines to fetch the pre-prepared data
# import utils
# filename2 = utils.fetch_data('4-reduction/energy_QENS_sample.h5')

qe_experiment = edyn.Experiment(display_name="QuasiElastic")
qe_experiment.load_hdf5(filename=filename2)
prepare_data(qe_experiment)

qe_experiment.plot_data(slicer=True)

Step 2: Fit the quasi-elastic sample at each Q#

We are now happy with the resolution and can turn to the quasi-elastic sample. Looking at the data we plotted earlier, it has a sharp elastic peak, a broader quasi-elastic peak, and a small flat background.

We describe it with a SampleModel containing:

a DeltaFunction for the elastic (immobile) scattering, and
a Lorentzian for the quasi-elastic broadening caused by motion.

There are two reasons for choosing a Lorentzian for the quasi-elastic scattering. The first is simple: the data looks like a Lorentzian. The second is that a Lorentzian, \(I=\frac{A}{\pi}\frac{\Gamma}{\Gamma^2+E^2}\), correspond to time correlations decaying exponentially, i.e. as \(e^{-t/\tau}\), where \(\tau=\hbar/\Gamma\) is the decay time. This is exactly how we expect correlations to decay in a diffusive process.

Exercise: create a new SampleModel for the quasi-elastic sample#

Follow the steps to create a SampleModel like we did above, but this time create both a DeltaFunction and a Lorentzian and append them. Give the DeltaFunction a reasonable start value for the area, and the Lorentzian both an area and width (half width at half max). Remember to set the y_unit to ‘counts’ everywhere.

Solution:

We build a new InstrumentModel, and this time we give it a resolution: the SampleModel from our elastic fit. All of its parameters are automatically fixed and the resolution is normalised to have area 1.

background_model_qe = edyn.BackgroundModel(y_unit="counts")

poly_qe = edyn.Polynomial(coefficients=[0.0], name="Background", y_unit="counts")
poly_qe.coefficients[0].min = 0.0
poly_qe.coefficients[0].fixed = True
background_model_qe.append_component(poly_qe)

instrument_model_qe = edyn.InstrumentModel(
    background_model=background_model_qe,
    resolution_model=elastic_analysis.sample_model,
    energy_offset=1e-3,
)

Exercise: collect the data, SampleModel and InstrumentModel in an Analysis object#

Solution:

Analysis handles the convolution of the sample_model with the resolution. The calculation is analytical where possible and numerical otherwise.

Before fitting, it is a good idea to check the start guesses by plotting the data together with the model. If the model does not at all resemble the data, then go back up and change the default values for the SampleModel components and try again. This process may take a few steps to get right.

qe_analysis.plot_data_and_model()

Adjust the start guesses until they look reasonable, then fit every Q independently and plot the result like we did above.

Solution:

First, inspect the fit. Does it look good at all Q? If not, go back and adjust the start guesses and try again. If yes, please continue.

The interesting parameters are the width and area of the Lorentzian. Let us plot them as a function of Q. xmin and xmax set the minimum and maximum of the x axis, while vmin and vmax set the minimum and maximum on the vertical axis. We can plot multiple parameters in the same figure, but they will share the same y axis, so we make two separate figures here.

qe_analysis.plot_parameters(
    names=["Lorentzian width"], vmin=0, vmax=0.03, xmin=0, xmax=2.1
)

qe_analysis.plot_parameters(
    names=["Lorentzian area"], vmin=0, vmax=2.0, xmin=0, xmax=2.1
)

Sometimes, the fitter produces no error bars on the fit parameters. We are working on fixing this, but have to work around it in the meantime. Run the cell below if needed.

def fix_missing_variances(analysis):
    """Fix missing variances on fitted parameters in an analysis."""
    for param in analysis.get_all_parameters():
        if param.variance == 0.0:
            param.variance = 0.00001 * abs(param.value)


fix_missing_variances(qe_analysis)

Before reading on, study the plot of the Lorentzian width as a function of Q.

display_quiz(quiz.q7)

display_quiz(quiz.q8)

Step 3: Fit a jump-diffusion model to the widths#

The width does not simply grow like \(Q^2\): it levels off at high \(Q\). This is the signature of jump diffusion, in which a particle sits still for a residence time \(\tau\) and then jumps to a new site. The half-width of the quasi-elastic Lorentzian is

\[ \Gamma(Q) = \frac{\hbar\,D\,Q^2}{1 + D\,\tau\,Q^2}, \]

where \(D\) is the diffusion coefficient and \(\tau\) is the residence (relaxation) time. At low \(Q\) this reduces to ordinary diffusion, \(\Gamma \approx \hbar D Q^2\), while at high \(Q\) it saturates at the plateau \(\hbar/\tau\).

Our width curve shows exactly this behaviour: a \(Q^2\) rise that bends over towards a plateau within the measured range, so the data constrain both parameters: the low-\(Q\) slope fixes \(D\), and the high-\(Q\) plateau fixes \(\tau\).

As a first step we fit the jump-diffusion model to the fitted Lorentzian parameters we obtained above. To do this, we create a JumpTranslationalDiffusion model.

Next, we need to tell EasyDynamics which parameter(s) to fit to this model. The JumpTranslationalDiffusion model predicts both the area and width of the Lorentzian. We therefore create a FitBinding, which we use to tell EasyDynamics that the ‘area’ of the model should model the parameter named Lorentzian area, while the width should model the parameter named Lorentzian width. We do this using a dictionary as shown below.

We then create a ParameterAnalysis object. It takes the parameters from our analysis and a single or a list of FitBindings to fit. We can simply give it the whole Analysis object, and it will extract the parameters itself.

jump_diffusion_model = edyn.JumpTranslationalDiffusion(
    name="Jump Translational Diffusion",
    diffusion_coefficient=4.6e-10,  # m^2/s
    relaxation_time=22.0,  # ps
    scale=0.5,
    y_unit="counts",
)

binding = edyn.FitBinding(
    model=jump_diffusion_model,
    targets={
        "area": "Lorentzian area",
        "width": "Lorentzian width",
    },
)

parameter_analysis = edyn.ParameterAnalysis(
    parameters=qe_analysis,
    bindings=binding,
)

We first plot the start guess to see if it is reasonable.

parameter_analysis.plot(names=["Lorentzian width"], xmin=0, xmax=2.1, vmin=0, vmax=0.02)

Then we fit and plot again.

parameter_analysis.fit()
parameter_analysis.plot(names=["Lorentzian width"], xmin=0, xmax=2.1, vmin=0, vmax=0.02)

parameter_analysis.plot(names=["Lorentzian area"], xmin=0, xmax=2.1, vmin=0, vmax=1.2)

We can read off the fitted jump-diffusion parameters, with uncertainties.

parameter_analysis.get_all_parameters()

[<Parameter 'diffusion_coefficient': 4.288e-10 ± 3.148e-11 m^2/s, bounds=[0.0:inf]>,
 <Parameter 'relaxation_time': 29.5564 ± 1.4116 ps, bounds=[0.0:inf]>,
 <Parameter 'scale': 0.8618 ± 0.0132 meV*counts, bounds=[0.0:inf]>]

A ParameterAnalysis can also be used to fit a Parameter to a more simple model. For example, we can fit the area of the Delta function to a polynomial. Here, we only fit a singe Parameter, so the targets is just the name of the Parameter as a string. There’s not really a physical reason to fit this; we’re just doing it to illustrate the method.

Notice that in this case, we have to be clear about the units of our fit function. For the x-axis, the unit is now ‘1/angstrom’ (or ‘1/Å’).

delta_poly = edyn.Polynomial(
    coefficients=[0.0, 0.1],
    name="delta_model",
    y_unit="meV*counts",
    x_unit="1/angstrom",
)
delta_binding = edyn.FitBinding(model=delta_poly, targets="DeltaFunction area")
delta_analysis = edyn.ParameterAnalysis(
    parameters=qe_analysis,
    bindings=delta_binding,
)
delta_analysis.plot()

delta_analysis.fit()
delta_analysis.plot()

Finally, if you need a model that is not (yet) built into EasyDynamics, you can always contact us, and we might implement it. In the meantime, you can define your own expression. It recognizes pi, e, hbar and kb, automatically extracts other parameters, and is unit-aware. You will not need it here, but we show it for completeness. You can optionally give parameters start values and give them units as shown below. You should get the same values as before.

jump_diffusion_width_expression = edyn.ExpressionComponent(
    expression="hbar * D * x**2 / (1 + D * x**2 * tau)",
    name="diffusion_expression",
    parameters={"D": 4.6e-10, "tau": 22.0},
    parameter_units={"D": "m^2/s", "tau": "ps"},
    x_unit="1/angstrom",
    y_unit="meV",
)

jump_diffusion_area = edyn.Polynomial(
    coefficients=[0.0], x_unit="1/angstrom", y_unit="meV * counts", name="scale"
)

width_binding = edyn.FitBinding(
    model=jump_diffusion_width_expression, targets="Lorentzian width"
)

area_binding = edyn.FitBinding(model=jump_diffusion_area, targets="Lorentzian area")

expr_parameter_analysis = edyn.ParameterAnalysis(
    parameters=qe_analysis,
    bindings=[width_binding, area_binding],
)

expr_parameter_analysis.fit()
expr_parameter_analysis.get_all_parameters()

[<Parameter 'D': 4.288e-10 ± 3.179e-11 m^2/s, bounds=[-inf:inf]>,
 <Parameter 'scale_c0': 0.8618 ± 0.0132, bounds=[-inf:inf]>,
 <Parameter 'tau': 29.5564 ± 1.4116 ps, bounds=[-inf:inf]>]

Optional exercise:#

Look up other diffusion models and try fitting them by using the ExpressionComponent.

Step 4: Fit the jump-diffusion model to all the data at once#

The two-step approach to find the diffusion constant and residence time shown above works, but we can do better. Now that we know the quasi-elastic scattering follows a jump-diffusion model, we can fit that model directly to the data, using all Q values simultaneously. In addition to the diffusion, we still describe the elastic incoherent scattering with a DeltaFunction.

We build a new SampleModel that has a DeltaFunction component and a JumpTranslationalDiffusion diffusion model, and new BackgroundModel and InstrumentModel objects.

You can put in the fitted values of scale, diffusion_coefficient and relaxation_time as starting guesses if you want.

delta_function_diff = edyn.DeltaFunction(
    name="DeltaFunction", area=0.15, y_unit="counts"
)

diffusion_model = edyn.JumpTranslationalDiffusion(
    name="Jump Translational Diffusion",
    diffusion_coefficient=2.5e-10,
    relaxation_time=18.0,  # ps
    scale=1.0,
    y_unit="counts",
)

sample_model_diff = edyn.SampleModel(
    components=delta_function_diff,
    diffusion_models=diffusion_model,
    y_unit="counts",
)

Exercise: Complete the creation of the analysis object for the Jump Diffusion#

Follow the steps above: define an InstrumentModel, where you pass the elastic analysis as the resolution model, consider adding a background, etc. Remember to use the right y_unit. Call the new analysis diffusion_analysis.

Solution:

As always, we check the start guess before fitting and iterate on the parameters until we are fairly close.

diffusion_analysis.plot_data_and_model()

Now we fit all the data simultaneously. If this takes a long time (>30 seconds), consider improving your starting guesses.

diffusion_analysis.fit(fit_method="simultaneous")

[FitResults(success=True
   n_pars=19, n_points=530
   chi2=188, reduced_chi2=0.3679
   n_evaluations=121
   iterations=121
   minimizer=LMFit
   message='Fit succeeded.'
   parameters={pParameter_341=0.122, pParameter_144=4.792e-10, pParameter_145=30.85, pParameter_143=0.8514, pParameter_356=0.001661, pParameter_343=0.0826, pParameter_357=0.00137, pParameter_345=0.07256, pParameter_358=0.001563, pParameter_347=0.08009, pParameter_359=0.001589, pParameter_349=0.07586, pParameter_360=0.001844, pParameter_351=0.06069, pParameter_361=0.001926, pParameter_353=0.06257, pParameter_362=0.001966, pParameter_355=0.07081, pParameter_363=0.001524}
 ),
 FitResults(success=True
   n_pars=19, n_points=530
   chi2=256.7, reduced_chi2=0.5023
   n_evaluations=121
   iterations=121
   minimizer=LMFit
   message='Fit succeeded.'
   parameters={pParameter_341=0.122, pParameter_144=4.792e-10, pParameter_145=30.85, pParameter_143=0.8514, pParameter_356=0.001661, pParameter_343=0.0826, pParameter_357=0.00137, pParameter_345=0.07256, pParameter_358=0.001563, pParameter_347=0.08009, pParameter_359=0.001589, pParameter_349=0.07586, pParameter_360=0.001844, pParameter_351=0.06069, pParameter_361=0.001926, pParameter_353=0.06257, pParameter_362=0.001966, pParameter_355=0.07081, pParameter_363=0.001524}
 ),
 FitResults(success=True
   n_pars=19, n_points=530
   chi2=347.1, reduced_chi2=0.6793
   n_evaluations=121
   iterations=121
   minimizer=LMFit
   message='Fit succeeded.'
   parameters={pParameter_341=0.122, pParameter_144=4.792e-10, pParameter_145=30.85, pParameter_143=0.8514, pParameter_356=0.001661, pParameter_343=0.0826, pParameter_357=0.00137, pParameter_345=0.07256, pParameter_358=0.001563, pParameter_347=0.08009, pParameter_359=0.001589, pParameter_349=0.07586, pParameter_360=0.001844, pParameter_351=0.06069, pParameter_361=0.001926, pParameter_353=0.06257, pParameter_362=0.001966, pParameter_355=0.07081, pParameter_363=0.001524}
 ),
 FitResults(success=True
   n_pars=19, n_points=530
   chi2=370.5, reduced_chi2=0.7251
   n_evaluations=121
   iterations=121
   minimizer=LMFit
   message='Fit succeeded.'
   parameters={pParameter_341=0.122, pParameter_144=4.792e-10, pParameter_145=30.85, pParameter_143=0.8514, pParameter_356=0.001661, pParameter_343=0.0826, pParameter_357=0.00137, pParameter_345=0.07256, pParameter_358=0.001563, pParameter_347=0.08009, pParameter_359=0.001589, pParameter_349=0.07586, pParameter_360=0.001844, pParameter_351=0.06069, pParameter_361=0.001926, pParameter_353=0.06257, pParameter_362=0.001966, pParameter_355=0.07081, pParameter_363=0.001524}
 ),
 FitResults(success=True
   n_pars=19, n_points=530
   chi2=426.7, reduced_chi2=0.835
   n_evaluations=121
   iterations=121
   minimizer=LMFit
   message='Fit succeeded.'
   parameters={pParameter_341=0.122, pParameter_144=4.792e-10, pParameter_145=30.85, pParameter_143=0.8514, pParameter_356=0.001661, pParameter_343=0.0826, pParameter_357=0.00137, pParameter_345=0.07256, pParameter_358=0.001563, pParameter_347=0.08009, pParameter_359=0.001589, pParameter_349=0.07586, pParameter_360=0.001844, pParameter_351=0.06069, pParameter_361=0.001926, pParameter_353=0.06257, pParameter_362=0.001966, pParameter_355=0.07081, pParameter_363=0.001524}
 ),
 FitResults(success=True
   n_pars=19, n_points=530
   chi2=422.2, reduced_chi2=0.8262
   n_evaluations=121
   iterations=121
   minimizer=LMFit
   message='Fit succeeded.'
   parameters={pParameter_341=0.122, pParameter_144=4.792e-10, pParameter_145=30.85, pParameter_143=0.8514, pParameter_356=0.001661, pParameter_343=0.0826, pParameter_357=0.00137, pParameter_345=0.07256, pParameter_358=0.001563, pParameter_347=0.08009, pParameter_359=0.001589, pParameter_349=0.07586, pParameter_360=0.001844, pParameter_351=0.06069, pParameter_361=0.001926, pParameter_353=0.06257, pParameter_362=0.001966, pParameter_355=0.07081, pParameter_363=0.001524}
 ),
 FitResults(success=True
   n_pars=19, n_points=530
   chi2=501.2, reduced_chi2=0.9808
   n_evaluations=121
   iterations=121
   minimizer=LMFit
   message='Fit succeeded.'
   parameters={pParameter_341=0.122, pParameter_144=4.792e-10, pParameter_145=30.85, pParameter_143=0.8514, pParameter_356=0.001661, pParameter_343=0.0826, pParameter_357=0.00137, pParameter_345=0.07256, pParameter_358=0.001563, pParameter_347=0.08009, pParameter_359=0.001589, pParameter_349=0.07586, pParameter_360=0.001844, pParameter_351=0.06069, pParameter_361=0.001926, pParameter_353=0.06257, pParameter_362=0.001966, pParameter_355=0.07081, pParameter_363=0.001524}
 ),
 FitResults(success=True
   n_pars=19, n_points=530
   chi2=543.7, reduced_chi2=1.064
   n_evaluations=121
   iterations=121
   minimizer=LMFit
   message='Fit succeeded.'
   parameters={pParameter_341=0.122, pParameter_144=4.792e-10, pParameter_145=30.85, pParameter_143=0.8514, pParameter_356=0.001661, pParameter_343=0.0826, pParameter_357=0.00137, pParameter_345=0.07256, pParameter_358=0.001563, pParameter_347=0.08009, pParameter_359=0.001589, pParameter_349=0.07586, pParameter_360=0.001844, pParameter_351=0.06069, pParameter_361=0.001926, pParameter_353=0.06257, pParameter_362=0.001966, pParameter_355=0.07081, pParameter_363=0.001524}
 )]

diffusion_analysis.plot_data_and_model(plot_residuals=True, autoscale=False)

The diffusion parameters are just a couple of numbers with uncertainties, so instead of plotting them we display them directly. Both the diffusion coefficient \(D\) and the residence time \(\tau\) are now determined by the data.

diffusion_model.get_global_variables()

[<Parameter 'diffusion_coefficient': 4.792e-10 ± 1.418e-11 m^2/s, bounds=[0.0:inf]>,
 <Parameter 'relaxation_time': 30.8458 ± 0.4605 ps, bounds=[0.0:inf]>,
 <Parameter 'scale': 0.8514 ± 0.0047 meV*counts, bounds=[0.0:inf]>]

For reference, here are the parameters from the two-step fit (fitting the Lorentzian widths and areas). Compare the two results.

parameter_analysis.get_all_parameters()

[<Parameter 'diffusion_coefficient': 4.288e-10 ± 3.148e-11 m^2/s, bounds=[0.0:inf]>,
 <Parameter 'relaxation_time': 29.5564 ± 1.4116 ps, bounds=[0.0:inf]>,
 <Parameter 'scale': 0.8618 ± 0.0132 meV*counts, bounds=[0.0:inf]>]

display_quiz(quiz.q9)

Finally, since we know what we put into McStas we can compare our answers to the true values: \(D = 4.6\times10^{-10}\) m\(^2\)/s and \(\tau = 22\) ps. Did you get similar values? Why/why not?

display_quiz(quiz.q10)

The numbers! What do they mean?#

We did not choose any particular material to simualate, just a material containing some static and some mobile hydrogen scattering according to the equations described here. Still, in a real analysis, we might continue as follows.

In water, \(D\approx 2.3 \times 10^{-9}\) m\(^2\)/s, so the diffusion here is about 5 times slower, indicating that the hydrogen is confined. In water, the residence time (i.e. the average time between jumps) is \(\tau \approx 1\) ps, so again, the hydrogen here is much more confined. This could for example be the case for clay or cement, where the strength of the material seems to correlate with how confined the water in it is.

The average jump length can be calculated as \(l \sim \sqrt{6D \tau} \sim 2.5\) \AA{}, which indicates that the hydrogen atoms are jumping between neighbouring sites. If this were a real material, we might know the crystal structure, we might have other measurements, and might therefore be able to draw some further conclusions about how the material functions and why.

The scale does not carry any particular meaning - it’s a complex combination of the neutron flux of the instrument, the instrument geometry and the sample composition and geometry.

This is the end of this analysis. EasyDynamics is work in progress, and we would very much like to hear your feedback - positive and negative. If you would like to fit your own data with EasyDynamics then please reach out. We would be delighted to help you get started.

Please write to henrik.jacobsen@ess.eu.

Fitting QENS data

Contents

Fitting QENS data#

Load and prepare the data#

Step 1: Determine the resolution from the elastic sample#

Normalise to the resolution area#

Step 2: Fit the quasi-elastic sample at each Q#

Exercise: create a new SampleModel for the quasi-elastic sample#

Exercise: collect the data, SampleModel and InstrumentModel in an Analysis object#

Step 3: Fit a jump-diffusion model to the widths#

Optional exercise:#

Step 4: Fit the jump-diffusion model to all the data at once#

Exercise: Complete the creation of the analysis object for the Jump Diffusion#

The numbers! What do they mean?#