Cumulants - cumulants.py

    """

DLS analysis by the method of Cumulants.

Definition
----------

THIS MODEL IS NOT INTENDED FOR THE ANALYSIS OF SAXS/SANS DATA!

This model is in part provided to illustrate the utility of SasView as a fitting
program for the analysis of input data that is *not* $I(q)$ vs $q$. But at the
same time it usefully provides a basic means of analysing Dynamic Light
Scattering (DLS) data using the Method of Cumulants (Koppel, 1972).

.. Note:: The Method of Cumulants is only valid if the particle size
          distribution is monomodal. However, the size distribution can be
          polydisperse.
      
Input Data
----------

The input data can be in any data format recognised by the SasView data loader,
though a format containing delimited two-column ASCII text will probably be the
most convenient.

SasView will look for the first pair of numerical values in the file, so it is
permissible to prefix the data with some lines of metadata and/or column
headers, but none of these will be imported. If the file contains multiple data
blocks it is likely that only the first one will be imported, so it is important
that the data you want to fit is the first data block in the file. This
reasonable degree of flexibility means that in some instances data files output
by commercial correlators can be read by SasView 'as is', without any need for
reformatting.

For example, this is a text format output by an LSi Correlator (see
https://lsinstruments.ch/en/) ::
    
    03/08/2020	17:30 PM
    Pseudo Cross Correlation
    Scattering angle:	110.0
    Duration (s):	60
    Wavelength (nm):	642.0
    Refractive index:	1.330
    Viscosity (mPas):	0.854
    Temperature (K):	309.0
    Laser intensity (mW):	0.0
    Average Count rate  A (kHz):	19388.1
    Average Count rate  B (kHz):	19388.1
    Intercept:	1.0000
    Cumulant 1st	-Inf
    Cumulant 2nd	-Inf	NaN
    Cumulant 3rd	-Inf	NaN

    Lag time (s)         g2-1
    0.000000e+00	2.153341e+02
    1.250000e-08	-1.000000e+00
    2.500000e-08	1.897621e-02
    3.750000e-08	1.454616e+00
    5.000000e-08	1.040391e+00
    6.250000e-08	9.113997e-01
    7.500000e-08	8.547216e-01
    8.750000e-08	8.269044e-01
    1.000000e-07	7.905431e-01
    1.125000e-07	7.863978e-01
    1.250000e-07	7.645010e-01
    1.375000e-07	7.728501e-01
    1.500000e-07	7.646473e-01
    1.625000e-07	7.516946e-01
    1.750000e-07	7.575662e-01
    1.875000e-07	7.556643e-01
    2.000000e-07	7.547962e-01
    ...
    5.033165e+01	-5.472073e-03

    Count Rate History (KHz)  CR CHA / CR CHB
    0.000000	18560.000000	18560.000000
    ...
    59.926118	16078.000000	16078.000000

In this example, SasView ignores all the metadata preceding the first (two-
column) data block, the correlation function, and ignores everything after the
two-column data block.

Whatever the file format it is imperative that the correct data is placed in
these the two columns::

    Column 1: what would normally be $q$:    Correlator time or $Tau$
                                             (in seconds)
    
    Column 2: what would normally be $I(q)$: Normalised Intensity Autocorrelation
                                             Function, $G2(Tau)$

Method of Analysis
------------------

For a monomodal dispersion of **monodisperse** scatterers

.. math::

    G2(	au) = 	ext{A} cdot exp(-2 Gamma 	au) + 	ext{background}

where $Gamma = D cdot q^2$ is the decay rate, $D$ is the mutual diffusion
coefficient and $q$ is the scattering vector.

For a monomodal dispersion of **polydisperse** scatterers Koppel showed

(FORMULA 1)

.. math::

    G2(	au) = 	ext{A} cdot exp(-2 Gamma_1 	au + Gamma_2 	au^2 -
    (1/3) Gamma_3 	au^3 + ...) + 	ext{background}

where each $Gamma$ is a cumulant of the distribution of decay rates arising
from the different sized scatterers. However, this function is unstable. A more
stable function is (Pusey et al, 1974; Frisken, 2001; Mailer et al, 2015)

(FORMULA 2)

.. math::

    G2(	au) = 	ext{A} cdot exp(-2 Gamma_1 	au) cdot
    (1 + (1/2) Gamma_2 	au^2 - (1/6) Gamma_3 	au^3 + ...)^2
    + 	ext{background}

Here $Gamma_1$ is now the *average* decay rate (representing a weighted average
of the different diffusion coefficients) and $Gamma_2$ is related to the
relative polydispersity index, $PDI = Gamma_2 / Gamma_1^2$.

The third cumulant is related to the *skewness* of the decay rate distribution
($Gamma_3 / Gamma_2^{3/2}$) and the fourth cumulant (not implemented in this
model) gives the *kurtosis* ($Gamma_4 / Gamma_2^2$). Further cumulants are
rarely, if ever, used.

FORMULA 2 is used by default in this model, but both formulae are provided
below. Simply comment/uncomment whichever is required.

Using this Model
----------------

The *scale* parameter will be what many commercial DLS software packages refer
to as the *Intercept* of the correlation function. This must be between 0 and 1,
but for the data to have any real meaning it should be between 0.6 and 1.

If the correlator has supplied its output as $G2(	au)-1$ (as in the data file
example above) then *background* should be essentially zero.

Do not attempt to fit all of the tail of the correlation function. Use the min
and max $q$-range for fitting boxes to limit the fit range to a sensible region
of the correlation function. You only need to fit the decay and the start of the
tail.

The best way to approach fitting the data is to gradually increase the number of
cumulants being fitted. So, initially, set *cumulant2* and *cumulant3* to zero,
and fit (check the boxes) the parameters *scale*, *background* and *cumulant1*.
Then fit *cumulant2*. Then (if required) fit *cumulant3*.

To view the fit it will probably be helpful to change the y-axis scale from the
SasView default of log(y) to y. To do this, right-click on the graph. It is
also possible to change the x and y axis labels from their SasView defaults of
$q$ and $intensity$.

Notes
-----

This model is inspired by a similar model in the SASfit model-fitting package
(Kohlbrecher, 2020).

Also see the model cumulants_dls in the SasView Marketplace.

References
----------

#.  D. E. Koppel. Analysis of macromolecular polydispersity in intensity
    correlation spectroscopy: The method of cumulants. The Journal of Chemical
    Physics, (1972), 57(11), 4814-4820. 
#.  P. N. Pusey, D. E. Koppel, D. E. Schaefer, R. D. Camerini-Otero, and S. H.
    Koenig. Intensity fluctuation spectroscopy of laser light scattered by
    solutions of spherical viruses. r17, q.beta., BSV, PM2, and t7. i. light-
    scattering technique. Biochemistry, (1974), 13(5), 952-960.
#.  B. J. Frisken. Revisiting the method of cumulants for the analysis of
    dynamic light-scattering data. Applied Optics, (2001), 40(24), 4087.
#.  A. G. Mailer, P. S. Clegg, and P. N. Pusey. Particle sizing by dynamic
    light scattering: nonlinear cumulant analysis. Journal of Physics: Condensed
    Matter, (2015), 27(14), 145102.
#.  J. Stetefeld, S. A. McKenna1 & T. R. Patel. Dynamic light scattering: a
    practical guide and applications in biomedical sciences. Biophys Rev (2016),
    8, 409-427.
#.  J. Kohlbrecher. User guide for the SASfit software package. Chapter 13.
    June 2020.

Verification
------------

The model has been tested by comparing its output to that from the commercial
LSi correlator software using the same input data and range of $	au$.

For a nominal 30 nm diameter PS latex standard::

    LSi    : intercept = 0.771, cumulant1 = 6700.0 /s,
             cumulant2 = 5.61E+06 /s^2, cumulant3 = -9.99E+08 /s^3
    SasView: intercept = 0.744, cumulant1 = 6031.1 /s,
             cumulant2 = 1.39E+05 /s^2, cumulant3 = -2.53E+10 /s^3
    
For a PNIPAM microgel::

    LSi    : intercept = 0.69, cumulant1 = 517.0 /s,
             cumulant2 = 6.45E+04 /s^2, cumulant3 = -5.0E+06 /s^3
    SasView: intercept = 0.73, cumulant1 = 511.4 /s,
             cumulant2 = 1.36E+05 /s^2, cumulant3 = -4.7E+07 /s^3
    
Visually, SasView appears to do a much better job of fitting the data than the
commercial software. This may reflect the fact that SasView has much better
optimisers. But which formula the commercial software is using is also unknown.

Authorship History
------------------

* **Author:** Steve King **Date:** 26/08/2020
* **Last Modified by:** **Date:**
* **Last Reviewed by:** **Date:**

"""
import numpy as np
from numpy import inf

name = "cumulants"
title = "Method of Cumulants"
description = """
    Computes the 1st, 2nd &
    3rd cumulants from DLS
    data using the method of
    Cumulants. READ THE DOCS!
"""
category = "shape-independent"
structure_factor = False
single = False

#   ["name", "units", default, [lower, upper], "type","description"],
parameters = [
    ["cumulant1", "/s", 6100.0, [-inf, inf], "", "1st cumulant"],
    ["cumulant2", "/s^2", 4.11e+06, [-inf, inf], "", "2nd cumulant"],
    ["cumulant3", "/s^3", 0.0, [-inf, inf], "", "3rd cumulant"],
]

def Iq(q, cumulant1, cumulant2, cumulant3):

#    FORMULA 1    
#    result = 
#    np.exp((-2.0*cumulant1*q)+(cumulant2*np.power(q,2.0))- 
#    ((1.0/3.0)*cumulant3*np.power(q,3.0)))

#   FORMULA 2
    result = 
    np.exp(-2.0*cumulant1*q)*np.power((1.0+(cumulant2*np.power(q,2.0)/2.0)- 
    (cumulant3*np.power(q,3.0)/6.0)),2.0)
    
    return result
	
Iq.vectorized = True  # Iq accepts an array of q values

tests = [
    [{'scale': 0.74552, 'background' : 0.0065923,
      'cumulant1' : 6105.5, 'cumulant2' : 4.1164e+06, 'cumulant3' : 0.0},
     [2.6e-06, 0.00983], [0.728831901773, 0.00659234803975]],
	]

Back to Model Download