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. Lamellar stack with paracrystal structure factor (SasView-aligned)¶

Stacked lamellar sheets with a paracrystal structure factor: Gaussian polydispersity in the repeat spacing. Matches lamellar_stack_paracrystal — SasView 6.1.3.

Reference: https://www.sasview.org/docs/user/models/lamellar_stack_paracrystal.html

Parameters (SasView defaults)¶

Parameter Variable Value
Scale (Γm, mass per area) scale 1
Background (cm⁻¹) background 0.001
Layer thickness t (Å) thickness 33
Number of layers N N 20 (integer, fixed at functor build time)
Mean d-spacing ⟨D⟩ (Å) d_spacing 250
Relative σ of d-spacing sigma_d 0.2
Contrast Δρ (10⁻⁶ Å⁻²) contrast 5.34 (= sld 1 − sld_solvent 6.34)

Form-factor × Structure factor (Bergstrom et al., 1999)¶

$$P_{bil}(q) = \left(\frac{\sin(qt/2)}{qt/2}\right)^2 = \mathrm{sinc}^2(qt/2)$$

$$w = \exp\!\left(-\frac{\sigma_D^2 q^2}{2}\right), \quad \sigma_D = \sigma_d\cdot\langle D\rangle$$

$$Z_N(q) = \frac{1-w^2}{1+w^2-2w\cos(q\langle D\rangle)} + \frac{a_N}{N[1+w^2-2w\cos(q\langle D\rangle)]^2}$$

$$a_N = 4w^2 - 2(w^3+w)\cos(q\langle D\rangle) - 4w^{N+2}\cos(Nq\langle D\rangle) + 2w^{N+3}\cos[(N-1)q\langle D\rangle] + 2w^{N+1}\cos[(N+1)q\langle D\rangle]$$

$$I(q) = 2\pi(\Delta\rho)^2\Gamma_m\frac{P_{bil}(q)}{q^2}Z_N(q) + \mathrm{background}$$

Note on N: The Bergstrom formula for $a_N$ is derived for integer $N$ via a telescoping sum. N is therefore a fixed Python int baked into the ESCAPE functor at build time.

Note on sigma_d: The denominator $1+w^2-2w\cos(qD)$ is zero when $w=1$ (i.e. sigma_d=0) at every Bragg peak $qD=2\pi k$. In that limit $Z_N \to N$. For sigma_d>0 the denominator is always positive and the formula is well-behaved. The default sigma_d=0.2 matches the SasView reference test.

In [2]:
# ── Fixed integer N (baked into the functor at build time) ─────────────────
N = 20   # integer — change and re-run this cell to rebuild the functor

# ── Variables ──────────────────────────────────────────────────────────────
q = esc.var("Q")

# ── Parameters ─────────────────────────────────────────────────────────────
scale      = esc.par("Scale",      1.0,  scale=1e8, fixed=True)
thickness  = esc.par("Thickness", 33.0,  units=esc.angstr)
d_spacing  = esc.par("d-spacing",250.0,  units=esc.angstr)
sigma_d    = esc.par("sigma_d",    0.2,  userlim=[0.0, 1.0])   # relative; 0.2 = SasView test default
contrast   = esc.par("Contrast",   5.34, scale=1e-6, units=f"{esc.angstr}^-2")
background = esc.par("Background", 0.001, userlim=[0.0, 0.03])

# ── Bilayer form factor: Pbil = sinc^2(qt/2) ──────────────────────────────
Pbil = esc.pow(esc.sinc(0.5 * q * thickness), 2)

# ── Paracrystal structure factor Z_N(q) ───────────────────────────────────
sigma_D = sigma_d * d_spacing                          # absolute σ of spacing
w       = esc.exp(-0.5 * esc.pow(sigma_D, 2) * esc.pow(q, 2))
cos_qD  = esc.cos(q * d_spacing)
denom   = 1.0 + esc.pow(w, 2) - 2.0 * w * cos_qD

# a_N: closed-form for integer N (Bergstrom eq. 5)
a_N = (4.0 * esc.pow(w, 2)
       - 2.0 * (esc.pow(w, 3) + w) * cos_qD
       - 4.0 * esc.pow(w, N + 2) * esc.cos(N       * q * d_spacing)
       + 2.0 * esc.pow(w, N + 3) * esc.cos((N - 1) * q * d_spacing)
       + 2.0 * esc.pow(w, N + 1) * esc.cos((N + 1) * q * d_spacing))

