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 with head groups (SasView-aligned)

Random lamellar phase with a head-tail-tail-head bilayer structure: separate SLDs for tail and head regions. Matches lamellar_hg — SasView 6.1.3.

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

Parameters (SasView defaults)

Parameter	Variable	Value
Scale	scale	1
Background (cm⁻¹)	background	0.001
Tail thickness δT (Å)	length_tail	15
Head thickness δH (Å)	length_head	10
Tail SLD (10⁻⁶ Å⁻²)	sld_tail	0.4
Head SLD (10⁻⁶ Å⁻²)	sld_head	3
Solvent SLD (10⁻⁶ Å⁻²)	sld_solvent	6

Form-factor

Total bilayer thickness: $\delta = 2(\delta_H + \delta_T)$.

$$P(q) = \frac{4}{q^2}\bigl\{\Delta\rho_H[\sin(q(\delta_H+\delta_T)) - \sin(q\delta_T)] + \Delta\rho_T\sin(q\delta_T)\bigr\}^2$$

where $\Delta\rho_H = \rho_{head} - \rho_{solvent}$ and $\Delta\rho_T = \rho_{tail} - \rho_{solvent}$.

$$I(q) = \frac{\mathrm{scale}\cdot 2\pi\cdot P(q)}{2(\delta_H+\delta_T)\cdot q^2} + \mathrm{background}$$

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

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

# ── Parameters ─────────────────────────────────────────────────────────────
scale        = esc.par("Scale",        1.0,  scale=1e8, fixed=True)
length_tail  = esc.par("Tail Length", 15.0,  units=esc.angstr)
length_head  = esc.par("Head Length", 10.0,  units=esc.angstr)
sld_tail     = esc.par("Tail SLD",    0.4,   scale=1e-6, units=f"{esc.angstr}^-2")
sld_head     = esc.par("Head SLD",    3.0,   scale=1e-6, units=f"{esc.angstr}^-2")
sld_solvent  = esc.par("Solvent SLD", 6.0,   scale=1e-6, units=f"{esc.angstr}^-2")
background   = esc.par("Background",  0.001, userlim=[0.0, 0.03])

# ── Contrasts ─────────────────────────────────────────────────────────────
drho_H = sld_head - sld_solvent
drho_T = sld_tail - sld_solvent

delta = 2.0 * (length_head + length_tail)   # total bilayer thickness

# ── Form factor ────────────────────────────────────────────────────────────
amp = (drho_H * (esc.sin(q * (length_head + length_tail)) - esc.sin(q * length_tail))
       + drho_T * esc.sin(q * length_tail))

P_q = 4.0 / esc.pow(q, 2) * esc.pow(amp, 2)

i1d = scale * 2.0 * np.pi * P_q / (delta * esc.pow(q, 2)) + background

i1d.device = "gpu"

qs = np.linspace(0.01, 1.0, 300)
i1d.show(coordinates=qs).config(
    title="Lamellar HG — 1D",
    xlog=True, ylog=True,
    xlabel=f"Q [{esc.angstr}^-1]", ylabel="I(q) [cm^-1]")

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}$.

qx = esc.var("qx")
qy = esc.var("qy")
q2d = esc.sqrt(esc.pow(qx, 2) + esc.pow(qy, 2))

amp_2d = (drho_H * (esc.sin(q2d * (length_head + length_tail)) - esc.sin(q2d * length_tail))
          + drho_T * esc.sin(q2d * length_tail))

P_2d = 4.0 / esc.pow(q2d, 2) * esc.pow(amp_2d, 2)
i2d  = scale * 2.0 * np.pi * P_2d / (delta * esc.pow(q2d, 2)) + 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="Lamellar HG — isotropic 2D (qx, qy)",
    xlabel=f"qx [{esc.angstr}^-1]", ylabel=f"qy [{esc.angstr}^-1]",
    cblog=True, colorscale="jet")

SasView reference model & comparison

ESCAPE parameter	SasView parameter	Notes
length_tail	length_tail	tail thickness δT (Å)
length_head	length_head	head thickness δH (Å)
sld_tail * 1e-6	sld	tail SLD (Å⁻²)
sld_head * 1e-6	sld_head	head SLD (Å⁻²)
sld_solvent * 1e-6	sld_solvent	solvent SLD (Å⁻²)

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

kernel = load_model("lamellar_hg")
f_sas  = DirectModel(empty_data1D(qs), kernel)
sas_pars = dict(scale=1.0, background=0.001,
                length_tail=15.0, length_head=10.0,
                sld=0.4, sld_head=3.0, sld_solvent=6.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"Lamellar HG I(q) — {len(qs)} pts",
    labels=["SASView", "ESCAPE GPU", "ESCAPE CPU"],
    line_styles=["solid", "dash", "dot"],
    line_widths=[2, 3, 3]
)

SASView GPU : 7 ms
ESCAPE GPU  : 0 ms
ESCAPE CPU  : 2 ms  (300 q-pts)
Max relative diff vs SasView: 0.00%

C:\Users\User\AppData\Local\Temp\ipykernel_60080\2514463710.py:16: UserWarning:

Input array does not own its data (e.g. it is a view or slice); data will be copied