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. Stacked disks (SasView-aligned)¶
Stacked core/shell discs (tactoids) with a Kratky-Porod structure factor for the inter-disc spacing.
Reference: https://www.sasview.org/docs/user/models/stacked_disks.html
Parameters (SasView defaults)¶
| Parameter | Variable | Value |
|---|---|---|
| Scale | scale |
1 |
| Background (cm⁻¹) | background |
0.001 |
| Core half-thickness (Å) | thick_core |
10 |
| Layer thickness (Å) | thick_layer |
10 |
| Radius (Å) | radius |
15 |
| Number of discs | n_stacking |
1 |
| σ_d (Å) | sigma_d |
0 |
| Core SLD (10⁻⁶ Å⁻²) | sld_core |
4 |
| Layer SLD (10⁻⁶ Å⁻²) | sld_layer |
0 |
| Solvent SLD (10⁻⁶ Å⁻²) | sld_solvent |
5 |
| Theta (deg), 2D only | theta |
0 |
| Phi (deg), 2D only | phi |
0 |
Form-factor (SasView stacked_disks.c)¶
Each single disc (layer/core/layer) amplitude:
$$t_1 = \pi R^2\cdot 2h\cdot(\rho_c-\rho_s)\,\mathrm{sinc}(h\,q_\parallel)\,\frac{2J_1(R\,q_\perp)}{R\,q_\perp}$$
$$t_2 = \pi R^2\cdot(\rho_l-\rho_s)\left[D\,\mathrm{sinc}((h+d)q_\parallel)-2h\,\mathrm{sinc}(h\,q_\parallel)\right]\frac{2J_1(R\,q_\perp)}{R\,q_\perp}$$
Structure factor for $n$ stacked discs (Kratky-Porod):
$$S(q,\alpha) = 1 + \frac{2}{n}\sum_{k=1}^{n-1}(n-k)\cos(kD\,q_\parallel)\exp\!\left(-\tfrac{k}{2}(D\,q_\parallel\,\sigma_d)^2\right)$$
In [2]:
# ── Variables ──────────────────────────────────────────────────────────────
q = esc.var("Q")
alpha = esc.var("alpha") # angle between disc normal and q
# ── Parameters ─────────────────────────────────────────────────────────────
scale = esc.par("Scale", 1.0, scale=1e8, fixed=True)
thick_core = esc.par("Thick core", 10.0, units=esc.angstr)
thick_layer = esc.par("Thick layer", 10.0, units=esc.angstr)
radius = esc.par("Radius", 15.0, units=esc.angstr)
n_stacking = esc.par("N stacking", 1.0, userlim=[1.0, 20.0], fixed=True)
sigma_d = esc.par("Sigma d", 0.0, units=esc.angstr)
sld_core = esc.par("SLD core", 4.0, scale=1e-6, units=f"{esc.angstr}^-2")
sld_layer = esc.par("SLD layer", 0.0, scale=1e-6, units=f"{esc.angstr}^-2")
sld_solvent = esc.par("SLD solvent", 5.0, scale=1e-6, units=f"{esc.angstr}^-2")
background = esc.par("Background", 0.001, userlim=[0.0, 0.03])
# ── Geometry ───────────────────────────────────────────────────────────────
h = 0.5 * thick_core # core half-thickness
d = thick_layer # layer thickness
D = 2.0 * (d + h) # centre-to-centre disc spacing
area = np.pi * esc.pow(radius, 2)
v_single = area * D # volume of one disc unit (layer/core/layer)
dr_core = sld_core - sld_solvent # core contrast
dr_layer = sld_layer - sld_solvent # layer contrast
# ── Oriented amplitude ─────────────────────────────────────────────────────
# q_axial = q * cos(alpha): component along disc normal
# q_radial = q * sin(alpha): component in disc plane
q_axial = q * esc.cos(alpha)
q_radial = q * esc.sin(alpha)
be = 2.0 * esc.j1_over_x(radius * q_radial) # Bessel factor
si1 = esc.sinc(h * q_axial) # sinc for core half-thickness
si2 = esc.sinc((h + d) * q_axial) # sinc for core + layer
t1 = area * 2.0 * h * dr_core * si1 * be
t2 = area * dr_layer * (D * si2 - 2.0 * h * si1) * be
pq = esc.pow(t1 + t2, 2)
# ── Structure factor (Python-level sum, n_stacking is fixed) ───────────────
# For n=1 the sum is empty and S(q)=1
n_int = 1 # default n_stacking value; update this when changing n_stacking
sq_sum = 0.0
for k in range(1, n_int):
qd_cos = D * q_axial
sq_sum = sq_sum + (n_int - k) * esc.cos(k * qd_cos) * esc.exp(-0.5 * k * esc.pow(qd_cos * sigma_d, 2))
sq = 1.0 + 2.0 * sq_sum / n_stacking
# ── Powder average ─────────────────────────────────────────────────────────
i1d = (scale / v_single
* esc.integral(pq * sq * esc.sin(alpha),
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, 1.0, 300)
i1d.show(coordinates=qs).config(
title="Stacked disks — 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)¶
For a fixed orientation $(\theta, \phi)$ the amplitude is evaluated directly at detector coordinates. For $n=1$ the structure factor is $S(q)=1$.
$$I_{\mathrm{2D}}(q_x,q_y) = \frac{\mathrm{scale}}{V_{\mathrm{single}}}\,(t_1+t_2)^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=[0.0, 360.0], units="deg")
deg = np.pi / 180.0
sin_t = esc.sin(theta * deg)
ux = sin_t * esc.cos(phi * deg)
uy = sin_t * esc.sin(phi * deg)
q_sq = esc.pow(qx, 2) + esc.pow(qy, 2)
q_par_2d = qx * ux + qy * uy
q_perp_2d = esc.sqrt(q_sq - esc.pow(q_par_2d, 2))
be_2d = 2.0 * esc.j1_over_x(radius * q_perp_2d)
si1_2d = esc.sinc(h * q_par_2d)
si2_2d = esc.sinc((h + d) * q_par_2d)
t1_2d = area * 2.0 * h * dr_core * si1_2d * be_2d
t2_2d = area * dr_layer * (D * si2_2d - 2.0 * h * si1_2d) * be_2d
pq_2d = esc.pow(t1_2d + t2_2d, 2)
# S(q) = 1 for n=1
i2d = scale / v_single * pq_2d + background
i2d.device = "gpu"
xs = np.linspace(-1.0, 1.0, 300); ys = np.linspace(-1.0, 1.0, 300)
xv, yv = np.meshgrid(xs, ys)
coords_2d = np.column_stack([xv.flatten(), yv.flatten()]).flatten()
i2d.show(coordinates=coords_2d).config(
title="Stacked disks — 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_layer * 1e-6 |
sld_layer |
layer SLD (Å⁻²) |
sld_solvent * 1e-6 |
sld_solvent |
solvent SLD (Å⁻²) |
thick_core |
thick_core |
core half-thickness (Å) |
thick_layer |
thick_layer |
layer thickness (Å) |
radius |
radius |
disc radius (Å) |
n_stacking |
n_stacking |
number of stacked discs |
sigma_d |
sigma_d |
spacing disorder (Å) |
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, 1.0, 300).copy()
kernel = load_model("stacked_disks")
f_sas = DirectModel(empty_data1D(qs), kernel)
sas_pars = dict(scale=1.0, background=0.001,
sld_core=4.0, sld_layer=0.0, sld_solvent=5.0,
thick_core=10.0, thick_layer=10.0, radius=15.0,
n_stacking=1, sigma_d=0.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"Stacked discs I(q) — {len(qs)} pts",
labels=["SASView", "ESCAPE GPU", "ESCAPE CPU"],
line_styles=["solid", "dash", "dot"],
line_widths=[2, 3, 3]
)
SASView GPU : 11 ms ESCAPE GPU : 1 ms ESCAPE CPU : 8 ms (300 q-pts) Max relative diff vs SasView: 0.00%
C:\Users\User\AppData\Local\Temp\ipykernel_47184\4133630627.py:17: UserWarning: Input array does not own its data (e.g. it is a view or slice); data will be copied
Out[5]:
In [ ]:
In [ ]: