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 we model it as a delta function. We do this by creating a SampleModel and adding a DeltaFunction component to it. The component acts as a template and gets copied to every Q when we attach the SampleModel to our Analysis object. Let's create the SampleModel.
We do not give the DeltaFunction a center value. In this case, the center will be fixed at 0 energy transfer. We set the start value of the area to 1.
delta_function = sm.DeltaFunction(name='DeltaFunction', area=1)
sample_model = sm.SampleModel(components=delta_function)
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.
We fix the area of the resolution to have value 1. If we did not do this, we would fit both the area of the delta function and of the resolution Gaussian, and the fit would never converge.
We finally insert the components in a ResolutionModel
resolution_components = sm.ComponentCollection()
res_gauss = sm.Gaussian(width=0.1, area=1, name='Res. Gauss')
res_gauss.area.fixed = True
resolution_components.append_component(res_gauss)
resolution_model = sm.ResolutionModel(components=resolution_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 did for the ResolutionModel
background_model = sm.BackgroundModel(components=sm.Polynomial(coefficients=[0.001]))
We combine the resolution abd background model into 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(
resolution_model=resolution_model,
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=sample_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. ])
- DeltaFunction 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]) - DeltaFunction 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.]) - 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)float64meV1.0, 1.0, ..., 1.0, 1.0σ = 0.0, 0.0, ..., 0.0, 0.0
Values:
array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])
Variances (σ²):
array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]) - 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.63407023e-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.33891037e-05, 1.37163920e-05, 1.12418484e-05, 1.18955355e-05, 1.14728274e-05, 9.35171956e-06, 1.83656779e-05, 1.03814854e-05, 1.29102630e-05, 1.47729559e-05, 1.26723273e-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=['DeltaFunction 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, giving it the resolution model determined in the vanadium analysis. We further fix all parameters in the resolution model and normalize it.
instrument_model = sm.InstrumentModel(
background_model=background_model,
resolution_model=vanadium_analysis.instrument_model.resolution_model,
)
instrument_model.resolution_model.fix_all_parameters()
instrument_model.normalize_resolution()
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 'diffusion_coefficient': 1.152e-08 ± 1.248e-09 m^2/s, bounds=[0.0:inf]>, <Parameter 'scale': 0.6062 ± 0.0426 meV, 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.instrument_model.resolution_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_all_parameters()
[<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_parameters()
[<Parameter 'diffusion_coefficient': 1.152e-08 ± 1.248e-09 m^2/s, bounds=[0.0:inf]>, <Parameter 'scale': 0.6062 ± 0.0426 meV, bounds=[0.0:inf]>]