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. Pearl necklace (SasView-aligned)¶
Colloidal spheres (pearls) connected by thin cylindrical strings, freely jointed with no preferential orientation. The model is isotropic.
Reference: https://www.sasview.org/docs/user/models/pearl_necklace.html
Parameters (SasView defaults)¶
| Parameter | Variable | Value |
|---|---|---|
| Scale | scale |
1 |
| Background (cm⁻¹) | background |
0.001 |
| Pearl radius (Å) | radius |
80 |
| Edge separation (Å) | edge_sep |
350 |
| String thickness (Å) | thick_string |
2.5 |
| Number of pearls | num_pearls |
3 |
| Pearl SLD (10⁻⁶ Å⁻²) | sld |
1 |
| String SLD (10⁻⁶ Å⁻²) | sld_string |
1 |
| Solvent SLD (10⁻⁶ Å⁻²) | sld_solvent |
6.3 |
Form-factor (Schweins 2004)¶
$$I(q) = \frac{\mathrm{scale}}{V}\cdot\frac{S_{ss}+S_{ff}+S_{fs}}{(N m_s + M m_r)^2} + \mathrm{background}$$
where $\psi(q) = 3j_1(qR)/(qR)$ (sphere form factor), $\Lambda(q) = \mathrm{Si}(ql)/(ql)$ (string factor), $A = 2R + l$ (centre-to-centre distance), $N$ = number of pearls, $M = N-1$ strings.
The 2D scattering pattern is the same as the 1D curve evaluated at $q = \sqrt{q_x^2+q_y^2}$ (isotropic model).
In [2]:
# ── Variables ──────────────────────────────────────────────────────────────
q = esc.var("Q")
t_si = esc.var("t_si") # integration variable for sine integral Si(x)
# ── Parameters ─────────────────────────────────────────────────────────────
scale = esc.par("Scale", 1.0, scale=1e8, fixed=True)
radius = esc.par("Radius", 80.0, units=esc.angstr)
edge_sep = esc.par("Edge sep", 350.0, units=esc.angstr)
thick_string = esc.par("Thick string", 2.5, units=esc.angstr)
num_pearls = esc.par("Num pearls", 3.0, userlim=[1.0, 20.0], fixed=True)
sld = esc.par("SLD pearl", 1.0, scale=1e-6, units=f"{esc.angstr}^-2")
sld_string = esc.par("SLD string", 1.0, scale=1e-6, units=f"{esc.angstr}^-2")
sld_solvent = esc.par("SLD solvent", 6.3, scale=1e-6, units=f"{esc.angstr}^-2")
background = esc.par("Background", 0.001, userlim=[0.0, 0.03])
# ── Derived quantities ─────────────────────────────────────────────────────
N = num_pearls
M = N - 1.0 # number of strings
l = edge_sep # string length (edge-to-edge)
A = 2.0 * radius + l # pearl centre-to-centre distance
v_pearl = (4.0 / 3.0) * np.pi * esc.pow(radius, 3)
v_string = np.pi * esc.pow(thick_string / 2.0, 2) * l
v_total = N * v_pearl + M * v_string
# Contrast-weighted masses
ms = (sld - sld_solvent) * v_pearl # pearl mass
mr = (sld_string - sld_solvent) * v_string # string mass
# ── Sphere form factor: psi = 3*j1(qR)/(qR) ────────────────────────────────
# j1(x)/x = (sin(x) - x*cos(x)) / x^3 → 3*j1(qR)/(qR) = 3*(sin-x*cos)/x^3
x_s = q * radius
psi = 3.0 * (esc.sin(x_s) - x_s * esc.cos(x_s)) / esc.pow(x_s, 3)
# ── String form factor: Lambda = Si(ql)/(ql) ───────────────────────────────
# Si(x) = integral_0^x sin(t)/t dt = integral_0^x sinc(t) dt
x_l = q * l
si_ql = esc.integral(esc.sinc(t_si), t_si, 0.0, x_l, numpoints=61, maxiter=5, epsabs=1e-5)
Lambda = si_ql / x_l
# ── beta = (Si(q*(A-R)) - Si(q*R)) / (q*l) ────────────────────────────────
x_AR = q * (A - radius)
x_R = q * radius
si_AR = esc.integral(esc.sinc(t_si), t_si, 0.0, x_AR, numpoints=61, maxiter=5, epsabs=1e-5)
si_R = esc.integral(esc.sinc(t_si), t_si, 0.0, x_R, numpoints=61, maxiter=5, epsabs=1e-5)
beta = (si_AR - si_R) / x_l
# sinc(qA) = sin(qA)/(qA)
sincA = esc.sinc(q * A)
# ── Schweins structure factors ──────────────────────────────────────────────
# Sss: pearl-pearl correlations
Sss = (2.0 * esc.pow(ms, 2) * esc.pow(psi, 2)
* (N / (1.0 - sincA) - N / 2.0
- (1.0 - esc.pow(sincA, N)) / esc.pow(1.0 - sincA, 2) * sincA))
# Sff: string-string correlations
Sff = (esc.pow(mr, 2)
* (M * (2.0 * Lambda - esc.sinc(0.5 * q * l))
+ 2.0 * M * esc.pow(beta, 2) / (1.0 - sincA)
- 2.0 * esc.pow(beta, 2) * (1.0 - esc.pow(sincA, M)) / esc.pow(1.0 - sincA, 2)))
# Sfs: pearl-string cross-correlations
Sfs = (mr * beta * ms * psi
* 4.0 * ((N - 1.0) / (1.0 - sincA)
- (1.0 - esc.pow(sincA, N - 1.0)) / esc.pow(1.0 - sincA, 2) * sincA))
i1d = scale / v_total * (Sss + Sff + Sfs) + background
In [3]:
i1d.device = "gpu"
qs = np.linspace(0.001, 0.2, 300)
i1d.show(coordinates=qs).config(
title="Pearl necklace — isotropic 1D",
xlog=True, ylog=True,
xlabel=f"Q [{esc.angstr}^-1]", ylabel="I(q) [cm^-1]")
Out[3]:
2D scattering (isotropic)¶
The pearl necklace model is isotropic: the 2D scattering pattern is identical to the 1D curve evaluated at $q = \sqrt{q_x^2+q_y^2}$. No orientation parameters are defined.
In [4]:
qx = esc.var("qx")
qy = esc.var("qy")
t_si2 = esc.var("t_si2") # separate integration variable for 2D
q_mag = esc.sqrt(esc.pow(qx, 2) + esc.pow(qy, 2))
x_s_2d = q_mag * radius
psi_2d = 3.0 * (esc.sin(x_s_2d) - x_s_2d * esc.cos(x_s_2d)) / esc.pow(x_s_2d, 3)
x_l_2d = q_mag * l
si_ql_2d = esc.integral(esc.sinc(t_si2), t_si2, 0.0, x_l_2d, numpoints=61, maxiter=0)
Lambda_2d = si_ql_2d / x_l_2d
x_AR_2d = q_mag * (A - radius)
x_R_2d = q_mag * radius
si_AR_2d = esc.integral(esc.sinc(t_si2), t_si2, 0.0, x_AR_2d, numpoints=61, maxiter=0)
si_R_2d = esc.integral(esc.sinc(t_si2), t_si2, 0.0, x_R_2d, numpoints=61, maxiter=0)
beta_2d = (si_AR_2d - si_R_2d) / x_l_2d
sincA_2d = esc.sinc(q_mag * A)
Sss_2d = (2.0 * esc.pow(ms, 2) * esc.pow(psi_2d, 2)
* (N / (1.0 - sincA_2d) - N / 2.0
- (1.0 - esc.pow(sincA_2d, N)) / esc.pow(1.0 - sincA_2d, 2) * sincA_2d))
Sff_2d = (esc.pow(mr, 2)
* (M * (2.0 * Lambda_2d - esc.sinc(0.5 * q_mag * l))
+ 2.0 * M * esc.pow(beta_2d, 2) / (1.0 - sincA_2d)
- 2.0 * esc.pow(beta_2d, 2) * (1.0 - esc.pow(sincA_2d, M)) / esc.pow(1.0 - sincA_2d, 2)))
Sfs_2d = (mr * beta_2d * ms * psi_2d
* 4.0 * ((N - 1.0) / (1.0 - sincA_2d)
- (1.0 - esc.pow(sincA_2d, N - 1.0)) / esc.pow(1.0 - sincA_2d, 2) * sincA_2d))
i2d = scale / v_total * (Sss_2d + Sff_2d + Sfs_2d) + background
i2d.device = "gpu"
xs = np.linspace(-0.2, 0.2, 300); ys = np.linspace(-0.2, 0.2, 300)
xv, yv = np.meshgrid(xs, ys)
coords_2d = np.column_stack([xv.flatten(), yv.flatten()]).flatten()
i2d.show(coordinates=coords_2d).config(
title="Pearl necklace — isotropic 2D (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 * 1e-6 |
sld |
pearl SLD (Å⁻²) |
sld_string * 1e-6 |
sld_string |
string SLD (Å⁻²) |
sld_solvent * 1e-6 |
sld_solvent |
solvent SLD (Å⁻²) |
radius |
radius |
pearl radius (Å) |
edge_sep |
edge_sep |
edge-to-edge separation (Å) |
thick_string |
thick_string |
string diameter (Å) |
num_pearls |
num_pearls |
number of pearls |
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.2, 300).copy()
kernel = load_model("pearl_necklace")
f_sas = DirectModel(empty_data1D(qs), kernel)
sas_pars = dict(scale=1.0, background=0.001,
sld=1.0, sld_string=1.0, sld_solvent=6.3,
radius=80.0, edge_sep=350.0, thick_string=2.5, num_pearls=3)
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"Pearl Necklace I(q) — {len(qs)} pts",
labels=["SASView", "ESCAPE GPU", "ESCAPE CPU"],
line_styles=["solid", "dash", "dot"],
line_widths=[2, 3, 3]
)
SASView GPU : 10 ms ESCAPE GPU : 0 ms ESCAPE CPU : 19 ms (300 q-pts) Max relative diff vs SasView: 5.19%
C:\Users\User\AppData\Local\Temp\ipykernel_56232\1150710002.py:16: UserWarning: Input array does not own its data (e.g. it is a view or slice); data will be copied
Out[5]:
In [ ]: