Aggregate and chain structure-factor models

[wpotrzebowski]

Submitting on behalf of @Janskovp

Summary

This pull request adds six new models under the structure-factor category, aimed at scattering from aggregates and related power-law behaviour. Most implementations use compiled C kernels with Python model definitions, following existing sasmodels conventions. A follow-up commit applies Ruff formatting to the new Python files.

New models

Model	Description
compact_polydisperse_cluster	(S(q)) for a polydisperse spherical cluster with internal correlations (Schulz polydispersity). Parameters include cluster radius, relative polydispersity, minimum scatterer separation, and weight-average aggregation number. Based on Pedersen et al. (micelles / “Babinet” particles) and related work.
fractal_aggregate	Teixeira-type fractal (S(q)) for fractal-like aggregates, reformulated so aggregation number is a fit parameter (distance between scatterers, fractal dimension, (N_\mathrm{agg})).
fractal_aggregate_discrete_chain	Same fractal core as above, with a high-(q) baseline modified using a discrete-chain / random-flight style term and a low-(q) crossover.
free_rotating_chain	Free-rotating chain (random flight) (S(q)) after Burchard & Kajiwara, with non-integer effective chain length via linear interpolation (Larsen, Pedersen & Arleth, eq. 23–24).
linear_aggregate	Linear aggregate (S(q)) with the same interpolation idea (Larsen, Pedersen & Arleth, eq. 21–22).
stabilized_power_law	Empirical (S(q) = 1 + \mathrm{amp},(0.01/q)^{\mathrm{pow}}). Listed under structure-factor; the module sets structure_factor = False—worth confirming intended behaviour in SasView vs other (S(q)) plugins.

Files

• compact_polydisperse_cluster.{py,c}
• fractal_aggregate.{py,c}
• fractal_aggregate_discrete_chain.{py,c}
• free_rotating_chain.{py,c}
• linear_aggregate.{py,c}
• stabilized_power_law.py (Python Iq only)

[wpotrzebowski]

There is some typo/nomenclature glitch. It should be compact not composite

[pkienzle]

The parameters are not following sasmodels naming conventions. Run $ python -m sasmodels.list_pars to see the parameter names used in existing models. See #305 for inconsistencies in the current models.

Existing structure factors expect radius_effective as the first parameter. This is populated from the form factor, averaged over the polydispersiy. The second parameter is volfraction, which pulls from the form factor volfraction if it is defined, with a correction for hollow shapes, wherein volfraction is the volume fraction of the shell.

Existing fractal models compute P*S at the end for P sphere or core-shell sphere. You may want to do the same for these models. The difference is in the handling of polydispersity. Keeping the structure factor separate gives you structure factor for the average sphere radius to a polydisperse form factor $S(q; \langle r \rangle) \int P(q; r) dr$, whereas the existing fractal models average over the structure factor and the sphere at that radius $\int P(q; r) S(q; r) dr$.

I implemented S(q) for a fractal so I could see the differences. Monodisperse gives identical results. Polydisperse simple P@S gives significantly different values at low Q. Those differences are much less if the β approximation (structure_factor_mode=1) is used.

$ python -m sasmodels.compare core_shell_sphere@explore/fractal_sq.py,fractal_core_shell -pars background=0 volfraction=0.05 -double -highq radius=200 cor_length=3*radius radius_pd=0.3 -lowq structure_factor_mode=0

fractal_sq.py

[Janskovp]

As a reply to @pkienzle :

The purpose of these aggregation structure factors is to make SASview more generally applicable so that the aggregation structure factors can be combined with any form factor, and hopefully in the future also with form factors of biomacromolecules calculated by AUSAXS within SASview including the scattering from the hydration layer.

The aggregation structure factors have been formulated in similar manner with aggregation number and distance between points as fit parameters. Usually, the purpose of applying then would be to determine the aggregation number, so it makes sense to have this as a fit parameter. Some structure factors have additional parameters like the fractals and compact aggregates (fractal dimension, and average radius and polydispersity, respectively). There are many systems where the distance between particles within an aggregate is not given by the radius of the constituting particles, for examples, micelles bridged by polymers/polyelectrolyte/proteins where the connecting string is often invisible in scattering and matrix mediated association, like precipitates in alloys, where lattice mismatch can lead to these effects. So to have them generally applicable, the distance between points should be a parameter.

The aggregation structure factors are different from those from liquid state theory and volume fraction is an irrelevant parameter.

There is nothing alarming in that the curves are different for different ways of including polydispersity of the constituting particles. I think the best way to preserve aggregation number as meaningful parameter, is to have the intensity as the product of the structure factor and an average form factor and not integrate the product over the size distribution. I agree that the beta/decoupling approximation is the best way to include them. I use that in my programs.

[pkienzle]

There are many systems where the distance between particles within an aggregate is not given by the radius of the constituting particles,

@Janskovp This case is covered by radius_effective_mode == 0. Then the average particle radius is ignored and the model uses the radius_effective entered by the user or fit by the optimizer. It does require that the parameter be called radius_effective, otherwise there is no way to check that the compute radial average is safely transmitted from a polydisperse form factor to the structure factor when radius_effective_mode != 0.

The aggregation structure factors are different from those from liquid state theory and volume fraction is an irrelevant parameter.

The current implementation in sasmodels/product.py requires that structure factor models have radius_effective and volfraction as the first two parameters. I put in ticket #720 to request that they be optional.

Meanwhile the structure factor models will have to include the volume fraction parameter but ignore it. This parameter is used within the structure factor to control the packing density, not as an overall scaling on P@S. That scaling happens in sasmodels/product.py after P@S is computed.

Currently there is no way to turn off the overall volume fraction scaling. For solid objects you can compensate by setting scale = 1/volfraction, but for hollow objects the combined scale is scale * volfraction / mean(shell_volume). The derived parameter mean(shell_volume) is not available to the constraints system, so there is no way compensate.

[Janskovp]

@pkienle:

I am still relatively new to SASview and I know little about what is hidden behind options that appear in the GUI. Anyway, I think that radius_effective_mode == 0. appears at 'unconstrained' in the fit panel meaning that the radius in the structure factor appears as an actual extra fit parameter, which is perfect for the versatility of the use of the cluster structure factors.

It would be good (less confusing) if the volume fraction does not appear under the structure factor when it does not depend on it. I think we could accept to keep radius_effective, although I think that the distance between points is the more natural parameter. It is easy to change the code to use radius_effective instead of the distance between points. I will talk with @wpotrzebowski Wojtek about modifying the codes to make them compatible with the requirements.

Regarding 'overall volume fraction scaling'. If I understand correctly, this is the 'scale' always appearing as the first fit parameter. For spheres, it is used to calculate number density as scale/V and there is an additional V^2 and Delta_rho^2 scaling from the form factor to give absolute intensities. I don't think that this 'interfers' with the aggregation structure factors, as it relates to the form factors. (By the way, there are similar 'issues' as with the shell structures for the polymer models (and possibly with microgel formfactors), where I_0 is used as fit parameter instead of expressing the model on absolute scale directly, mainly because the model does not directly involve the polymer volume.)

[pkienzle]

I added a note to #720 to remove the unused parameters that we just introduced.

Regarding radius_effective, I think it should be defined as half the distance between points. For example, if the radius is set by hard spheres, the distance between the centers is twice the radius. Inside the models I'm seeing code like
r = 0.5 * radius_effective; I believe this should just be r = radius_effective.

Regarding volfraction, it is applied during P@S with overall scale equal to scale/mean(volume)*volfraction. So even when it is unused inside S(q) it should still have a description like "volume fraction of material" since it applies to the P@S product.

Regarding the other parameters, here's what existing convention suggests for the names:

R_clust => radius_cluster (convention is to spell out the name)
sig_rel_R => radius_cluster_pd   (this is what I use in sasmodels scripting for the polydispersity parameter)
N_agg => n_aggreg   (used in the polymer micelle model) 
D_fract => fractal_dim (used in fractal and gelfit models)
pow => power (used in power law)
amp => peak_scale (see below)

Scale parameters don't follow the usual "math" convention of e.g., radius_core. Instead they use forms like porod_scale. Here's a list of names used in various models:

• coefficent_1: two_power_law
• i_zero: mono_gauss_coil, poly_gauss_coil
• gauss_scale: gauss_lorentz_gel
• guinier_scale: gel_fit
• lorentz_scale: correlation_length, gauss_lorentz_gel, gel_fit
• lorentz_scale_1, lorentz_scale_2: two_lorentzian
• peak_scale: broad_peak
• porod_scale: broad_peak, correlation_length

[Janskovp]

@pkienzle @wpotrzebowski

I agree with most of the renaminings. They are meaningful. But I have some comments:

"r = 0.5 * radius_effective; I believe this should just be r = radius_effective." yes, I think you are right. I will check the code but perhaps it is the other way around that in the code, one should have dist_points = 2.0 * radius_effective.

"sig_rel_R => radius_cluster_pd (this is what I use in sasmodels scripting for the polydispersity parameter)" This is a bad choice since 'polydispersity' is not well-defined (Different in SAS (sigma/average), DLS (sigma/av)^2, polymer science (1+(sigma/av)^2)). So I suggest 'radius_cluster_sigma_relative', which should define it fully.

"amp => peak_scale (see below)" Since it does not relate to a peak, I think that 'scale_power_law' is better, although it is not directly this due to the stabilization.

[pkienzle]

Yes, it should be scale_power_law except that we already have porod_scale, lorentz_scale etc. so I recommend power_law_scale instead.

Regarding the polydispersity parameter, we are using the SAS definition of 1-σ relative weights. Your implementation of SCHULZ matches that in sasmodels, so the "unmarked form" for the polydispersity parameter name radius_cluster_pd is appropriate. The other automatic polydispersity parameters are *_pd_n and *_pd_nsigma, corresponding to NPOI and hard-coded 6.0 in your P_POLY function.

The following distributions are defined in sasmodels. Here r is the nominal value of the parameter, p is the relative polydispersity and n is the number of sigmas for the distribution:

• Gaussian: exp(–½ (x-r)²/σ²) over [r – nσ, r + nσ] for σ = r p
• Uniform: 1 over [r – w, r + w] for w = r p
• Rectangle: 1 over [r – √3 w, r + √3 w] for w = r p [from σ² = (b-a)²/12 for uniform]
• LogNormal: exp(–½ (ln x – ln r)²/(σ/r)²) / (x (σ/r)) over [r – nσ, r + nσ] for σ = r p
• Schulz: $z^z R^{z-1} / [ e^{Rz} r Γ(z) ]$ over [r – nσ, r + nσ] for R = x/r, z = (r/σ)², σ = r p

Your definition of Schulz is the same after some rearrangement and change of variables. Original from SCHULZ:

$((Z+1)/R_A)^{Z+1} R^Z / [e^{(Z+1) R/R_A} Γ(Z+1)]$

with $Z$ = 1/SIGL^2 - 1, $R_Α$ = radius_cluster - radius_effective/2, and $R$ in [0, radius_effective/2 + 6 radius_cluster SIGL].

After substitution using z = Z+1, R = R/R_A and r = R_A the formulas match:

$(z/r)^z (rR)^{z-1} / [e^{z R} Γ(z)] = z^z R^{z-1} / [e^{z R} r Γ(z)]$

So the notion of pd width in sasmodels matches SIGL in your code.

The center and limits of the distribution are surprising. I would have expected R_A = radius_cluster with R over [radius_effective, radius_cluster + 6 radius_cluster SIGL].

[Janskovp]

Yes, my Schulz follows the most common definition.

The structure factor is based on the analytical solution of polydisperse spheres with a Schulz, so it is very fast. But it required some approximations to apply it in this context: R_L is the radius of the embedding sphere and R_L has relative sigma_L, but the distribution of the centers of the substructures should only be from R_L to R_L- R_S, where R_s is the same as radius_effective. There are no subparticles in spheres of radisu R_L < R_S and since R_L is the radius of the embedding sphere, we have to have R_L- R_S as the maximum value.

To use the analytical solution we neglect the problems at 0 radius, since this is anyway contributing very little, and we have to make the approximation that R_L- R_S has the same sigma as R_L (this can probably be improved for example using sigma_L * R_L/(R_L - R_S), but I did not introduce that).

There is a numerical calculation of the effective structure factor that gives the internal correlation of the subparticles. This has to be done numerically, as the spheres at the periphery (R_L- 2 * R_S) to R_L- R_S has fewer neighbors, and that there are no subparticles in a large sphere with radius R_L = R_S due to definition of R_L being the radius of the embedding sphere. Therefore the integration starts at R_S. I think/hope that for reasonable polydispersities, the upper limit of the integration is large enough. Anyway, all of this is done in order to have a reasonable estimate of local correlations with slowing the calculations down too much. And note that this is done on the 'N level' and not on the 'N^2 level' as it gives the high q behavior, so it is less justified to neglect the discrepancies at lower bound. Anyway, S_HS(q) cannot be averaged over the Schulz distribution in any case...

I am working on a document that gives these more subtle details of the model and we could include it in the documentation in the PY file later on.