Brownian Diffusion¶
We here show how to set up an Analysis object and use it to first fit an artificial vanadium measurement to obtain the resolution. Next, we use the fitted resolution to fit an artificial measurement of a model with diffusion and some elastic scattering.
We extract and plot the relevant parameters and fit them to a diffusion model. Finally, we show how to fit all the data simultaneously to the diffusion model.
# Imports
import pooch
import easydynamics as edyn
import easydynamics.sample_model as sm
# Make the plots interactive
%matplotlib widget
We first create an Experiment object to contain the data. The data must either be a hdf5 file or a scipp.DataArray; in both cases it must have coordinates Q and energy. We here use Pooch to download an example vanadium data set.
The data can be rebinned if needed, but we will show how to do that in a different tutorial.
# Load the vanadium data
vanadium_experiment = edyn.Experiment('Vanadium')
file_path = pooch.retrieve(
url='https://github.com/easyscience/dynamics-lib/raw/refs/heads/master/docs/docs/tutorials/data/vanadium_data_example.h5',
known_hash='16cc1b327c303feeb88fb9dda5390dc4880b62396b1793f98c6fef0b27c7b873',
)
vanadium_experiment.load_hdf5(filename=file_path)
We can visualize the data in multiple ways, relying on plopp: https://scipp.github.io/plopp/
We here show two ways to look at the data: as a 2d colormap with intensity as function of Q and energy, and as a slicer with intensity as function of energy for various Q.
If you want $Q$ on the x axis, then set transpose_axes=True
vanadium_experiment.plot_data(slicer=False, transpose_axes=False)
vanadium_experiment.plot_data(slicer=True)
We now want to fit the vanadium data to determine our resolution. The scattering from vanadium is almost exclusively incoherent elastic, so it can be modelled as a delta function. There are two ways to do this. The first is to add a DeltaFunction component to a SampleModel and next create a ResolutionModel. The second method, which we here shall use, is to fit the data to a SampleModel and then convert it to a ResolutionModel.
We now want to define our resolution function. We will here model it as a Gaussian. We create a ComponentCollection and append the Gaussian to it. We can add as many components to our resolution as we like; sometimes you need several Gaussians and other functions to accurately describe the resolution.
vanadium_components = sm.ComponentCollection()
res_gauss = sm.Gaussian(width=0.1, area=1, name='Res. Gauss')
vanadium_components.append_component(res_gauss)
vanadium_model = sm.SampleModel(components=vanadium_components)
The background intensity was not 0, so we also create a background model. We use a Polynomial with a single coefficient, i.e. a flat background. We here show how to create the BackgroundModel and add the background in a single line. We could of course also add it like we did for the SampleModel or first create a ComponentCollection like we just did for the resolution.
background_model = sm.BackgroundModel(components=sm.Polynomial(coefficients=[0.001]))
We add the background model to an InstrumentModel. This model also contains a fittable energy offset to account for instrument misalignment. All components are centered at this energy offset.
instrument_model = sm.InstrumentModel(
background_model=background_model,
)
We are now ready to collect everything in an analysis object. We give it a display name, the experiment, the sample model and the instrument model. It will then automatically generate a model for each Q using the templates given in the SampleModel, ResolutionModel and BackgroundModel.
vanadium_analysis = edyn.Analysis(
display_name='Vanadium Full Analysis',
experiment=vanadium_experiment,
sample_model=vanadium_model,
instrument_model=instrument_model,
)
Let us first fit a single Q index and plot the data and model to see how it looks. For this, we use the independent fit method and choose an arbitrary Q index
fit_result_independent_single_Q = vanadium_analysis.fit(fit_method='independent', Q_index=5)
vanadium_analysis.plot_data_and_model(Q_index=5)
The fit looks good, so let us fit all Q indices independently and plot the results.
fit_result_independent_all_Q = vanadium_analysis.fit(fit_method='independent')
vanadium_analysis.plot_data_and_model()
It can be nice to inspect the fit parameters, and sometimes continue working with them. To do this, we can convert them to a scipp dataset.
We can also plot the parameters as a function of Q using the plot_parameters method.
# Inspect the Parameters as a scipp Dataset
vanadium_analysis.parameters_to_dataset()
- Q: 16
- Q(Q)float641/Å0.1, 0.227, ..., 1.873, 2.0
Values:
array([0.1 , 0.22666667, 0.35333333, 0.48 , 0.60666667, 0.73333333, 0.86 , 0.98666667, 1.11333333, 1.24 , 1.36666667, 1.49333333, 1.62 , 1.74666667, 1.87333333, 2. ])
- Polynomial_c0(Q)float64𝟙0.099, 0.095, ..., 0.098, 0.101σ = 0.002, 0.002, ..., 0.002, 0.002
Values:
array([0.09942061, 0.09482336, 0.10250515, 0.09830358, 0.10212681, 0.1023704 , 0.09957118, 0.1020655 , 0.09643967, 0.10180256, 0.10191825, 0.09198844, 0.09427405, 0.09806923, 0.09764543, 0.10051938])
Variances (σ²):
array([5.53666458e-06, 5.69780006e-06, 4.81032728e-06, 4.89897712e-06, 4.99355959e-06, 4.08027374e-06, 7.91272766e-06, 4.56234563e-06, 5.13312615e-06, 6.29261176e-06, 5.53177991e-06, 7.73519893e-06, 6.33375780e-06, 5.95523123e-06, 4.94577573e-06, 4.98191241e-06]) - Res. Gauss area(Q)float64meV0.522, 0.529, ..., 0.538, 0.533σ = 0.018, 0.018, ..., 0.017, 0.017
Values:
array([0.52209408, 0.5285229 , 0.51741732, 0.53239306, 0.52764215, 0.5299243 , 0.52765768, 0.51477418, 0.52913757, 0.53258438, 0.5324246 , 0.53238945, 0.52271039, 0.52427006, 0.53758902, 0.53293296])
Variances (σ²):
array([0.00030823, 0.00033506, 0.00025861, 0.00028047, 0.00027407, 0.0002244 , 0.00044408, 0.00024497, 0.00029737, 0.00034951, 0.00030663, 0.00047018, 0.00037023, 0.00033697, 0.00028746, 0.00028007]) - Res. Gauss center(Q)float64meV0.0, 0.0, ..., 0.0, 0.0σ = 0.0, 0.0, ..., 0.0, 0.0
Values:
array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])
Variances (σ²):
array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]) - Res. Gauss width(Q)float64meV0.102, 0.100, ..., 0.103, 0.102σ = 0.003, 0.003, ..., 0.003, 0.003
Values:
array([0.10222849, 0.09980219, 0.10094395, 0.10308798, 0.10108169, 0.10126272, 0.10063122, 0.09902007, 0.10372129, 0.10254205, 0.10142844, 0.10315841, 0.10076283, 0.10044881, 0.10346177, 0.10187953])
Variances (σ²):
array([9.25121501e-06, 8.81462902e-06, 7.61444384e-06, 8.02505303e-06, 7.63407022e-06, 6.32632931e-06, 1.23202070e-05, 6.87361362e-06, 8.78013502e-06, 9.96635847e-06, 8.46568221e-06, 1.33366848e-05, 1.06131452e-05, 9.28744245e-06, 8.02343194e-06, 7.89725747e-06]) - energy_offset(Q)float64meV0.001, -0.001, ..., -0.001, 0.001σ = 0.004, 0.004, ..., 0.003, 0.003
Values:
array([ 0.00135435, -0.00052522, 0.00093028, -0.00192069, -0.00366331, 0.00128853, 0.000512 , 0.00090986, -0.00100558, 0.00036124, 0.00054189, -0.00216345, -0.00137654, 0.00129853, -0.001399 , 0.00055349])
Variances (σ²):
array([1.33891038e-05, 1.37163919e-05, 1.12418484e-05, 1.18955355e-05, 1.14728274e-05, 9.35171954e-06, 1.83656779e-05, 1.03814854e-05, 1.29102630e-05, 1.47729560e-05, 1.26723274e-05, 1.99321937e-05, 1.56095496e-05, 1.41010062e-05, 1.20573093e-05, 1.16915411e-05])
# Plot some of fitted parameters as a function of Q
vanadium_analysis.plot_parameters(names=['Res. Gauss area'])
vanadium_analysis.plot_parameters(names=['Res. Gauss width'])
vanadium_analysis.plot_parameters(names=['energy_offset'])
We are now happy with our resolution function and can start looking at the data we want to fit. We first load and inspect the data in the same way as before.
diffusion_experiment = edyn.Experiment('Diffusion')
file_path = pooch.retrieve(
url='https://github.com/easyscience/dynamics-lib/raw/refs/heads/master/docs/docs/tutorials/data/diffusion_data_example.h5',
known_hash='5fe846b19aacbda8b8b936eb2e5310d025dc56c25b0b353521e7d6b921f229ab',
)
diffusion_experiment.load_hdf5(filename=file_path)
diffusion_experiment.plot_data(slicer=True, ymax=4)
The data seems to have a sharp elastic peak, a quasielastic peak and a non-zero background. We set up the corresponding SampleModel and BackgroundModel just like before.
delta_function = sm.DeltaFunction(name='DeltaFunction', area=0.2)
lorentzian = sm.Lorentzian(name='Lorentzian', area=0.5, width=0.3)
component_collection = sm.ComponentCollection(
components=[delta_function, lorentzian],
)
sample_model = sm.SampleModel(
components=component_collection,
)
background_model = sm.BackgroundModel(components=sm.Polynomial(coefficients=[0.001]))
We also create a new instrument_model and attach it to our analysis. We now give it the sample model from the vanadium fit. All parameters are automatically fixed and the resolution model is normalized to have area 1.
instrument_model = sm.InstrumentModel(
background_model=background_model,
resolution_model=vanadium_analysis.sample_model,
)
diffusion_analysis = edyn.Analysis(
display_name='Diffusion Analysis',
experiment=diffusion_experiment,
sample_model=sample_model,
instrument_model=instrument_model,
)
We don't want to fit our resolution anymore, so we fix all the parameters in it.
Analysis handles the convolution of the sample_model with the resolution_model. The calculation is analytical when possible, and otherwise numerical. For numerical convolution, it will give warnings if the result might be inaccurate due to various numerical errors. In that case, two settings can be varied to improve the result.
upsample_factor improves accuracy for narrow signals. The default is 5.
extension_factor improves accuracy for broad signals. Here, the default is 0.2.
It is furthermore possible to toggle whether detailed balance correction is normalized or not (see next tutorial)
diffusion_analysis.convolution_settings.upsample_factor = 6
diffusion_analysis.convolution_settings.extension_factor = 0.2
diffusion_analysis.convolution_settings.normalize_detailed_balance = True
Before we start fitting it is a good idea to check how good our start guesses are. We do this by plotting the data and the model using plot_data_and_model:
diffusion_analysis.plot_data_and_model()
The start guesses are not perfect, but they look good enough to start fitting. Let's give it a try!
diffusion_analysis.fit(fit_method='independent')
diffusion_analysis.plot_data_and_model()
The fit looks good, so now we want to look at the most interesting fit parameters: the width and area of the Lorentzian.
diffusion_analysis.plot_parameters(names=['Lorentzian width', 'Lorentzian area'])
There is a clear trend: the area is more or less constant, while the width seems to increase with Q^2. We therefore try and fit a Brownian translational diffusion model to the data. In this model, the scattering is given by
$$
I(Q,E) = S \frac{\Gamma(Q)}{\Gamma(Q)^2 + E^2},
$$
where $\Gamma(Q) = D Q^2$ and $D$ is the diffusion coefficient. $S$ is an overall scale.
To fit the Brownian translational diffusion model to the data, we use the ParameterAnalysis. This time, we wish to fit the Lorentzian, and we wish to fit both its area (scale, $S$) and width to the diffusion model. We do this by creating a FitBinding, saying we want to fit both the area an width of the component called Lorentzian
brownian_diffusion_model = sm.BrownianTranslationalDiffusion(
name='Brownian Translational Diffusion', diffusion_coefficient=2.4e-9, scale=0.5
)
binding = edyn.FitBinding(
parameter_name='Lorentzian',
model=brownian_diffusion_model,
modes=['area', 'width'],
)
parameter_analysis = edyn.ParameterAnalysis(
parameters=diffusion_analysis,
bindings=[binding],
)
We first plot the start guess to see if it's reasonable.
parameter_analysis.plot()
parameter_analysis.plot(names=['Polynomial_c0'])
This looks pretty good. Let's fit and plot again:
parameter_analysis.fit()
parameter_analysis.plot()
And we can get the parameters of the diffusion model:
parameter_analysis.get_all_parameters()
[<Parameter 'scale': 0.6062 ± 0.0426 meV, bounds=[0.0:inf]>, <Parameter 'diffusion_coefficient': 1.152e-08 ± 1.248e-09 m^2/s, bounds=[0.0:inf]>]
This is a good result, but we can do better. Now that we know that the quasielastic scattering is well described by a model of diffusion, we can fit this model directly to the data, fitting all $Q$ simultaneously.
In addition to this diffusion model, we will still fit the elastic incoherent scattering.
We create a new SampleModel which as a DeltaFunction component for the elastic incoherent scattering and a BrownianTranslationalDiffusion diffusion model to describe the rest. We also create a new BackgroundModel and InstrumentModel.
delta_function = sm.DeltaFunction(name='DeltaFunction', area=0.2)
component_collection = sm.ComponentCollection(
components=[delta_function],
)
diffusion_model = sm.BrownianTranslationalDiffusion(
name='Brownian Translational Diffusion', diffusion_coefficient=2.4e-9, scale=0.5
)
sample_model = sm.SampleModel(
components=component_collection,
diffusion_models=diffusion_model,
)
background_model = sm.BackgroundModel(components=sm.Polynomial(coefficients=[0.001]))
instrument_model = sm.InstrumentModel(
background_model=background_model,
resolution_model=vanadium_analysis.sample_model,
)
We attach all our models to a new Analysis object.
diffusion_model_analysis = edyn.Analysis(
display_name='Diffusion Full Analysis',
experiment=diffusion_experiment,
sample_model=sample_model,
instrument_model=instrument_model,
)
diffusion_model_analysis.instrument_model.resolution_model.fix_all_parameters()
As always, we first check the start parameters before we fit
diffusion_model_analysis.plot_data_and_model()
We can now fit all the data simultaneously to our model. It looks good!
diffusion_model_analysis.fit(fit_method='simultaneous')
diffusion_model_analysis.plot_data_and_model(ymax=2)
It does not make sense to plot the diffusion parameters, as they are just two numbers with uncertainties, but we can display them (with uncertainties) like this.
diffusion_model.get_global_variables()
[<Parameter 'diffusion_coefficient': 1.126e-08 ± 9.799e-11 m^2/s, bounds=[0.0:inf]>, <Parameter 'scale': 0.6938 ± 0.0043 meV, bounds=[0.0:inf]>]
For reference, here are the diffusion parameters found from fitting the width and scale of the fitted Lorentzians. Notice that the error bars when fitting everything simultaneously are a lot lower than the two-step fit.
parameter_analysis.get_all_variables()
[<Parameter 'scale': 0.6062 ± 0.0426 meV, bounds=[0.0:inf]>, <Parameter 'diffusion_coefficient': 1.152e-08 ± 1.248e-09 m^2/s, bounds=[0.0:inf]>]