.. include:: substitutions.txt .. currentmodule:: jscatter Fitting 1D, 2D 3D ... --------------------- 1D fits with attributes ^^^^^^^^^^^^^^^^^^^^^^^ Some Sinusoidal fits with different kinds to use dataArray attributes. **We use dataArray for fitting**. .. literalinclude:: ../../examples/example_BG_SinusoidalFitting.py :language: python :lines: 1-28 .. image:: ../../examples/SinusoidalFit0.png :align: center :height: 300px :alt: SinusoidalFit 2D fit with attributes ^^^^^^^^^^^^^^^^^^^^^^ A **2D** fit using the attribute *B* stored in the dataArray of a dataList as second dimension. **We use dataList for fitting**. This might be extended to several attributes allowing **multidimensional fitting**. See also :ref:`Simple diffusion fit of not so simple diffusion case` .. literalinclude:: ../../examples/example_BG_SinusoidalFitting.py :language: python :lines: 29- .. image:: ../../examples/SinusoidalFit.png :align: center :height: 300px :alt: SinusoidalFit 2D fit with attributes and linked parameter ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ A **2D** fit as above but the fit parameter is linked to the attribute ``B`` that apears several times with the same value. In experimental data this might be a repeated measurement with a different intensity for the same process. We get only as many fit parameters as we have unique values for ``B`` (here 3 different values). Important is to ``link`` by adding 'B' in ``{..., 'p':[0, 'B'], ...}`` to the list of starting values. .. literalinclude:: ../../examples/example_BG_SinusoidalFitting_linkedParameter.py :language: python :lines: 2- .. image:: ../../examples/SinusoidalFit_linkedParameter.png :align: center :height: 300px :alt: SinusoidalFit 2D fitting using .X, .Z, .W ^^^^^^^^^^^^^^^^^^^^^^^^^^^ Unlike the previous we use here data with two dimensions in .X,.Z coordinates (optional .W for 3D). **We use dataArray for fitting**. Additional one could use again attributes to increase dimesion. This mainly depends on the data. Another example is shown in `Fitting the 2D scattering of a lattice`_. .. literalinclude:: ../../examples/example_BG_2DFitting.py :language: python .. image:: ../../examples/Sinusoidal2D.png :align: center :height: 300px :alt: Sinusoidal3D Simple diffusion fit of not so simple diffusion case ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Here the long part with description from the first example. This is the diffusion of a protein in solution. This is NOT constant as for Einstein diffusion. These simulated data are similar to data measured by Neutron Spinecho Spectroscopy, which measures on the length scale of the protein and therefore also rotational diffusion contributes to the signal. At low wavevectors additional the influence of the structure factor leads to an upturn, which is neglected in the simulated data. To include the correction :math:`D_T(q)=D_{T0} H(q)/S(q)` look at :func:`~.structurefactor.fluid.hydrodynamicFunct`. For details see this tutorial review `Biehl et al. Softmatter 7,1299 (2011) `_ A key point for a simultaneous combined fit is that each dataArray has an individual 'q' which is used in the fit model as 'wavevector'. .. literalinclude:: ../../examples/example_BG_simple_diffusion.py :language: python :end-before: # This is the basis of the simulated data above .. image:: ../../examples/DiffusionFit.jpg :align: center :height: 300px :alt: Picture about diffusion fit .. literalinclude:: ../../examples/example_BG_simple_diffusion.py :language: python :start-after: p.save('DiffusionFit.png') .. image:: ../../examples/effectiveDiffusion.jpg :align: center :height: 300px :alt: diffusion fit result A long example for diffusion and how to analyze step by step ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ This is a long example to show possibilities. A main feature of the fit is that we can change from a constant fit parameters to a parameter dependent one by simply changing A to [A]. .. literalinclude:: ../../examples/example_BG_fit_diffusion.py :language: python :lines: 2- Fitting a multiShellcylinder in various ways ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. literalinclude:: ../../examples/example_BG_fit_multicylinder.py :language: python :lines: 2- Fitting the 2D scattering of a lattice ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ This example shows how to fit a 2D scattering pattern as e.g. measured by small angle scattering. As a first step one should fit radial averaged data to limit the lattice constants to reasonable ranges and deduce reasonable values for background and other parameters. Because of the topology of the problem with a large plateau and small minima most fit algorithm (using a gradient) will fail. One has to take care that the first estimate of the starting parameters (see above) is already close to a very good guess or one may use a method to find a global minimum as *differential_evolution* or brute force scanning. This requires some time. In the following we limit the model to a few parameters. One needs basically to include more as e.g. the orientation of the crystal in 3D and more parameters influencing the measurement. Prior knowledge about the experiment as e.g. a preferred orientation during a shear experiment is helpful information. Another possibility is to normalize the data e.g. to peak maximum=1 with high q data also fixed. As a good practice it is useful to first fit 1D data to fix some parameters and add e.g. orientation in a second step with 2D fitting. .. literalinclude:: ../../examples/example_BG_2DFitting_CrystalPeaks.py :language: python :lines: 1- .. image:: ../../examples/2Dhex_org.png :align: center :width: 300px :alt: Ewald2 .. image:: ../../examples/2Dhex_fit.png :align: center :width: 300px :alt: Ewald2 Using cloudscattering as fit model ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ At the end a complex shaped object: A cube decorated with spheres of different scattering length. .. literalinclude:: ../../examples/example_BG_cloudscattering.py :lines: 1- .. image:: ../../examples/cubeWithSpheres.png :align: center :height: 300px :alt: cubeWithSpheres Bayesian inference for fitting ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ method='bayes' uses Bayesian inference for modeling and the MCMC algorithms for sampling (we use `emcee `_). Requires a larger amount of function evaluations but *returns errors* from Bayesian statistical analysis and allows inclusion of prior knowledge about parameters and restrictions. For Bayesian analysis the log probability :math:`\ln\,p(y\,|\,x,\sigma,a) + ln(p(a_j))` is maximized with the log likelihood :math:`ln (p(y\,|\,x,\sigma,a))` and the prior :math:`p(a_j)`. Asuming Gaussian statistics we use .. math:: \ln\,p(y\,|\,x,\sigma,a) = -\frac{1}{2} \sum_i \left[\frac{(Y_i-f(X_i,a_1,..a_p))^2}{\sigma_i^2} + \ln \left ( p(a_j) \right )\right] By default an uniform (so-called "uninformative") prior is used .. math:: log(p(a_j)) = \left\{ \begin{array}{ll} 0 & \mbox{if $a_j$ in limits};\\ -inf & \mbox{otherwise}.\end{array} \right. which can be changed using the parameter ``ln_prior`` to be more informative. Assuming for parameter :math:`a_j` a Gaussian distribution of deviations from a mean with width sigma .. math:: \ln \left ( p(a_j) \right ) = - \lambda \sum_i w_j^2 where for :math:`\lambda=0` the prior vanishes. :math:`\lambda` quantifies of by how much we believe that e.g. :math:`w_j = a_j-a_{mean}` should be close to zero and can be choosen as :math:`\lambda=n 0.5/\sigma^2` assuming that :math:`a_j` shoud be within a Gaussian of width sigma, maybe determined from a different measurement. The scale :math:`n` scales the Gaussian that the weight is comparable to the log likelihood. Using :math:`\lambda=0.5/\sigma^2` weights the constrain equal to one datapoint. In general the :math:`p(a_j)` depend only on fit parameters :math:`a_j` and not on :math:`Y_i` or :math:`f(X_i,...)`. To allow more complex constraints, e.g. fitting a protein structure with the constrain that the backbone is not overstretched, the modelValues can be accessed if *log_prior* uses the keyword 'modelValues'. Also the weigthed squared errors :math:`\chi^2_i=\frac{[Y_i-f(X_i,a_1,..a_p)]^2}{\sigma_i^2}` can be accessed using the keyword 'x2i' to regularise in some way different meassurements according to their errors. The example shows how some knowledge about parameters A and D that can be used for the prior. We assume here that the amplitude A is distributed around 1 within a Gaussian distribution and D around 0.09. The evaluation might take somewhat longer for bayes (here 100s on my desktop). .. literalinclude:: ../../examples/example_BG_fit_bayes.py :lines: 1- .. image:: ../../examples/images/bayescorner_withprior.jpg :align: center :height: 300px :alt: bayescorner_withprior If we neglect the normalisation terms -ln_prior is the regularisation term used below. Regularisation for fitting ^^^^^^^^^^^^^^^^^^^^^^^^^^ Regularisation in least squares offers the posibility to use prior knowledge similar to 'bayes' to constrain the solution. Different to 'bayes' it uses the faster :math:`\chi^2 minimization` algorithms with the additional constrain :math:`\propto ln(p(w))`. For a Gaussian prior :math:`p(w)` corresponds to `ridge regression `_ and for a Laplace prior to `Lasso regression `_ . You can define whatever you want (some things are more justified than others) as long as there is a minimum and :math:`ln(p(w))>0`. For :math:`\chi^2 minimization` we minimise :math:`\chi_{red}^2 - ln(p(a_j))`. For a Gaussian prior (`Ridge regression `_) this is .. math:: \chi^2 = \frac{1}{n-p} \sum_i^n \frac{[Y_i-f(X_i,a_1,..a_p)]^2}{\sigma_i^2} + \lambda \sum_j w_j^2 where for :math:`\lambda=0` the conventional :math:`\chi_{red}^2 minimization` is retrieved. :math:`\lambda` quantifies of by how much we believe that e.g. :math:`w_i = A-A_{mean}` should be close to zero and can be choosen as :math:`\lambda=\frac{1}{n}0.5/\sigma^2` similar to the above ln_prior assuming that A shoud be within a Gaussian of width sigma, maybe determined from a different measurement. Selects the weight you want to apply : - Using :math:`\lambda=0.5/\sigma^2` for a single variable in Gaussian distribution would be weighted like all measurement values in :math:`\chi_{red}^2` together. - A divisor 1/n scales the Gaussian that the weight is comparable to one data point. - Weighing with :math:`1/n_{j}` weight each :math:`w_j` like :math:`\chi_{red}^2`. - If the :math:`\lambda \sum_j w_j^2` are to strong weigthed, this will dominate the result. - For the below example with i5 containing 16 Q values with about 20 points a weight equal to 20 points corresponding to one Q might be reasonable. In general the :math:`p(a_j)` depend only on fit parameters :math:`a_j` and not on :math:`Y_i` or :math:`f(X_i,...)`. The *best* :math:`\lambda` may be selected from a series of tries as the one with best fit results. To allow more complex constraints the modelValues can be accessed if *log_prior* uses the keyword ``modelValues``. The weigthed squared errors :math:`\chi^2_i=\frac{[Y_i-f(X_i,a_1,..a_p)]^2}{\sigma_i^2}` can be accessed as ``modelValues.chi2i`` and degrees of freedom as ``modelValues.dof``. Anything calculated in the model can be accessed in ``modelValues``. This can be used to constrain by something calculated in the model or to weight e.g. different meassurements according to their errors, number of datapoints or relative contribution to :math:`\chi^2`. The example shows how some knowledge about parameters A and D can be used for the regularisation. We assume here that the amplitude A is distributed around 1 within a Gaussian distribution and D around 0.09. !!! The example is an example for to strong constraints. Its more how to do it and here without advantage. .. literalinclude:: ../../examples/example_BG_fit_Regularisation.py :lines: 1- Iterative reweighting of multiple measurements ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ In a combined fit of :math:`K` datasets of different types, a naive summation of :math:`\chi^2` contributions - assuming identical and independent distribution of deviations - causes large datasets or dataset with systematically smaller error to dominate the minimization. Reweighting corrects this by assigning each dataset a weight :math:`w_k` like :math:`\chi^2= \sum_k w_k\chi^2_k`. Application : As an example we might think about a NSE measurement with good defined error bars and a dynamic light scattering (DLS) measurement which typically gives no error at all for the correlation function. Here we can estimate an error for DLS and refine this estimate. Another application might be to identify dataset that need larger reweighting because of suspicious data or error estimates, e.g a large number of zero counts in bins or wrong counting during the measurement, and identify outliers. Using different weights was introduced by [1]_ and others as `Iteratively Reweighted Least-Squares `_ (IRLS). IRLS iteratively calculates weights of individual data points and simultaneously tries to find model parameters that fit the bulk of the data using a least-squares approach, to produce estimates robust against outliers. Reweighting of datasets from different measurements was used by Larsen and Pedersen for combined fit of SAXS and SANS measyrements [2]_, [3]_. Here we describe a *iterative* `Expectation-Maximization (EM) algorithm `_ used with the command :py:func:`~jscatter.dataarray.dataList.regularizedEM` as an *iterative reweighting* of different mesurements like SAXS/SANS with different contrast or combining masurements of different instruments like neutron time of flight, backscattering and neutron spin echo as combined in [4]_ or [5]_. *Iterative reweigthing* avoids time consuming Bayesian inference by faster iteration using standard fit routines. This is in particular useful if the model evaluation takes some time. .. collapse:: Excurs: Bayesian Analysis with precision hyperparameter **Bayesian Analysis** In the Bayesian framework the combined fit is formulated as a posterior distribution over the model parameters :math:`\theta` and per-dataset precision hyperparameters :math:`\alpha_k`: .. math:: P(\theta, \{\alpha_k\} \mid \mathrm{data}) \propto P(\theta) \prod_{k=1}^{K} \mathcal{L}_k(\theta, \alpha_k)\, P(\alpha_k), where the Gaussian likelihood for dataset :math:`k` is .. math:: \ln \mathcal{L}_k(\theta, \alpha_k) = \frac{N_k}{2}\ln\alpha_k - \frac{\alpha_k}{2}\,\chi^2_k(\theta), and :math:`\alpha_k` rescales the trustworthiness of dataset :math:`k` — large :math:`\alpha_k` means the dataset is weighted heavily, small :math:`\alpha_k` means it is down-weighted. Without rescaling the likelihood is only :math:`\ln \mathcal{L}_k(\theta)=-\frac{1}{2}\chi^2_k(\theta)`. The conjugate prior for :math:`\alpha_k` - our prior knowledge - is the Gamma distribution .. math:: P(\alpha_k) = \mathrm{Gamma}(\alpha_k \mid a_k, b_k) \propto \alpha_k^{a_k-1}\,e^{-b_k \alpha_k}, whose kernel matches the likelihood in :math:`\alpha_k` exactly, so the posterior on :math:`\alpha_k` is also Gamma with updated parameters: .. math:: P(\alpha_k \mid \theta, \mathrm{data}) = \mathrm{Gamma}\!\left(a_k + \tfrac{N_k}{2},\; b_k + \tfrac{\chi^2_k(\theta)}{2}\right), with posterior mean .. math:: \langle\alpha_k\rangle = \frac{a_k + N_k/2}{b_k + \chi^2_k(\theta)/2}. Marginalizing :math:`\alpha_k` out of the joint posterior yields a marginal posterior over :math:`\theta` alone, in which the uncertainty of each dataset's error calibration is propagated into the parameter estimates rather than ignored. The prior parameters :math:`a_k` and :math:`b_k` control the degree of regularisation: - :math:`a_k = b_k = 0` (Jeffreys): no prior knowledge; weights adapt freely to the data. - :math:`a_k = b_k = 2` (weakly informative): prior mean :math:`\langle\alpha_k\rangle_0 = 1`, expressing the belief that stated errors are approximately correct while allowing adaptation. - :math:`a_k = b_k = N_k/2` (strongly informative): prior mean 1 with variance :math:`\propto 1/N_k`, strongly constraining :math:`\alpha_k \approx 1`. The **EM-algorithm** applied to a Bayesian model in which each dataset carries a precision hyperparameter :math:`\alpha_k \sim \mathrm{Gamma}(a_k, b_k)` is an iterative procedure with an *Expectation* and a *Maximization* step. The general weight update at each iteration :math:`\nu` using the posterior mean is .. math:: w_k^{(\nu+1)} = \langle\alpha_k\rangle = \frac{a_k + N_k/2}{b_k + {\chi^2_k}^{(\nu)}/2} = \frac{N_k + 2a_k}{{\chi^2_k}^{(\nu)} + 2b_k}, where :math:`N_k^{\mathrm{eff}} = N_k + 2a_k` is an effective sample size inflated by the prior, and :math:`\chi^2_{\mathrm{reg},k} = 2b_k` is a regularisation floor that prevents :math:`w_k` from diverging when :math:`{\chi^2_k}^{(\nu)} \to 0`. Convergence is guaranteed because each EM iteration increases the marginal posterior :math:`P(\theta \mid \mathrm{data})`. The iterative reweighting algorithm is the EM-algorithm applied to this model: the E-step replaces :math:`\alpha_k` by its posterior mean :math:`\langle\alpha_k\rangle` given the current parameters, and the M-step maximises the resulting weighted cost over :math:`\theta`. This is equivalent to empirical Bayes — the hyperparameters :math:`\alpha_k` are estimated from the data itself rather than fixed in advance. The *regularized EM-algorithm* uses for :math:`\chi^2 minimization` the total cost function :math:`C(\theta)` (with parameters :math:`\theta = X_i,a_1,..a_p`) .. math:: C(\theta) = \frac{1}{K} \sum_{k=1}^{K} w_k \, \chi^2_{\mathrm{red},k} , \qquad \chi^2_{\mathrm{red},k} = \frac{\chi^2_k(\theta)} {N_k}, \qquad \chi^2_k(\theta) = \sum_{i \in k} \frac{\bigl(Y_i - f_i(\theta)\bigr)^2}{\sigma_i^2}, where :math:`N_k` is the number of datapoints in dataset :math:`k`. The total cost function is defined as the mean of per-dataset reduced :math:`\chi^2` values, normalised so that a good combined fit gives :math:`C \approx 1`, directly analogous to the single-dataset criterion :math:`\chi^2_{\mathrm{red}} \approx 1`: The iteration is initialised at the prior mean :math:`w_k^{(0)} = a_k/b_k`. For the Jeffreys prior the prior mean is undefined and :math:`w_k^{(0)} = 1` is used as a neutral starting point. An initial M-step with these weights produces the first residuals :math:`{\chi^2_k}^{(0)}` from which the E-step computes the first data-informed weights. At convergence, the fixed point condition .. math:: w_k^{(\infty)} \cdot {\chi^2_{\mathrm{red},k}}^{(\infty)} = 1 \qquad \forall\, k holds for every dataset individually, and the total cost satisfies :math:`\langle C \rangle = 1`. The reweighting implies that :math:`\chi^2_{red} \to 1` that it is not anymore a measure of fit quality but individual unweighted :math:`\chi^2_{\mathrm{red},k}` still describe the fit quality for each dataset and can be used to determine a better error estimate (DLS above) or to identify suspicious data with high reweighting. A short example: Zimm dynamics of a polymer in theta-solvent with DLS and NSE. *Iterative reweighting* is used to detect the wrong errors in synthetic data. .. literalinclude:: ../../examples/example_BG_fit_regularizedEM.py :language: python .. image:: ../../examples/images/iterativeReweight.jpg :align: center :height: 300px :alt: iterativeReweight - Regularisation with a non-informative prior With a `Jeffreys prior `_ (non-informative prior) using :math:`\mathrm{Gamma}(a_k \to 0, b_k \to 0)`, the posterior mean of :math:`\alpha_k` at fixed :math:`\theta` is .. math:: \langle\alpha_k\rangle = \frac{a_k + N_k/2}{b_k + \chi^2_k/2} \xrightarrow{\,a_k,\,b_k\,\to\,0\,} \frac{N_k}{\chi^2_k} = \frac{1}{\chi^2_{\mathrm{red},k}}, The weights are updated at each iteration :math:`\nu` as :math:`w_k^{(\nu+1)} = 1 / {\chi^2_{\mathrm{red},k}}^{(\nu)}`. The initial weights :math:`w_k^{(0)} = 1` correspond to equal weighting of all datasets in reduced :math:`\chi^2` units, which is the prior expectation that all datasets are correctly normalised. The Jeffreys case adds no information beyond what the data contain. So in the strict sense the Jeffreys EM iteration is not regularisation. - Regularisation with an informative prior The iterative reweighting scheme can be stabilised by replacing the Jeffreys prior on :math:`\alpha_k` with an informative :math:`\mathrm{Gamma}(a_k >0, b_k>0)` prior. This is particularly useful for small datasets where :math:`\chi^2_{\mathrm{red},k}` is noisy and the unregularised weights may fluctuate strongly between iterations. Three natural choices of prior are summarised below. Remember that for :math:`\mathrm{Gamma}(a_k, b_k)` the :math:`mean=a/b` and :math:`\sigma^2=a/b^2` . .. list-table:: :header-rows: 1 :widths: 30 15 15 40 * - Prior - :math:`a_k` - :math:`b_k` - Effect * - Jeffreys (no regularisation) - :math:`0` - :math:`0` - Weights vary freely; exact recovery of :math:`N_k/\chi^2_k`; can be noisy for small :math:`N_k` * - Weakly informative - :math:`2` - :math:`2` - :math:`\mu=1`; :math:`\sigma^2=1/2`; floor :math:`2b_k=4`; update :math:`(N_k+4)/(\chi^2_k+4)` * - Strongly informative - :math:`N_k/2` - :math:`N_k/2` - :math:`\mu=1`; :math:`\sigma^2 \propto 1/N_k`; enforces :math:`\chi^2_{\mathrm{red},k}\approx 1` as hard constraint The weakly informative choice :math:`a_k = b_k = 2` is recommended as a default. It centres :math:`\alpha_k` near 1, expressing the belief that the stated errors are approximately correct, while still allowing the weights to adapt if a dataset is systematically over- or under-weighted. The resulting update .. math:: w_k^{(\nu+1)} = \frac{N_k + 4}{{\chi^2_k}^{(\nu)} + 4} differs negligibly from the Jeffreys rule for large :math:`N_k`, but suppresses erratic weight changes when :math:`N_k` is small. The convergence criterion .. math:: w_k^{(\infty)} \cdot {\chi^2_{\mathrm{red},k}}^{(\infty)} \approx 1 \qquad \forall\, k still holds to within corrections of order :math:`2 a_k/N_k`, which vanish as :math:`N_k \to \infty`. .. note:: For the strongly informative prior :math:`a_k = b_k = N_k/2`, the weight update is .. math:: w_k^{(\nu+1)} = \frac{2N_k}{{\chi^2_k}^{(\nu)} + N_k}, and the convergence criterion evaluates to .. math:: w_k^{(\infty)} \cdot {\chi^2_{\mathrm{red},k}}^{(\infty)} = \frac{2{\chi^2_{\mathrm{red},k}}^{(\infty)}} {1 + {\chi^2_{\mathrm{red},k}}^{(\infty)}}, which equals 1 only when :math:`{\chi^2_{\mathrm{red},k}}^{(\infty)} = 1`, confirming that this prior enforces correct error calibration as a hard constraint. .. warning:: If a dataset shows :math:`\chi^2_{\mathrm{red},k} \gg 1` at convergence even with the strongly informative prior, this indicates model inadequacy rather than error misestimation. Per-dataset residuals should always be inspected after convergence. References ---------- .. [1] Robust Regression Using Iteratively Reweighted Least-Squares. Holland, P. W.; Welsch, R. E. Communications in Statistics - Theory and Methods 1977, 6 (9), 813–827. https://doi.org/10.1080/03610927708827533. .. [2] Optimal Weights and Priors in Simultaneous Fitting of Multiple Small-Angle Scattering Datasets. Larsen, A. H. J Appl Cryst 2025, 58 (3), 934–947. https://doi.org/10.1107/S1600576725002390. .. [3] Experimental Noise in Small-Angle Scattering Can Be Assessed Using the Bayesian Indirect Fourier Transformation. Larsen, A. H.; Pedersen, M. C. Journal of Applied Crystallography 2021, 54 (5), 1281–1289. https://doi.org/10.1107/S1600576721006877. .. [4] Dynamic fluctuations in a highly cross-linked polybutadiene rubber Rosi et al J. Chem. Phys. 162, 214902 (2025) https://doi.org/10.1063/5.0265589 .. [5] Fast internal dynamics in alcohol dehydrogenase. Monkenbusch, M., Stadler, A., Biehl, R., Ollivier, J., Zamponi, M. & Richter, D. Journal of Chemical Physics, 143, 075101 (2015), https://doi.org/10.1063/1.4928512