In [1]:
import os
os.environ["SAS_OPENCL"] = "cuda" # use CUDA GPU backend for sasmodels
import escape as esc
import numpy as np
esc.require("0.9.8")
Loading material database from C:\dev\escape-core\python\src\escape\scattering\..\data\mdb\materials.db
SAXS. Form-factors. Core-shell parallelepiped (SasView-aligned)¶
Rectangular solid with a core-shell structure: a core of dimensions $A \times B \times C$ with three independent rim SLDs on each pair of faces. Matches core_shell_parallelepiped — SasView 6.1.3.
Reference: https://www.sasview.org/docs/user/models/core_shell_parallelepiped.html
Parameters (SasView defaults)¶
| Parameter | Variable | Value |
|---|---|---|
| Scale | scale |
1 |
| Background (cm⁻¹) | background |
0.001 |
| Core SLD (10⁻⁶ Å⁻²) | sld_core |
1 |
| A-rim SLD (10⁻⁶ Å⁻²) | sld_a |
2 |
| B-rim SLD (10⁻⁶ Å⁻²) | sld_b |
4 |
| C-rim SLD (10⁻⁶ Å⁻²) | sld_c |
2 |
| Solvent SLD (10⁻⁶ Å⁻²) | sld_solvent |
6 |
| Core side A (Å) | length_a |
35 |
| Core side B (Å) | length_b |
75 |
| Core side C (Å) | length_c |
400 |
| A-rim thickness (Å) | thick_rim_a |
10 |
| B-rim thickness (Å) | thick_rim_b |
10 |
| C-rim thickness (Å) | thick_rim_c |
10 |
| Theta (deg), 2D only | theta |
0 |
| Phi (deg), 2D only | phi |
0 |
| Psi (deg), 2D only | psi |
0 |
Form-factor¶
Using the sinc building block $S(Q_X, L) = L\,\mathrm{sinc}(Q_X L/2)$, the total amplitude is:
$$F = (\rho_{core}-\rho_{solv})\,S_A\,S_B\,S_C$$ $$\quad + (\rho_A-\rho_{solv})\,[S_{A+2t_A}-S_A]\,S_B\,S_C$$ $$\quad + (\rho_B-\rho_{solv})\,S_A\,[S_{B+2t_B}-S_B]\,S_C$$ $$\quad + (\rho_C-\rho_{solv})\,S_A\,S_B\,[S_{C+2t_C}-S_C]$$
Powder average (SasView convention — first-octant integral, normalised by $\pi/2$):
$$I(q) = \frac{\mathrm{scale}}{V} \cdot \frac{2}{\pi} \int_0^{\pi/2}\int_0^{\pi/2}F^2\sin\alpha\,d\alpha\,d\beta + \mathrm{background}$$
where $Q_A = q\sin\alpha\sin\beta$, $Q_B = q\sin\alpha\cos\beta$, $Q_C = q\cos\alpha$.
Normalisation note: SasView integrates over the first octant and normalises by $\int_0^{\pi/2}\sin\alpha\,d\alpha\cdot\int_0^{\pi/2}d\beta = 1\cdot\pi/2$. In the SasView C code this is implemented by substituting $\beta = (\pi/2)u$ and dropping the Jacobian $\pi/2$, giving an effective factor of $2/\pi$ on the double integral.
In [2]:
# ── Variables ──────────────────────────────────────────────────────────────
q = esc.var("Q")
alpha = esc.var("alpha") # polar angle
beta = esc.var("beta") # azimuthal angle
# ── Parameters ─────────────────────────────────────────────────────────────
scale = esc.par("Scale", 1.0, scale=1e8, fixed=True)
sld_core = esc.par("Core SLD", 1.0, scale=1e-6, units=f"{esc.angstr}^-2")
sld_a = esc.par("A-rim SLD", 2.0, scale=1e-6, units=f"{esc.angstr}^-2")
sld_b = esc.par("B-rim SLD", 4.0, scale=1e-6, units=f"{esc.angstr}^-2")
sld_c = esc.par("C-rim SLD", 2.0, scale=1e-6, units=f"{esc.angstr}^-2")
sld_solvent = esc.par("Solvent SLD", 6.0, scale=1e-6, units=f"{esc.angstr}^-2")
length_a = esc.par("Length A", 35.0, units=esc.angstr)
length_b = esc.par("Length B", 75.0, units=esc.angstr)
length_c = esc.par("Length C", 400.0, units=esc.angstr)
thick_rim_a = esc.par("Thick Rim A", 10.0, units=esc.angstr)
thick_rim_b = esc.par("Thick Rim B", 10.0, units=esc.angstr)
thick_rim_c = esc.par("Thick Rim C", 10.0, units=esc.angstr)
background = esc.par("Background", 0.001, userlim=[0.0, 0.03])
# ── Geometry ───────────────────────────────────────────────────────────────
# Total volume (core + 6 rim slabs, corners not double-counted)
V = (length_a * length_b * length_c
+ 2.0 * thick_rim_a * length_b * length_c
+ 2.0 * thick_rim_b * length_a * length_c
+ 2.0 * thick_rim_c * length_a * length_b)
# ── q-components ──────────────────────────────────────────────────────────
QA = q * esc.sin(alpha) * esc.sin(beta)
QB = q * esc.sin(alpha) * esc.cos(beta)
QC = q * esc.cos(alpha)
# S(Q, L) = L * sinc(Q*L/2)
SA = length_a * esc.sinc(0.5 * QA * length_a)
SB = length_b * esc.sinc(0.5 * QB * length_b)
SC = length_c * esc.sinc(0.5 * QC * length_c)
SA2 = (length_a + 2.0 * thick_rim_a) * esc.sinc(0.5 * QA * (length_a + 2.0 * thick_rim_a))
SB2 = (length_b + 2.0 * thick_rim_b) * esc.sinc(0.5 * QB * (length_b + 2.0 * thick_rim_b))
SC2 = (length_c + 2.0 * thick_rim_c) * esc.sinc(0.5 * QC * (length_c + 2.0 * thick_rim_c))
# ── Amplitude ─────────────────────────────────────────────────────────────
F = ((sld_core - sld_solvent) * SA * SB * SC
+ (sld_a - sld_solvent) * (SA2 - SA) * SB * SC
+ (sld_b - sld_solvent) * SA * (SB2 - SB) * SC
+ (sld_c - sld_solvent) * SA * SB * (SC2 - SC))
# ── Powder average ─────────────────────────────────────────────────────────
# SasView normalises the first-octant integral by pi/2 (= int_0^{pi/2} d_beta).
# Equivalently: multiply the raw double-integral by 2/pi.
inner_beta = esc.integral(
esc.pow(F, 2) * esc.sin(alpha),
beta, 0.0, np.pi / 2.0,
numpoints=61, maxiter=5, epsabs=1e-5
)
i1d = (scale / V * (2.0 / np.pi)
* esc.integral(inner_beta, alpha, 0.0, np.pi / 2.0, numpoints=61, maxiter=5, epsabs=1e-5)
+ background)
In [3]:
i1d.device = "gpu"
qs = np.linspace(0.001, 0.7, 300)
i1d.show(coordinates=qs).config(
title="Core-shell parallelepiped — powder average (1D)",
xlog=True, ylog=True,
xlabel=f"Q [{esc.angstr}^-1]", ylabel="I(q) [cm^-1]")
Out[3]:
2D oriented scattering (qx, qy)¶
Same SasView rotation convention as parallelepiped: $\theta$ in $xz$-plane, $\phi$ around $z$, $\psi$ around $C$-axis.
$$I_{\mathrm{2D}}(q_x,q_y) = \frac{\mathrm{scale}}{V}\,F^2 + \mathrm{background}$$
In [4]:
qx = esc.var("qx")
qy = esc.var("qy")
theta = esc.par("Theta", 0.0, userlim=[0.0, 180.0], units="deg")
phi = esc.par("Phi", 0.0, userlim=[-180.0, 180.0], units="deg")
psi = esc.par("Psi", 0.0, userlim=[-180.0, 180.0], units="deg")
deg = np.pi / 180.0
th = theta * deg; ph = phi * deg; ps = psi * deg
ct = esc.cos(th); st = esc.sin(th)
cp = esc.cos(ph); sp = esc.sin(ph)
ck = esc.cos(ps); sk = esc.sin(ps)
Cx, Cy = st * cp, st * sp
a0x, a0y = cp * ct, sp * ct
b0x, b0y = -sp, cp
Ax = ck * a0x + sk * b0x; Ay = ck * a0y + sk * b0y
Bx = -sk * a0x + ck * b0x; By = -sk * a0y + ck * b0y
QA2 = qx * Ax + qy * Ay
QB2 = qx * Bx + qy * By
QC2 = qx * Cx + qy * Cy
SA_2d = length_a * esc.sinc(0.5 * QA2 * length_a)
SB_2d = length_b * esc.sinc(0.5 * QB2 * length_b)
SC_2d = length_c * esc.sinc(0.5 * QC2 * length_c)
SA2_2d = (length_a + 2.0 * thick_rim_a) * esc.sinc(0.5 * QA2 * (length_a + 2.0 * thick_rim_a))
SB2_2d = (length_b + 2.0 * thick_rim_b) * esc.sinc(0.5 * QB2 * (length_b + 2.0 * thick_rim_b))
SC2_2d = (length_c + 2.0 * thick_rim_c) * esc.sinc(0.5 * QC2 * (length_c + 2.0 * thick_rim_c))
F_2d = ((sld_core - sld_solvent) * SA_2d * SB_2d * SC_2d
+ (sld_a - sld_solvent) * (SA2_2d - SA_2d) * SB_2d * SC_2d
+ (sld_b - sld_solvent) * SA_2d * (SB2_2d - SB_2d) * SC_2d
+ (sld_c - sld_solvent) * SA_2d * SB_2d * (SC2_2d - SC_2d))
i2d = scale / V * esc.pow(F_2d, 2) + background
i2d.device = "gpu"
xs = np.linspace(-0.7, 0.7, 300); ys = np.linspace(-0.7, 0.7, 300)
xv, yv = np.meshgrid(xs, ys)
coords_2d = np.column_stack([xv.flatten(), yv.flatten()]).flatten()
i2d.show(coordinates=coords_2d).config(
title="Core-shell parallelepiped — oriented 2D SAXS (qx, qy)",
xlabel=f"qx [{esc.angstr}^-1]", ylabel=f"qy [{esc.angstr}^-1]",
cblog=True, colorscale="jet")
Out[4]:
SasView reference model & comparison¶
| ESCAPE parameter | SasView parameter | Notes |
|---|---|---|
sld_core * 1e-6 |
sld_core |
core SLD (Å⁻²) |
sld_a * 1e-6 |
sld_a |
A-face rim SLD (Å⁻²) |
sld_b * 1e-6 |
sld_b |
B-face rim SLD (Å⁻²) |
sld_c * 1e-6 |
sld_c |
C-face rim SLD (Å⁻²) |
sld_solvent * 1e-6 |
sld_solvent |
solvent SLD (Å⁻²) |
length_a |
length_a |
core side A (Å) |
length_b |
length_b |
core side B (Å) |
length_c |
length_c |
core side C (Å) |
thick_rim_a |
thick_rim_a |
A-rim thickness (Å) |
thick_rim_b |
thick_rim_b |
B-rim thickness (Å) |
thick_rim_c |
thick_rim_c |
C-rim thickness (Å) |
In [5]:
import time
import matplotlib.pyplot as plt
from sasmodels.core import load_model
from sasmodels.data import empty_data1D
from sasmodels.direct_model import DirectModel
qs = np.linspace(0.001, 0.7, 300).copy()
kernel = load_model("core_shell_parallelepiped")
f_sas = DirectModel(empty_data1D(qs), kernel)
sas_pars = dict(scale=1.0, background=0.001,
sld_core=1.0, sld_a=2.0, sld_b=4.0, sld_c=2.0, sld_solvent=6.0,
length_a=35.0, length_b=75.0, length_c=400.0,
thick_rim_a=10.0, thick_rim_b=10.0, thick_rim_c=10.0)
f_sas(**sas_pars)
i1d.device = "gpu"; i1d(qs[:5])
def timeit(fn, n=5):
t0 = time.perf_counter()
for _ in range(n): result = fn()
return (time.perf_counter() - t0) / n * 1e3, result
t_sas, Iq_sas = timeit(lambda: f_sas(**sas_pars))
i1d.device = "gpu"
t_gpu, Iq_gpu = timeit(lambda: i1d(qs), n=3)
i1d.device = "cpu"
t_cpu, Iq_cpu = timeit(lambda: i1d(qs))
i1d.device = "gpu"
print(f"SASView GPU : {t_sas:.0f} ms")
print(f"ESCAPE GPU : {t_gpu:.0f} ms")
print(f"ESCAPE CPU : {t_cpu:.0f} ms ({len(qs)} q-pts)")
rel = np.max(np.abs((Iq_gpu - Iq_sas) / Iq_sas)) * 100
print(f"Max relative diff vs SasView: {rel:.2f}%")
esc.overlay(Iq_sas, Iq_gpu, Iq_cpu, coordinates=qs).config(
xlabel="Q (1/A)", ylabel="I(q) (1/cm)",
xlog=True, ylog=True,
fig_title=f"Core-shell parallelepiped I(q) — {len(qs)} pts",
labels=["SASView", "ESCAPE GPU", "ESCAPE CPU"],
line_styles=["solid", "dash", "dot"],
line_widths=[2, 3, 3]
)
C:\Users\User\AppData\Local\Temp\ipykernel_42060\2646649914.py:17: UserWarning: Input array does not own its data (e.g. it is a view or slice); data will be copied
SASView GPU : 11 ms ESCAPE GPU : 534 ms ESCAPE CPU : 407 ms (300 q-pts) Max relative diff vs SasView: 0.06%
Out[5]:
In [ ]: