diff --git a/molgroups/refl1d_interface/__init__.py b/molgroups/refl1d_interface/__init__.py index 738a45c..dd3b473 100644 --- a/molgroups/refl1d_interface/__init__.py +++ b/molgroups/refl1d_interface/__init__.py @@ -2,6 +2,10 @@ MolgroupsMixedExperiment, make_samples) +from molgroups.refl1d_interface.sas import (SASReflectivityModel, + SASReflectivityExperiment, + SASReflectivityMolgroupsExperiment) + from molgroups.refl1d_interface.layers import (MolgroupsStack, MolgroupsLayer) diff --git a/molgroups/refl1d_interface/examples/tiox_dopc/tiox_dopc_sas.py b/molgroups/refl1d_interface/examples/tiox_dopc/tiox_dopc_sas.py new file mode 100644 index 0000000..3b8618a --- /dev/null +++ b/molgroups/refl1d_interface/examples/tiox_dopc/tiox_dopc_sas.py @@ -0,0 +1,158 @@ +"""Example script using molgroups Refl1D interface objects""" + +import numpy as np +from refl1d.names import Parameter, SLD, Slab, FitProblem, load4 +from refl1d.probe.resolution import divergence +from molgroups import components as cmp +from molgroups.refl1d_interface import (SolidSupportedBilayer, + MolgroupsLayer, + MolgroupsStack, + MolgroupsExperiment) + +from molgroups.refl1d_interface import SASReflectivityMolgroupsExperiment, SASReflectivityModel + +## === Probes/data files === +probe_d2o = load4('ch061_d2o_ph7.refl', back_reflectivity=True, name='D2O') +probe_h2o = load4('ch060_h2o_ph7.refl', back_reflectivity=True, name='H2O') + +# Probe parameters +probes = [probe_d2o, probe_h2o] + +# Probe parameters +intensity = Parameter(name='intensity', value=0.8).range(0.65, 1.0) +sample_broadening = Parameter(name='sample broadening', value=0.0).range(-0.003, 0.02) +theta_offset = Parameter(name='theta offset', value=0.0).range(-0.02, 0.02) + +# apply background and intensity to all probes +for probe in probes: + probe.background.limits = (-np.inf, np.inf) + probe.background.range(-1e-6, 1e-5) + probe.intensity = intensity + + # if probes support these + probe.sample_broadening = sample_broadening + probe.theta_offset = theta_offset + +## === Structural parameters === + +vf_bilayer = Parameter(name='volume fraction bilayer', value=0.9).range(0.0, 1.0) +l_lipid1 = Parameter(name='inner acyl chain thickness', value=10.0).range(8, 30) +l_lipid2 = Parameter(name='outer acyl chain thickness', value=10.0).range(8, 18) +l_submembrane = Parameter(name='submembrane thickness', value=10.0).range(0, 50) +sigma = Parameter(name='bilayer roughness', value=5).range(0.5, 9) +global_rough = Parameter(name ='substrate roughness', value=5).range(2, 9) +tiox_rough = Parameter(name='titanium oxide roughness', value=4).range(2, 9) +d_oxide = Parameter(name='silicon oxide layer thickness', value=10).range(5, 30) +d_tiox = Parameter(name='titanium oxide layer thickness', value=110).range(100, 200) + +## === Materials === + +# Material definitions +d2o = SLD(name='d2o', rho=6.3000, irho=0.0000) +h2o = SLD(name='h2o', rho=-0.56, irho=0.0000) +tiox = SLD(name='tiox', rho=2.1630, irho=0.0000) +siox = SLD(name='siox', rho=4.1000, irho=0.0000) +silicon = SLD(name='silicon', rho=2.0690, irho=0.0000) + +# Material SLD parameters +d2o.rho.range(5.3000, 6.36) +h2o.rho.range(-0.56, 0.6) +tiox.rho.range(1.2, 3.2) +siox.rho.range(2.8, 4.8) + +## === Molecular groups === + +# overlap between substrate and molgroups layer +overlap = 30.0 + +# thickness of molgroups layer +thickness = 200.0 + +# define lipids and number fractions +DOPC = cmp.Lipid(name='DOPC', headgroup=cmp.pc, tails=2 * [cmp.oleoyl], methyls=[cmp.methyl]) +lipidlist = [DOPC] +lipid_nf = [1.0] + +def bilayer(substrate, contrast): + + blm = SolidSupportedBilayer(name='bilayer', + overlap=overlap, + lipids=lipidlist, + inner_lipid_nf=lipid_nf, + outer_lipid_nf=lipid_nf, + rho_substrate=tiox.rho, + l_siox=0.0, + vf_bilayer=vf_bilayer, + l_lipid1=l_lipid1, + l_lipid2=l_lipid2, + l_submembrane=l_submembrane, + substrate_rough=tiox_rough, + sigma=sigma) + + mollayer = MolgroupsLayer(base_group=blm, + thickness=thickness, + contrast=contrast, + name='bilayer layer ' + contrast.name) + + return MolgroupsStack(substrate=substrate, + molgroups_layer=mollayer, + name=mollayer.name) + +## == Sample layer stack == + +layer_silicon = Slab(material=silicon, thickness=0.0000, interface=global_rough) +layer_siox = Slab(material=siox, thickness=d_oxide, interface=global_rough) +layer_tiox = Slab(material=tiox, thickness=d_tiox - overlap, interface=0.00) + +substrate = layer_silicon | layer_siox | layer_tiox + +# Use the bilayer definition function to generate the bilayer SLD profile, passing in the relevant parameters. +sample_d2o, sample_h2o = [bilayer(substrate, contrast) for contrast in [d2o, h2o]] + +## === Problem definition === +## step = True corresponds to a calculation of the reflectivity from an actual profile +## with microslabbed interfaces. When step = False, the Nevot-Croce +## approximation is used to account for roughness. This approximation speeds up +## the calculation tremendously, and is reasonably accurate as long as the +## roughness is much less than the layer thickness +step = False +STEPSIZE=0.5 + +# calculate full transverse divergence (2 * FWHM) for MAGIK reflectometer +S1_transverse = 150.0 +S2_transverse = 25.0 +L2 = 330.0 +L1 = 1403.0 + 330.0 +dTl = 2 * np.ones_like(probe_d2o.Q) * divergence(0, (S1_transverse, S2_transverse), (L1, L2)) + +bilayer_thickness = Parameter(name='SANS bilayer thickness', value=40.0).range(20.0, 80.0) +bilayer_spacing = Parameter(name='SANS bilayer spacing', value=20.0).range(5.0, 50.0) +bilayer_sld = Parameter(name='SANS bilayer sld', value=0.4).range(-0.5, 0.5) +volume_fraction_bilayer = Parameter(name='SANS volume fraction', value=0.00001).range(0.0, 0.0001) +n_bilayers = Parameter(name='SANS number of bilayers', value=4).range(1, 20) + +sans_parameters = {'scale': 1.0, + 'background': 0.0, + 'sld': bilayer_sld, + 'thick_shell': bilayer_thickness, + 'thick_solvent': bilayer_spacing, + 'radius': 100.0, + 'volfraction': volume_fraction_bilayer, + 'n_shells': n_bilayers * 2 - 1 + } + +sasmodel_d = SASReflectivityModel(sas_model_name='multilayer_vesicle', + dtheta_l=dTl, + parameters=sans_parameters | {'sld_solvent': d2o.rho}) + +sasmodel_h = SASReflectivityModel(sas_model_name='multilayer_vesicle', + dtheta_l=dTl, + parameters=sans_parameters | {'sld_solvent': h2o.rho}) + + +model_d2o = SASReflectivityMolgroupsExperiment(sas_model=sasmodel_d, sample=sample_d2o, probe=probe_d2o, dz=STEPSIZE, step_interfaces = step) +model_h2o = SASReflectivityMolgroupsExperiment(sas_model=sasmodel_h, sample=sample_h2o, probe=probe_h2o, dz=STEPSIZE, step_interfaces = step) + +problem = FitProblem([model_d2o, model_h2o]) + +problem.name = "tiox_dopc_d2o_h2o" diff --git a/molgroups/refl1d_interface/sas.py b/molgroups/refl1d_interface/sas.py new file mode 100644 index 0000000..a6a279d --- /dev/null +++ b/molgroups/refl1d_interface/sas.py @@ -0,0 +1,233 @@ +""" Module for interfacing combined reflectivity and small angle scattering models +""" + +from dataclasses import dataclass + +import numpy as np +import plotly.graph_objs as go + +from bumps.parameter import Parameter +from bumps.webview.server.custom_plot import CustomWebviewPlot +from refl1d.experiment import Experiment +from refl1d.probe.resolution import dTdL2dQ, sigma2FWHM +from refl1d.webview.server.colors import COLORS + +from sasmodels.core import load_model +from sasmodels.direct_model import DirectModel +from sasmodels.data import Data1D + +from .experiment import MolgroupsExperiment + +@dataclass +class SASReflectivityModel: + """ Class to hold sasmodels model information + """ + sas_model_name: str | None = None + dtheta_l: float | None = None + parameters: dict[str, float | Parameter] | None = None + +class SASReflectivityMixin: + """ + Mixin class that adds SAS capabilities to ANY Refl1D Experiment. + It overrides reflectivity(), parameters(), and registers the SAS plot. + + Requires a probe object with Q, T, L, dL attributes, due to the + need to calculate resolution. + """ + + # Type hinting for the mixin (expects these to exist on the child) + sas_model: 'SASReflectivityModel' + _sasmodel: object + _cache: dict + probe: object + name: str + + def _init_sas(self, sas_model): + """ Helper to initialize SAS components. Call this from the child __init__ """ + self.sas_model = sas_model + + # Initialize Parameters + if sas_model is not None: + for k, p in self.sas_model.parameters.items(): + if not isinstance(p, Parameter): + self.sas_model.parameters[k] = Parameter.default(p, name=k) + + # Initialize Kernel + if sas_model.sas_model_name is not None: + self._sasmodel = load_model(sas_model.sas_model_name) + + # Register the Decomposition Plot + # Note: 'self' here will be the full Experiment instance + self.register_webview_plot( + plot_title='SAS/Refl Decomposition', + plot_function=sas_decomposition_plot, + change_with='parameter' + ) + + def parameters(self): + # Merge parent parameters (Experiment/Molgroups) with SAS parameters + return super().parameters() | {'sas': self.sas_model.parameters if self.sas_model is not None else {}} + + def sas(self): + """ Calculate the small angle scattering I(q) """ + key = ("small_angle_scattering") + if key not in self._cache: + data = Data1D(x=self.probe.Q) + + # calculate Q-transformed slit widths. dQ is assumed to be sigma, while dtheta_l is 2 * FWHM + data.dxl = dTdL2dQ(np.zeros_like(self.probe.T), self.sas_model.dtheta_l, self.probe.L, self.probe.dL) + data.dxw = 2 * sigma2FWHM(self.probe.dQ) + + pars = {k: float(p) for k, p in self.sas_model.parameters.items()} + + # Calculate + sasmodel = DirectModel(data=data, model=self._sasmodel) + Iq = sasmodel(**pars) + self._cache[key] = Iq + return self._cache[key] + + def reflectivity(self, resolution=True, interpolation=0): + # 1. Get base reflectivity (Calculated by Experiment or MolgroupsExperiment) + Q, Rq = super().reflectivity(resolution, interpolation) + + # 2. Add SAS signal + Iq = self.sas() if self.sas_model is not None else 0.0 + return Q, Rq + Iq + + +# --- 2. CONCRETE CLASSES --- + +@dataclass(init=False) +class SASReflectivityExperiment(SASReflectivityMixin, Experiment): + """ + Standard SAS + Reflectivity Experiment. + Inherits from Experiment. + """ + sas_model: 'SASReflectivityModel' = None + + def __init__(self, sas_model=None, sample=None, probe=None, name=None, **kwargs): + # 1. Initialize Parent (Experiment) + super().__init__(sample, probe, name, **kwargs) + + # 2. Initialize Mixin + self._init_sas(sas_model) + + +@dataclass(init=False) +class SASReflectivityMolgroupsExperiment(SASReflectivityMixin, MolgroupsExperiment): + """ + Molgroups-Enabled SAS + Reflectivity Experiment. + Inherits from MolgroupsExperiment. + """ + sas_model: 'SASReflectivityModel' = None + + def __init__(self, sas_model=None, sample=None, probe=None, name=None, **kwargs): + # 1. Initialize Parent (MolgroupsExperiment) + # This automatically registers the CVO plots, Table plots, etc. + super().__init__(sample, probe, name, **kwargs) + + # 2. Initialize Mixin + self._init_sas(sas_model) + +def sas_decomposition_plot(model: SASReflectivityExperiment, problem=None) -> CustomWebviewPlot: + """ + Webview plot that shows the decomposition of the signal into + Reflectivity (Rq) and SAS (Iq) components. + + Args: + model: The SASReflectivityExperiment instance (passed by webview) + problem: The Bumps FitProblem (passed by webview) + """ + + # 2. Calculate Components + # Total Theory = R(q) + I(q) + # We call reflectivity() which returns the sum + Q, total_theory = model.reflectivity() + + # SANS Component = I(q) + # We call sas() directly. Handle case where sas_model might be None + if model.sas_model is not None: + Iq = model.sas() + else: + Iq = np.zeros_like(Q) + + # Reflectivity Component = R(q) + # Derived by subtraction to ensure consistency + Rq = total_theory - Iq + + # 3. Get Data for comparison + data_y = model.probe.R + data_dy = model.probe.dR + + # 4. Construct Plotly Figure + fig = go.Figure() + + # -- Trace: Data -- + fig.add_trace(go.Scatter( + x=Q, y=data_y, + error_y=dict( + type='data', + array=data_dy, + visible=True, + color='rgba(0, 0, 0, 0.25)' # <--- Explicit opacity for error bars + ), + mode='markers', + name='Data', + marker=dict( + color='rgba(0, 0, 0, 0.25)', # <--- Explicit opacity for markers (0.25 = 25% visible) + size=6 + ) + )) + + # -- Trace: Total Theory -- + fig.add_trace(go.Scatter( + x=Q, y=total_theory, + mode='lines', + name='Total Model (R+I)', + line=dict(color=COLORS[0], width=3) + )) + + # -- Trace: Reflectivity Component -- + fig.add_trace(go.Scatter( + x=Q, y=Rq, + mode='lines', + name='Reflectivity R(q)', + line=dict(color=COLORS[1], width=2, dash='dash') + )) + + # -- Trace: SANS Component -- + fig.add_trace(go.Scatter( + x=Q, y=Iq, + mode='lines', + name='SANS I(q)', + line=dict(color=COLORS[2], width=2, dash='dot') + )) + + # 5. Styling + fig.update_layout( + title=f'Signal Decomposition: {model.name}', + xaxis_title='Q (Å⁻¹)', + xaxis_type='linear', + template='plotly_white', + yaxis=dict( + title='Intensity (R + I)', + type='log', + exponentformat='power', # <--- This forces 10^x notation + showexponent='all' # Ensures exponents are shown for all ticks + ), + legend=dict(x=0.01, y=0.01, xanchor='left', yanchor='bottom', bgcolor='rgba(255,255,255,0.8)') + ) + + # 6. Prepare CSV Export Data + # Simple CSV format: Q, Data, Error, Total, Rq, Iq + csv_header = "Q,R,dR,Theory,Rq,Iq\n" + csv_rows = [] + for i in range(len(Q)): + row = f"{float(Q[i]):.6e},{float(data_y[i]):.6e},{float(data_dy[i]):.6e},{float(total_theory[i]):.6e},{float(Rq[i]):.6e},{float(Iq[i]):.6e}" + csv_rows.append(row) + + export_data = csv_header + "\n".join(csv_rows) + + return CustomWebviewPlot(fig_type='plotly', + plotdata=fig, + exportdata=export_data) diff --git a/pyproject.toml b/pyproject.toml index 4740084..ba83de4 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -30,6 +30,7 @@ dependencies = [ [project.optional-dependencies] examples = ["sasdata", "pandas", "sasmodels"] refl1d = ["refl1d"] +sas = ["sasmodels"] [project.urls] repository = "https://github.com/reflectometry/molgroups"