# Guard denom against zero.
# When sigma_d=0 and qD=2*pi*k: w=1, denom=0, and Z_N -> N analytically.
# Threshold 1e-6 is safe: with sigma_d>=0.01 the minimum denom is ~1e-4.
Z_N = esc.conditional(
    denom < 1e-6,
    float(N),
    (1.0 - esc.pow(w, 2)) / denom + a_N / (N * esc.pow(denom, 2))
)

# ── Combined intensity ─────────────────────────────────────────────────────
i1d = (2.0 * np.pi * esc.pow(contrast, 2) * scale
       * Pbil / esc.pow(q, 2) * Z_N
       + background)
In [3]:
i1d.device = "gpu"

qs = np.linspace(0.001, 0.5, 300)
i1d.show(coordinates=qs).config(
    title=f"Lamellar stack paracrystal — powder average (1D), N={N}",
    xlog=True, ylog=True,
    xlabel=f"Q [{esc.angstr}^-1]", ylabel="I(q) [cm^-1]")
Out[3]:

2D isotropic scattering (qx, qy)¶

The model is isotropic. The 2D intensity is the same as 1D with $q = \sqrt{q_x^2+q_y^2}$.

In [4]:
qx = esc.var("qx")
qy = esc.var("qy")
q2d = esc.sqrt(esc.pow(qx, 2) + esc.pow(qy, 2))

Pbil_2d   = esc.pow(esc.sinc(0.5 * q2d * thickness), 2)
w_2d      = esc.exp(-0.5 * esc.pow(sigma_D, 2) * esc.pow(q2d, 2))
cos_qD_2d = esc.cos(q2d * d_spacing)
denom_2d  = 1.0 + esc.pow(w_2d, 2) - 2.0 * w_2d * cos_qD_2d

a_N_2d = (4.0 * esc.pow(w_2d, 2)
          - 2.0 * (esc.pow(w_2d, 3) + w_2d) * cos_qD_2d
          - 4.0 * esc.pow(w_2d, N + 2) * esc.cos(N       * q2d * d_spacing)
          + 2.0 * esc.pow(w_2d, N + 3) * esc.cos((N - 1) * q2d * d_spacing)
          + 2.0 * esc.pow(w_2d, N + 1) * esc.cos((N + 1) * q2d * d_spacing))

Z_N_2d = esc.conditional(
    denom_2d < 1e-6,
    float(N),
    (1.0 - esc.pow(w_2d, 2)) / denom_2d + a_N_2d / (N * esc.pow(denom_2d, 2))
)

i2d = (2.0 * np.pi * esc.pow(contrast, 2) * scale
       * Pbil_2d / esc.pow(q2d, 2) * Z_N_2d
       + background)

i2d.device = "gpu"

xs = np.linspace(-0.5, 0.5, 300); ys = np.linspace(-0.5, 0.5, 300)
xv, yv = np.meshgrid(xs, ys)
coords_2d = np.column_stack([xv.flatten(), yv.flatten()]).flatten()
i2d.show(coordinates=coords_2d).config(
    title=f"Lamellar stack paracrystal — isotropic 2D, N={N}",
    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 Å⁻²
thickness thickness layer thickness t (Å)
N (Python int) Nlayers number of layers — fixed at build time
d_spacing d_spacing mean d-spacing ⟨D⟩ (Å)
sigma_d sigma_d relative σ of d-spacing (use >0 for physical results)
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.5, 300).copy()

kernel = load_model("lamellar_stack_paracrystal")
f_sas  = DirectModel(empty_data1D(qs), kernel)
sas_pars = dict(scale=1.0, background=0.001, sld=1.0, sld_solvent=6.34,
                thickness=33.0, Nlayers=N, d_spacing=250.0, sigma_d=0.2)

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"Lamellar stack paracrystal I(q) — N={N}, {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  : 0 ms
ESCAPE CPU  : 16 ms  (300 q-pts)
Max relative diff vs SasView: 0.00%

C:\Users\User\AppData\Local\Temp\ipykernel_34320\752763899.py:15: UserWarning:

Input array does not own its data (e.g. it is a view or slice); data will be copied
Out[5]:
In [ ]: