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. Parallelepiped (SasView-aligned)¶
Uniform rectangular parallelepiped with side lengths A, B, C. Matches parallelepiped — SasView 6.1.3.
Reference: https://www.sasview.org/docs/user/models/parallelepiped.html
Parameters (SasView defaults)¶
| Parameter | Variable | Value |
|---|---|---|
| Scale | scale |
1 |
| Background (cm⁻¹) | background |
0.001 |
| Contrast Δρ (10⁻⁶ Å⁻²) | contrast |
3 (= sld 4 − sld_solvent 1) |
| Side A (Å) | length_a |
35 |
| Side B (Å) | length_b |
75 |
| Side C (Å) | length_c |
400 |
| Theta (deg), 2D only | theta |
60 |
| Phi (deg), 2D only | phi |
60 |
| Psi (deg), 2D only | psi |
60 |
Form-factor¶
The amplitude is a product of stable sinc factors:
$$F(\mathbf{q}) = \Delta\rho\,V\,\mathrm{sinc}\!\left(\frac{A\,q_A}{2}\right)\mathrm{sinc}\!\left(\frac{B\,q_B}{2}\right)\mathrm{sinc}\!\left(\frac{C\,q_C}{2}\right)$$
with $V = ABC$ and $\mathrm{sinc}(x) = \sin(x)/x$.
Powder average over all orientations:
$$I(q) = \frac{\mathrm{scale}}{V}\frac{1}{4\pi}\int_0^{\pi}\int_0^{2\pi}|F|^2\sin\Theta\,d\Phi\,d\Theta + \mathrm{background}$$
In [2]:
# ── Variables ──────────────────────────────────────────────────────────────
q = esc.var("Q")
Theta = esc.var("Theta") # polar angle of q in particle frame
Phi = esc.var("Phi") # azimuthal angle of q in particle frame
# ── Parameters ─────────────────────────────────────────────────────────────
scale = esc.par("Scale", 1.0, scale=1e8, fixed=True)
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)
contrast = esc.par("Contrast", 3.0, scale=1e-6, units=f"{esc.angstr}^-2")
background = esc.par("Background", 0.001, userlim=[0.0, 0.03])
# ── Geometry ───────────────────────────────────────────────────────────────
volume = length_a * length_b * length_c
ha, hb, hc = 0.5 * length_a, 0.5 * length_b, 0.5 * length_c
# ── q-components in particle frame ────────────────────────────────────────
qA = q * esc.sin(Theta) * esc.cos(Phi)
qB = q * esc.sin(Theta) * esc.sin(Phi)
qC = q * esc.cos(Theta)
F1d = contrast * volume * esc.sinc(ha * qA) * esc.sinc(hb * qB) * esc.sinc(hc * qC)
# ── Powder average ─────────────────────────────────────────────────────────
inner_phi = esc.integral(
esc.pow(F1d, 2) * esc.sin(Theta),
Phi, 0.0, 2.0 * np.pi,
numpoints=61, maxiter=0,
)
i1d = (scale / volume * (1.0 / (4.0 * np.pi))
* esc.integral(inner_phi, Theta, 0.0, np.pi, numpoints=61, maxiter=50, epsabs=1e-5)
+ background)
In [3]:
i1d.device = "gpu"
qs = np.linspace(0.001, 0.7, 300)
i1d.show(coordinates=qs).config(
title="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)¶
For a fixed orientation $(\theta, \phi, \psi)$ the amplitude is evaluated directly at detector coordinates $(q_x, q_y)$. The rotation sequence follows SasView: first $\theta$ in the $xz$-plane, then $\phi$ around $z$, then $\psi$ around the $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", 60.0, userlim=[0.0, 180.0], units="deg")
phi = esc.par("Phi", 60.0, userlim=[-180.0, 180.0], units="deg")
psi = esc.par("Psi", 60.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)
# C-axis direction
Cx, Cy = st * cp, st * sp
# Transverse axes before psi roll
a0x, a0y = cp * ct, sp * ct
b0x, b0y = -sp, cp
# Roll around C by psi
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
F2d = contrast * volume * esc.sinc(ha * qA2) * esc.sinc(hb * qB2) * esc.sinc(hc * qC2)
i2d = scale / volume * esc.pow(F2d, 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="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 |
|---|---|---|
contrast * 1e-6 |
sld - sld_solvent |
contrast in Å⁻² |
length_a |
length_a |
side A (Å) |
length_b |
length_b |
side B (Å) |
length_c |
length_c |
side C (Å) |
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("parallelepiped")
f_sas = DirectModel(empty_data1D(qs), kernel)
sas_pars = dict(scale=1.0, background=0.001, sld=4.0, sld_solvent=1.0,
length_a=35.0, length_b=75.0, length_c=400.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"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_50632\528777119.py:15: 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 : 42 ms ESCAPE CPU : 159 ms (300 q-pts) Max relative diff vs SasView: 13.96%
Out[5]:
In [ ]: