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. Finally, we show how to fit directly to the diffusion model.
In the near future, it will be possible to fit the width and area of the Lorentzian to the diffusion model as well.
# Imports
import pooch
from easydynamics.analysis.analysis import Analysis
from easydynamics.experiment import Experiment
from easydynamics.sample_model import BrownianTranslationalDiffusion
from easydynamics.sample_model import ComponentCollection
from easydynamics.sample_model import DeltaFunction
from easydynamics.sample_model import Gaussian
from easydynamics.sample_model import Lorentzian
from easydynamics.sample_model import Polynomial
from easydynamics.sample_model.background_model import BackgroundModel
from easydynamics.sample_model.instrument_model import InstrumentModel
from easydynamics.sample_model.resolution_model import ResolutionModel
from easydynamics.sample_model.sample_model import SampleModel
# 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 = Experiment('Vanadium')
file_path = pooch.retrieve(
url='https://github.com/easyscience/dynamics-lib/raw/refs/heads/docs/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.
vanadium_experiment.plot_data(slicer=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 = DeltaFunction(display_name='DeltaFunction', area=1)
sample_model = 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 = ComponentCollection()
res_gauss = Gaussian(width=0.1, area=1, display_name='Res. Gauss')
res_gauss.area.fixed = True
resolution_components.append_component(res_gauss)
resolution_model = 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 = BackgroundModel(components=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 = 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 = 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 = Experiment('Diffusion')
file_path = pooch.retrieve(
url='https://github.com/easyscience/dynamics-lib/raw/refs/heads/docs/docs/docs/tutorials/data/diffusion_data_example.h5',
)
diffusion_experiment.load_hdf5(filename=file_path)
diffusion_experiment.plot_data(slicer=True)
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 = DeltaFunction(display_name='DeltaFunction', area=0.2)
lorentzian = Lorentzian(display_name='Lorentzian', area=0.5, width=0.3)
component_collection = ComponentCollection(
components=[delta_function, lorentzian],
)
sample_model = SampleModel(
components=component_collection,
)
background_model = BackgroundModel(components=Polynomial(coefficients=[0.001]))
We also create a new instrument_model and attach it to our analysis. We want to use the resolutin that we determined from the vanadium data. Ideally, we'd just attach it to the instrument model like before. However, this is currently not possible due to a small issue in EasyDynamics that will be fixed asap. For now, we need to hack it a bit.
instrument_model = InstrumentModel(
background_model=background_model,
)
diffusion_analysis = Analysis(
display_name='Diffusion Analysis',
experiment=diffusion_experiment,
sample_model=sample_model,
instrument_model=instrument_model,
)
# We need to hack in the resolution model from the vanadium analysis,
# since the setters and getters overwrite the model. This will be fixed
# asap.
diffusion_analysis.instrument_model._resolution_model = (
vanadium_analysis.instrument_model.resolution_model
)
We don't want to fit our resolution anymore, so we fix all the parameters in it.
# We fix all parameters of the resolution model.
diffusion_analysis.instrument_model.resolution_model.fix_all_parameters()
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. In later versions of EasyDynamics it will be possible to fit them to e.g. a DiffusionModel.
# Let us look at the most interesting fit parameters
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.
In addition to this diffusion model, there is still 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 = DeltaFunction(display_name='DeltaFunction', area=0.2)
component_collection = ComponentCollection(
components=[delta_function],
)
diffusion_model = BrownianTranslationalDiffusion(
display_name='Brownian Translational Diffusion', diffusion_coefficient=2.4e-9, scale=0.5
)
sample_model = SampleModel(
components=component_collection,
diffusion_models=diffusion_model,
)
background_model = BackgroundModel(components=Polynomial(coefficients=[0.001]))
instrument_model = InstrumentModel(
background_model=background_model,
)
We attach all our models to a new Analysis object, again being mindful that we need to hack in the resolution. I promise this will be fixed soon!
diffusion_model_analysis = Analysis(
display_name='Diffusion Full Analysis',
experiment=diffusion_experiment,
sample_model=sample_model,
instrument_model=instrument_model,
)
# We again need to hack in the resolution model from the vanadium
# analysis, since the setters and getters overwrite the model. This will
# be fixed asap.
diffusion_model_analysis.instrument_model._resolution_model = (
vanadium_analysis.instrument_model.resolution_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()
It does not make sense to plot the diffusion parameters, 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, bounds=[0.0:inf]>]