improve structure factor calculations with beta approximation

[pkienzle]

Separate calculation of F so we can return <F> and <F*F> from the I(q) calculator. These will be used for an improved structure factor calculation with non-spherical particles. The average <.> is over size dispersity and orientation.

[richardh]

The "beta(Q)" correction or adjustment to S(Q) is needed for any polydisperse particles (including spheres) and any monodipserse (or indeed polydisperse) non-spherical particle. The most useful reference is M.Kotlarchyk & S-H Chen, J.Chem.Phys 79(1983)2461-2469. Ask Richard if you would like a copy.

Their equations (14) & (15) replace the original S(Q) by S'(Q) = 1 + beta(Q)( S(Q) -1) where

beta(Q) = <F>.<F>/<F.F>  tends to "damp out" the oscillations in S(Q).

<F> will only be needed in addition to the usual <F.F> if an S(Q) is involved, so we will need a flag to show this.

Note that as usual we will have to be careful with where the scaling to get P(Q) is included, the Np.(Del rho)^2, number of particles per unit volume times scattering contrast in P(Q) will cancel in beta(Q).

Note some subtle differences in convention as to what is included in P(Q) depending on authors and what is convenient, in Kotlarchyk & Chen P(Q) is per particle.

It may also be the case for various reasons that users want to turn the beta(Q) term on/off so we will need controls to do this.

Users will also want to save beta(Q) along with at least S'(Q) and perhaps the original S(Q), so we ought somehow to make clear in sasview whether a given "S(Q)" work space includes the beta(Q) correction or not.

Note that this is all still within the "decoupling approximation" that the sizes and/or orientations of particles are totally uncorrelated with their positions.

In the real world this is no longer true as concentration increases if particles are say very polydisperse (the small ones tend to be in the gaps between the large ones) or when particles are anisotropic (they can no longer rotate freely in all directions without feeling the influcence of their neighbours). For strongly interacting particles, such as charged ones in low ionic strength solvent, the concentrations at which the theory breaks down can be relativley low.

[richardh]

Whilst we are modifying the S(Q) multiplication here for beta(Q) we might also implement some or all of #780 to allow users further controls to set or fit the effective volume and/or effective radius in S(Q). The present system that forces both of these parameters to values based on those in P(Q) is, being polite, in my opinion rather too restrictive.

[pkienzle]

Some models will be more complicated since the kernel itself is defined as using a P*S.

For example, stacked disk uses
P(q)^2 (N + 2 S(q))
where
S(q) = sum_k (N-k)cos(k d q_c) e^(k (d q_c sigma)^2)

The current software avoids this by not allowing P*S when P is stacked_disk.

[pkienzle]

Ticket #1091 (1D speed improvements) suggests unrolling the 1D integration loop to increase the exploitable parallelism on modern graphics cards. Currently parallelism is limited to the number of q points in the curve, which means a GPU with 1000 processors will be mostly idle when computing a sans curve with 100 points. While retooling the models for beta approximation, we can set them up to compute the different q_a/q_b/q_c for the given q in parallel.

We can also make progress on ticket #392 (adaptive integration). Since the q points will be evaluated independently we do not need a uniform mesh, and can instead evaluate more q_a/q_b/q_c points when evaluating larger q. If we tie the point density to the largest dimension of the shape then we can probably get away without adaptive integration.

[richardh]

I am unsure whether beta(Q) either could or should be applied to 2d data from oriented particles, perhaps others might like to comment. I suspect that we should for now restrict ourselves to 1d data.

Does anyone know of sasview users who use 2d oriented particle models with an S(Q)?

[richardh]

Replying to pkienzle:

The three lammellar stack models are similar (are they also not allowed an S(Q) ?), I suspect that since they all have an included 1d S(Q) we might need a beta(Q) if the sheet or disc thickness is polydisperse ??? However I would not make dealing with this a high priority - perhaps a separate ticket?

As discussed later (23May2018), these stacks have an internal 1d S(Q), but might also have an external inter- randomly oriented stack 3d S(Q). So we need to decide whether the internal S(Q) is one for which beta(Q) is needed, and then whether we actually can calculate <F> of the stack to use with an inter-stack beta(Q).

Some models will be more complicated since the kernel itself is defined as using a P*S.

For example, stacked disk uses
P(q)^2 (N + 2 S(q))
where
S(q) = sum_k (N-k)cos(k d q_c) e^(k (d q_c sigma)^2)

The current software avoids this by not allowing P*S when P is stacked_disk.

[richardh]

S(Q) in sasfit manula pages

[richardh]

Later we really ought to implement at least the "local monodisperse approximation" as per section 4.1.3 of the sasfit manual (to be attached), and no doubt other methods, so we should anticipate menus and controls for several possible choices.

[richardh]

Moving this ticket to beta approximation project - this is the primary ticket for the project, but needs many sub-tickets adding.

[gregsuczewski]

Above is the data points for the sphere and ellipsoid models with polydispersion using schulz distribution. I have included polydispersion without structure factor (IQD), polydispersion with simple hardsphere structure factor (IQSD), and polydispersion with beta approximation (IQBD) for both my models and Sasfit's.
I have not directly included the data from sasview because my models for IQD and IQSD are basically identical to it.

The parameters for the sphere are radius=20, sld=4,sld_solvent=1,volfraction=0.3,radius_pd=0.1, radius_pd_nsigmas=3, radius_effective=20.
The parameters for the ellipsoid are radius_polar=20, radius_equatorial=10,sld=4,sld_solvent=1,volfraction=0.3,radius_polar_pd=0.1,radius_polar_pd_nsigmas=8,radius_equatorial_pd=0, effective_radius = 13.1353356684.
There is a factor of 1e-4/volume*volfraction that I have multiplied to sasfits data. Paul and I are still figuring out why sasfit does not include this in their calculation of I(q).
For IQ and IQSD I used number density normalization (to match sasview) and for IQBD I used volume fraction normalization (to more closely match sasfits calculations). Paul, Yun and I are still discussing these differences.

[gregsuczewski]

data points for comparison between beta approximation and sasfit calculations

[richardh]

Sorry, due to other matters I have not had time to look at the test data, or generate further test data from FISH.

In terms of where we are going with the overall project (as this is the top level ticket) my personal view is that we head to the simplest possible goals to start with, all in the QT gui 5.0 version.

Default P(Q)S(Q) to:

Add S(Q) tab with controls for (a) effective radius tied to P(Q) parameters (default) or free  [posssibly later to allow further choices of ER calculation relevant to the P(Q), initial prototype in ESS_GUI_iss959 branch ]

(b) beta(Q) calc Off (default) or On [but use menus that anticipate further   choices such as local monodisperse]

FitPage Parameter tab to:

© Always show SQ_radius as a separate parameter even when it is "constrained to ER from P(Q)"

(d) Have "+" drop down on S(Q) volume fraction to reveal setttings in similar way to those for polydispersity.

Note that © has already been implemented in qt but has some issues to resolve, including getting back the value of ER to display (see #1126).

LATER we look at the sum/multi model generator and fix that to do likewise. [I believe that the "scale" factor there is currently wrong in 4.2 for P(Q)S(Q) see #1134]

[richardh]

On Jul 20, 2018, at 6:28 AM, torin.cooper-bennun@… wrote: (with replies from Paul K)

Hi Paul (et al.),

I have some things I’d like clarified regarding implementation of the beta approximation (both for my information and so I know what direction to head in, UI-wise). Am I right in the following?

1. We are only concerned (for now…?) with implementing changes in ProductModel? (and surrounding code) as opposed to MixtureModel?.

Yes.

I'm hoping that mixture is something simple like beta = (<Fa> + <Fb>)²/(<Fa²> + <Fb²>) which gives:

I = (<Fa²> + <Fb²>).(1 + (beta(S - 1))               (1)

That would lets modify mixture model to return <F> as <Fa> + <Fb> and <F²> as <Fa²> + <Fb²>, and we would be able to use the result in product. I don't know if that is a useful approximation though.

[I think it is OK, just have to accumulate <F> and <F²> as we go along. I rewrote equation (1) which was wrong in the email, Richard ]

2. The flag indicating the use of the beta approximation, or some other approximation in the future, e.g. local monodisperse, will be specific to products, therefore needs to be added as an argument to themake_product_info function (which the UI will call).

The simplest plan is to add a choice variable to the product model parameter table to select between beta/non-beta. Therefore no change in calling sequence from the rest of the program.

3. The ModelInfo? class will need to store this flag as an optional attribute (i.e. it is None when irrelevant, but has meaning when the composition is ‘product’); make_product_info will simply set the flag as necessary.

No changes necessary to ModelInfo? if beta is a parameter to the product model.

4. The implementation of the calculation itself will reside in ProductKernel?.call, where F(Q) and F2(Q) will be retrieved from the P(Q) model, etc.; it will just look for the ModelInfo? flag for indication of what calculation to perform.

Yes, except it will look at the beta flag.

Also, when it comes to choosing whether ER and VR values for S(Q) are ‘free’ or pulled from P(Q), am I right in saying the following:

· if one/both are to be free parameters, it should be specified in ModelInfo? and therefore also in the call to make_product_info, so that ProductKernel?.call will set them to the provided values instead of using calc_er_vr
· in the case of pulling them from P(Q), these ‘intermediate’ values are calculated in calc_er_vr, which can just be called directly from the UI (I assume; it takes a call_details parameter for polydispersity but I think there’s a layer of abstraction between that and the UI that I’m missing…)

The plan here is to add a new type of parameter to the model which computes a value if asked, otherwise compute nothing. This function would accept the equivalent of a Details object indicating the parameter values and dispersion points. Probably pass it in as a named tuple so the function can reference p.radius, etc. as needed, and thus simplify the function signature (should have done this for the python-based Iq function as well). The existing ER/VR would be translated into these. The availability of the Structure Factor option in the GUI would be triggered by the presence of an "F_F2" function in the model info, which can return both F and F2 to the caller, which product model will use to compute beta. If R_eff of S is tied via constraints to P.R_from_volume or P.R_from_curvature, only then will the value be substituted in, and only when the constraints are calculated.

I suggest that the GUI not show these function variables except when they are used

This may require some overhaul of the constraints processing system to get right, since now they are likely only applied when computed from the constraints tab. The UI ought to allow constraint expressions in lieu of parameter values on each model screen, and these ought to be processed in the course of evaluating that model tab alone. This is likely too ambitious for the summer.

Even more ambitious is the next level of approximation, where S*P is averaged over the dispersion (the locally-monodisperse approximation), will require a more significant restructuring. Likely we will need to generate the model on the fly, with S, P and ER all available inside OpenCL kernel. Alternatively, we may be able to pull the dispersion loops out of the kernel, and instead schedule the 10s or 100s of kernel calls that are needed either from a wrapper kernel, or perhaps from python.

A feature Richard has been asking for is to allow parameters to be tied to other parameters within the dispersion loops, for example to preserve the volume of a vesicle while allowing it to compress or expand its polar radius. Another example would be tying the cap radius to the cylinder radius in a capped cylinder. This could be done by generating kernel source on the fly with the constraints built in.