Corrections to parallelepiped docs

[smk78]

This issue serves two purposes. The first is to record a useful and detailed exchange regarding the derivation and implementation of the parallelepiped model. The second is to note some potential clarifications to the parallelepiped docs that arose from that discussion.

Background
User RongchaoC emailed the developers list with the following query:

some classical SAXS books do not record the parallelepiped or there are certain printing errors in their statements. I have placed in the attachments the formulas in two books written by Professor D. I. Svergun in 1987 and 2013. The parallelepiped formula in the 1987 book is different from the parallelepiped formula in the 2013. I guess that both the formulas recorded in these two books have some printing errors, but the formula derivation is very difficult for me. I would like to know if you understand the theoretical scattering formula of parallelepiped.

@smk78 responded:

I also do not know if the formula in those Svergun references were derived by Svergun or 'copied' from another source. I suspect the latter. And I suspect that source may be the same one that SasView quotes as Reference [1] on https://www.sasview.org/docs/user/models/parallelepiped.html. But I do not have access to that reference.

However, many of the SasView models are Python translations of the models in Steve Kline's Igor-based model library for data analysis at NIST (https://www.nist.gov/ncnr/data-reduction-analysis/sans-software). Jan Ilavsky's IRENA program (https://www.aps.anl.gov/Science/Scientific-Software/Irena) also re-uses some of these Igor models too and I note that in the documentation for IRENA it says (https://saxs-igorcodedocs.readthedocs.io/en/stable/Irena/FormStructureFactors.html):

"In the case of rectangular Parallelpiped see NIST form factor description. It seems they had to go to original manuscript and recreate the form factor from the German original, Mittelbach and Porod, Acta Phys. Austriaca 14 (1961) 185-211, equations (1), (13), and (14) (in German!). Most publications citing this form factor seem to be wrong (I think there is error in Pedersen 1997 manuscript I was working with, Steven cites other manuscripts which seem to have bugs in them)."

So clearly people have been concerned about the equations for parallelepipeds in the literature for a long time!

In trying to answer your question I did a quick search to see if there was any more recent literature. I did find this paper: https://doi.org/10.1016/S0022-4073(01)00133-9 I do not know if it is helpful to you or not?

Rongchao responded:

I also noticed that the original text of the derivation of the scattering function of the parallelepiped theory is the Mittelbach & Porod 1961 (written in German) and the journal has been closed in 1985, so it is really difficult to obtain the original text. Fortunately, I have read the scattering function P(q, a, b, c) in the formula (63) on page 191 in the Pedersen's 1997 article. However, formula (63) also has a printing error. One should change qasinαcosβ, and other similar intermediate terms, into q(a/2)sinαcosβ, Because the scattering function is derived from the center (also known as the center of mass), so the upper limit of size should be half of the edge length. It is gratifying to see such a change in an article in 2017 Yimin Mao et al. J. Phys. Chem. B 2017, 121, 1340-1351. In formula (6) on page 1343 of this article, qasin#cos# is replaced by q(a/2)sin#cos#.

@pkienzle then responded thus:

The calculation in sasmodels for the parallelepiped is easiest to understand using the 3D Fourier transform.

Like all models, Iqabc() is computed for the particle in canonical orientation, with a along x, b along y and c along z. The calculator rotates (qx, qy, qz) into (qa, qb, qc) according to rotation angles (θ, φ, ψ) rather than repeating the rotation in every model.

A parallelepiped in canonical orientation is simply a product of rectangular functions in the individual dimensions. That is, ρ(x,y,z) = Π(x/a) Π(y/b) Π(z/c) where dimensions are length (a, b, c) and the rectangular function is Π(w) = 1 in [–½, +½]. Since the function is separable, the Fourier transform of the product is the product of the Fourier transform. Since the Fourier transform of Π(w) is just sinc(q) = sin(q)/(q), then the Fourier transform or ρ is the product of sinc functions. There is some bookkeeping to move the constants through the transform correctly, but that is the heart of the derivation. This is the formula given in Feigin & Svergun (1987) and used in sasmodels Iqabc().

For the 1-D pattern, it requires integrating the q vectors through (θ,φ) in [0, π/2] × [0, π/2] and relying on 8-fold symmetry for the full set of rotations. Looking at the code there is the usual u = cos θ substitution so that the outer integral is over [0,1], and a μ = q b / 2 substitution to normalize to shape width.

Note: the code and documentation for parallelepiped appear to be inconsistent, with the code using μ = qB/2 but the docs using μ = qB. However, the code is dropping a number of ½ scale factors, which appears to compensate. I think both code and docs would be better off without this substitution.

The docs should say ⟨P(q,α,β)⟩ = … rather than P(q, α) = …. The first is correct since the equations include the integration implied by ⟨ … ⟩.

I would prefer to express P in terms of θ and φ since that's what we use for 2D patterns.

⟨P(q)⟩ = ∫ ∫ F²(q,θ,φ) sin θ dφ dθ over θ=[0,π/2], φ=[0,π/2]

with

F = V Δρ sinc(a qa) sinc(b qb) sinc(c qc)
sinc(w) = sin(w)/w
(qa, qb, qc) = (q sin φ sin θ, q cos φ sin θ, q cos θ) [PAK corrected: I had written qb = q cos φ, sin θ]
(a, b, c) = (A/2, B/2, C/2)
V = volume = A B C
Δρ = contrast

After expanding F we get

⟨P(q)⟩ = V² Δρ² ∫ sin θ sinc²(c qc) ∫ sinc²(a qa) sinc²(b qb) dφ dθ over θ=[0,π/2], φ=[0,π/2]

Then the u = cos θ substitution results in sin θ = √1–u²}, and we get

⟨P(q)⟩ = V² Δρ² ∫ sinc²(c q u) ∫ sinc²(a q sin φ √1–u²} ) sinc²(b q cos φ √1–u²} ) dφ du over u=[0,1], φ=[0,π/2]

The code scales the integral by 1/4. IIRC, this is for the substitution from Gauss-Legendre integration points [-1, 1] into u = [0, 1] and φ = [0, π/2].

There may be a 4π kicking around from the number of steradians in the sphere. I don't remember why we are not scaling by 8 since we are only going over an eighth of the sphere with θ and φ.

We are scaling by V rather than V² so that we can do number density averaging over a distribution of lengths.

Actions
As noted by @pkienzle:

1. The docs should say ⟨P(q,α,β)⟩ = … rather than P(q, α) = …. The first is correct since the equations include the integration implied by ⟨ … ⟩.
2. The code and documentation for parallelepiped appear to be inconsistent, with the code using μ = qB/2 but the docs using μ = qB. However, the code is dropping a number of ½ scale factors, which appears to compensate. I think both code and docs would be better off without this substitution.