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 multi-shell sphere (SasView-aligned)¶

Spherical core with 1 to 10 concentric shell structures, each with an independently specified SLD. Matches core_multi_shell — SasView 6.1.3.

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

Parameters (SasView defaults, 4 shells)¶

Parameter Variable Value
Scale scale 1
Background (cm⁻¹) background 0.001
Core radius (Å) radius 200
Core SLD (10⁻⁶ Å⁻²) sld_core 1
Solvent SLD (10⁻⁶ Å⁻²) sld_solvent 6.4
Shell k SLD (10⁻⁶ Å⁻²) sld[k] 1.7
Shell k thickness (Å) thickness[k] 40

Form-factor¶

The amplitude is a sum over the core and $N$ shells:

$$F(q) = \sum_{k=0}^{N} 3V_k(\rho_k - \rho_{k+1})\frac{\sin(qR_k) - qR_k\cos(qR_k)}{(qR_k)^3}$$

where $R_0 = R_{core}$, $R_k = R_{core} + \sum_{j=1}^{k} t_j$, $\rho_0 = \rho_{core}$, and $\rho_{N+1} = \rho_{solvent}$.

$$I(q) = \frac{\mathrm{scale}}{V_{outer}}F^2(q) + \mathrm{background}$$

The model is isotropic: 2D uses $q = \sqrt{q_x^2+q_y^2}$.

In [2]:
# ── Variables ──────────────────────────────────────────────────────────────
q = esc.var("Q")

# ── Parameters ─────────────────────────────────────────────────────────────
scale       = esc.par("Scale",       1.0,  scale=1e8, fixed=True)
radius      = esc.par("Core Radius",200.0, units=esc.angstr)
sld_core    = esc.par("Core SLD",    1.0,  scale=1e-6, units=f"{esc.angstr}^-2")
sld_solvent = esc.par("Solvent SLD", 6.4,  scale=1e-6, units=f"{esc.angstr}^-2")
background  = esc.par("Background",  0.001, userlim=[0.0, 0.03])

# ── Shell parameters (4 shells) ────────────────────────────────────────────
n_shells = 4
sld_shells = [
    esc.par(f"Shell {k+1} SLD",       1.7, scale=1e-6, units=f"{esc.angstr}^-2")
    for k in range(n_shells)
]
thick_shells = [
    esc.par(f"Shell {k+1} Thickness", 40.0, units=esc.angstr)
    for k in range(n_shells)
]

# ── Build radii and SLD list ───────────────────────────────────────────────
# R[0] = core radius, R[k] = R[k-1] + thickness[k-1]
radii = [radius]
for k in range(n_shells):
    radii.append(radii[-1] + thick_shells[k])

# SLD list: sld[0]=core, sld[1..n_shells]=shells, sld[n_shells+1]=solvent
slds = [sld_core] + sld_shells + [sld_solvent]

# ── Spherical kernel helper ────────────────────────────────────────────────
def sphere_kern(qr):
    return esc.conditional(
        esc.abs(qr) < 1e-10,
        1.0 / 3.0,
        (esc.sin(qr) - qr * esc.cos(qr)) / esc.pow(qr, 3),
    )

# ── Amplitude: sum over interfaces ────────────────────────────────────────
F = 0.0
for k in range(n_shells + 1):
    Rk   = radii[k]
    Vk   = 4.0 / 3.0 * np.pi * esc.pow(Rk, 3)
    drho = slds[k] - slds[k + 1]
    F    = F + 3.0 * Vk * drho * sphere_kern(q * Rk)

V_outer = 4.0 / 3.0 * np.pi * esc.pow(radii[-1], 3)
i1d     = scale / V_outer * esc.pow(F, 2) + background
In [3]:
i1d.device = "gpu"

qs = np.linspace(0.001, 0.3, 300)
i1d.show(coordinates=qs).config(
    title="Core multi-shell sphere — 1D",
    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))

F2d = 0.0
for k in range(n_shells + 1):
    Rk   = radii[k]
    Vk   = 4.0 / 3.0 * np.pi * esc.pow(Rk, 3)
    drho = slds[k] - slds[k + 1]
    F2d  = F2d + 3.0 * Vk * drho * sphere_kern(q2d * Rk)

i2d = scale / V_outer * esc.pow(F2d, 2) + background

i2d.device = "gpu"

xs = np.linspace(-0.3, 0.3, 300); ys = np.linspace(-0.3, 0.3, 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 multi-shell sphere — 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
radius radius core radius (Å)
sld_core * 1e-6 sld_core core SLD (Å⁻²)
sld_solvent * 1e-6 sld_solvent solvent SLD (Å⁻²)
sld_shells[k] * 1e-6 sld[k] shell k SLD (Å⁻²)
thick_shells[k] thickness[k] shell k thickness (Å)

Note: SasView core_multi_shell uses n as a runtime integer for the number of shells. The ESCAPE implementation above uses a fixed n_shells=4 at notebook definition time. To change the number of shells, modify n_shells and re-run the cell.

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.3, 300).copy()

kernel = load_model("core_multi_shell")
f_sas  = DirectModel(empty_data1D(qs), kernel)
sas_pars = dict(
    scale=1.0, background=0.001,
    radius=200.0, sld_core=1.0, sld_solvent=6.4,
    n=4,
    **{f"sld{k+1}": 1.7 for k in range(4)},
    **{f"thickness{k+1}": 40.0 for k in range(4)},
)

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 multi-shell sphere I(q) — {len(qs)} pts",
    labels=["SASView", "ESCAPE GPU", "ESCAPE CPU"],
    line_styles=["solid", "dash", "dot"],
    line_widths=[2, 3, 3]
)

SASView GPU : 9 ms
ESCAPE GPU  : 0 ms
ESCAPE CPU  : 14 ms  (300 q-pts)
Max relative diff vs SasView: 0.03%

C:\Users\User\AppData\Local\Temp\ipykernel_19120\3165804179.py:20: UserWarning:

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