From 68dc3d637dfeda12d9dc291c128378dd121db38b Mon Sep 17 00:00:00 2001 From: Artur Glavic Date: Wed, 13 May 2026 12:30:18 +0200 Subject: [PATCH 01/22] Allow new curly bracket syntax in stack for SubStackType objects --- orsopy/fileio/model_language.py | 44 ++++++++++++++++++++++++++++++--- 1 file changed, 41 insertions(+), 3 deletions(-) diff --git a/orsopy/fileio/model_language.py b/orsopy/fileio/model_language.py index e859d62..4880a88 100644 --- a/orsopy/fileio/model_language.py +++ b/orsopy/fileio/model_language.py @@ -10,6 +10,8 @@ from dataclasses import dataclass from typing import Dict, List, Optional, Union +import yaml + from ..utils.chemical_formula import Formula from . import model_complex from .base import Header, Literal @@ -26,13 +28,13 @@ def find_idx(string, start, value): return next_idx -def find_closing(string, start): +def find_closing(string, start, brackets="()"): open_brackets = 1 idx = start while idx < len(string): - if string[idx] == "(": + if string[idx] == brackets[0]: open_brackets += 1 - if string[idx] == ")": + if string[idx] == brackets[1]: open_brackets -= 1 if open_brackets == 0: return idx @@ -87,6 +89,24 @@ def resolve_names(self, resolvable_items): # if there is a higher level envirnment, it is kept if not overwritten environment = self.environment obj = SubStack(repetitions=rep, stack=sub_stack.strip(), environment=environment) + elif "{" in stack[idx:next_idx]: + # allow supplying SubStackType classes to be defined in the stack string directly + # with the syntac ClassName{p1: v1, p2: v2} + open_idx = find_idx(stack, idx, "{") + close_idx = find_closing(stack, open_idx + 1, brackets="{}") + class_name = stack[idx:open_idx].strip() + for T in SubStackType.__subclasses__(): + if T.sub_stack_class == class_name: + ssclass = T + break + obj = ssclass.from_dict(yaml.load(stack[open_idx : close_idx + 1], Loader=yaml.FullLoader)) + try: + thickness = float(stack[close_idx + 1 : next_idx].strip()) + except ValueError: + pass + else: + if getattr(obj, "thickness", "ignore") is None: + obj.thickness = thickness else: items = stack[idx:next_idx].strip().rsplit(None, 1) item = items[0].strip() @@ -279,6 +299,24 @@ def resolve_stack(self): else: environment = None obj = SubStack(repetitions=rep, stack=sub_stack.strip(), environment=environment) + elif "{" in stack[idx:next_idx]: + # allow supplying SubStackType classes to be defined in the stack string directly + # with the syntac ClassName{p1: v1, p2: v2} + open_idx = find_idx(stack, idx, "{") + close_idx = find_closing(stack, open_idx + 1, brackets="{}") + class_name = stack[idx:open_idx].strip() + for T in SubStackType.__subclasses__(): + if T.sub_stack_class == class_name: + ssclass = T + break + obj = ssclass.from_dict(yaml.load(stack[open_idx : close_idx + 1], Loader=yaml.FullLoader)) + try: + thickness = float(stack[close_idx + 1 : next_idx].strip()) + except ValueError: + pass + else: + if getattr(obj, "thickness", "ignore") is None: + obj.thickness = thickness else: items = stack[idx:next_idx].strip().rsplit(None, 1) item = items[0].strip() From 7ae893940885756cc9cdbc6743285804d2a017b4 Mon Sep 17 00:00:00 2001 From: Artur Glavic Date: Mon, 18 May 2026 12:47:30 +0200 Subject: [PATCH 02/22] Start creating class for lipid Bilayer using build-in lipids or two materials, still missing correct layer calculation --- orsopy/fileio/model_complex.py | 80 +++++++++++++++++++++++++++++++++- 1 file changed, 78 insertions(+), 2 deletions(-) diff --git a/orsopy/fileio/model_complex.py b/orsopy/fileio/model_complex.py index b82cccc..0ef7c0b 100644 --- a/orsopy/fileio/model_complex.py +++ b/orsopy/fileio/model_complex.py @@ -5,10 +5,12 @@ have a common "sub_stack_class" attribute that has to be set to the class name. """ -from dataclasses import dataclass +from dataclasses import dataclass, field from typing import List, Optional, Union -from .base import ComplexValue, Header, Literal, Value +from ..slddb.constants import u2g +from ..utils.chemical_formula import Formula +from .base import Header, Literal, Value from .model_building_blocks import SPECIAL_MATERIALS, Composit, Layer, Material, ModelParameters, SubStackType @@ -108,3 +110,77 @@ def resolve_to_layers(self) -> List[Layer]: output.append(Layer(material=composition, thickness=thickness, roughness=roughness)) output[0].roughness = self.roughness return output + + +@dataclass +class Bilayer(Header, SubStackType): + """ + Building block corresponding to a bilayer of lipids with outer (heads) and inner (tails) + molecule definition. The bilayer is calculated from molecular volume, area per molecule (apm) and + the level of hydration. + The solvent used is either the environment set by a SubStack above or the default_solvent attribute + of the ModelParameters. + """ + + _known_lipids = { + "DMPC": ( + Material(formula="C10H18O8NP", number_density=Value(1.0 / 3.19, "nm^2")), + Material(formula="C26H54", number_density=Value(1.0 / 7.82, "nm^2")), + ), + } + + lipid: Optional[Literal["DMPC"]] = None + outer: Optional[Union[Material, str]] = None + inner: Optional[Union[Material, str]] = None + apm: Optional[Union[float, Value]] = field(default_factory=lambda: Value(0.7, unit="nm^2")) + outer_hydration: Optional[float] = 0.3 + inner_hydration: Optional[float] = 0.3 + outer_hydration_2: Optional[float] = None + inner_hydration_2: Optional[float] = None + sub_stack_class: Literal["Bilayer"] = "Bilayer" + + _environment = None + + def resolve_names(self, resolvable_items): + self._materials = [] + for i, mi in enumerate([self.outer, self.inner]): + if isinstance(mi, Material): + material = mi + elif mi in resolvable_items: + material = resolvable_items[mi] + elif mi in SPECIAL_MATERIALS: + material = SPECIAL_MATERIALS[mi] + else: + material = Material(formula=mi) + self._materials.append(material) + if "environment" in resolvable_items: + self._environment = resolvable_items["environment"] + if self._environment in resolvable_items: + self._environment = resolvable_items[self._environment] + elif self._environment in SPECIAL_MATERIALS: + self._environment = SPECIAL_MATERIALS[self._environment] + elif isinstance(self._environment, str): + self._environment = Material(formula=self._environment) + + def resolve_defaults(self, defaults: ModelParameters) -> None: + if not self._environment: + self._environment = defaults.default_solvent + else: + self._environment.resolve_defaults(defaults) + + def resolve_to_layers(self) -> List[Layer]: + for material in [self._environment] + self._materials: + material.generate_density() + if material.number_density is None: + # need to generate number density from mass density and formula + formula = Formula(material.formula, strict=True) + fu_mass = 0.0 + for element, number in formula.elements: + if element.mass is None: + raise ValueError(f"No mass known for element {element}") + fu_mass += number * element.mass + material.number_density = Value( + material.mass_density.as_unit("g/cm^3") / fu_mass / u2g * 1e21, unit="1/nm^3" + ) + + return [] From f353dae6c42697e3b6bdaa7194f6403ab9055c29 Mon Sep 17 00:00:00 2001 From: Artur Glavic Date: Mon, 18 May 2026 13:48:55 +0200 Subject: [PATCH 03/22] Minor improvements for resolution of default_solvent and lipid name implementation --- orsopy/fileio/model_building_blocks.py | 35 +++++++++++++++----------- orsopy/fileio/model_complex.py | 23 ++++++++++------- 2 files changed, 34 insertions(+), 24 deletions(-) diff --git a/orsopy/fileio/model_building_blocks.py b/orsopy/fileio/model_building_blocks.py index f4e431b..9ad00a0 100644 --- a/orsopy/fileio/model_building_blocks.py +++ b/orsopy/fileio/model_building_blocks.py @@ -147,20 +147,6 @@ def get_sld(self, xray_energy=None) -> complex: return 0.0j -@dataclass -class ModelParameters(Header): - roughness: Optional[Value] = field(default_factory=lambda: Value(0.3, "nm")) - length_unit: Optional[str] = "nm" - mass_density_unit: Optional[str] = "g/cm^3" - number_density_unit: Optional[str] = "1/nm^3" - sld_unit: Optional[str] = "1/angstrom^2" - magnetic_moment_unit: Optional[str] = "muB" - slice_resolution: Optional[Value] = field(default_factory=lambda: Value(1.0, "nm")) - default_solvent: Optional[Material] = field( - default_factory=lambda: Material(formula="H2O", mass_density=Value(1.0, "g/cm^3")) - ) - - @dataclass class Composit(Header): composition: Dict[str, float] @@ -184,7 +170,7 @@ def resolve_names(self, resolvable_items): self._composition_materials[key] = material - def resolve_defaults(self, defaults: ModelParameters): + def resolve_defaults(self, defaults: "ModelParameters"): for mat in self._composition_materials.values(): mat.resolve_defaults(defaults) @@ -209,6 +195,25 @@ def get_sld(self, xray_energy=None): return material.get_sld(xray_energy=xray_energy) +@dataclass +class ModelParameters(Header): + roughness: Optional[Value] = field(default_factory=lambda: Value(0.3, "nm")) + length_unit: Optional[str] = "nm" + mass_density_unit: Optional[str] = "g/cm^3" + number_density_unit: Optional[str] = "1/nm^3" + sld_unit: Optional[str] = "1/angstrom^2" + magnetic_moment_unit: Optional[str] = "muB" + slice_resolution: Optional[Value] = field(default_factory=lambda: Value(1.0, "nm")) + default_solvent: Optional[Union[Material, Composit, str]] = field( + default_factory=lambda: Material(formula="H2O", mass_density=Value(1.0, "g/cm^3")) + ) + + def __post_init__(self): + super().__post_init__() + if isinstance(self.default_solvent, str): + self.default_solvent = Material(formula=self.default_solvent) + + SPECIAL_MATERIALS = { "vacuum": Material(sld=ComplexValue(real=0.0, imag=0.0, unit="1/angstrom^2")), "air": Material(formula="N8O2", mass_density=Value(1.225, unit="kg/m^3")), diff --git a/orsopy/fileio/model_complex.py b/orsopy/fileio/model_complex.py index 0ef7c0b..a8862cb 100644 --- a/orsopy/fileio/model_complex.py +++ b/orsopy/fileio/model_complex.py @@ -124,8 +124,8 @@ class Bilayer(Header, SubStackType): _known_lipids = { "DMPC": ( - Material(formula="C10H18O8NP", number_density=Value(1.0 / 3.19, "nm^2")), - Material(formula="C26H54", number_density=Value(1.0 / 7.82, "nm^2")), + Material(formula="C10H18O8NP", number_density=Value(1.0 / 3.19, "1/nm^2")), + Material(formula="C26H54", number_density=Value(1.0 / 7.82, "1/nm^2")), ), } @@ -142,8 +142,12 @@ class Bilayer(Header, SubStackType): _environment = None def resolve_names(self, resolvable_items): + if self.lipid is None: + oi_mats = (self.outer, self.inner) + else: + oi_mats = self._known_lipids[self.lipid] self._materials = [] - for i, mi in enumerate([self.outer, self.inner]): + for i, mi in enumerate(oi_mats): if isinstance(mi, Material): material = mi elif mi in resolvable_items: @@ -155,12 +159,13 @@ def resolve_names(self, resolvable_items): self._materials.append(material) if "environment" in resolvable_items: self._environment = resolvable_items["environment"] - if self._environment in resolvable_items: - self._environment = resolvable_items[self._environment] - elif self._environment in SPECIAL_MATERIALS: - self._environment = SPECIAL_MATERIALS[self._environment] - elif isinstance(self._environment, str): - self._environment = Material(formula=self._environment) + if not isinstance(self._environment, Material): + if self._environment in resolvable_items: + self._environment = resolvable_items[self._environment] + elif self._environment in SPECIAL_MATERIALS: + self._environment = SPECIAL_MATERIALS[self._environment] + elif isinstance(self._environment, str): + self._environment = Material(formula=self._environment) def resolve_defaults(self, defaults: ModelParameters) -> None: if not self._environment: From a186df6c57f56163b18df444e548af182f90a6b9 Mon Sep 17 00:00:00 2001 From: Artur Glavic Date: Mon, 18 May 2026 13:52:31 +0200 Subject: [PATCH 04/22] default apm unit resolution --- orsopy/fileio/model_complex.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/orsopy/fileio/model_complex.py b/orsopy/fileio/model_complex.py index a8862cb..d4d7906 100644 --- a/orsopy/fileio/model_complex.py +++ b/orsopy/fileio/model_complex.py @@ -172,6 +172,10 @@ def resolve_defaults(self, defaults: ModelParameters) -> None: self._environment = defaults.default_solvent else: self._environment.resolve_defaults(defaults) + if isinstance(self.apm, Value) and self.apm.unit is None: + self.apm.unit = defaults.length_unit + "^2" + elif not isinstance(self.apm, Value): + self.apm = Value(self.apm, unit=defaults.length_unit + "^2") def resolve_to_layers(self) -> List[Layer]: for material in [self._environment] + self._materials: From ad78a81a527baaab7684b5a10e4e619713dc2179 Mon Sep 17 00:00:00 2001 From: Artur Glavic Date: Tue, 19 May 2026 10:37:18 +0200 Subject: [PATCH 05/22] Add example for Bilayer and first version resolving to layers with thickness and density, math needs checking --- examples/resolve_and_plot_model.py | 2 +- examples/sample_model_example_10.yml | 33 +++++++++++++++++++++++ orsopy/fileio/model_complex.py | 39 +++++++++++++++++++++++++--- 3 files changed, 70 insertions(+), 4 deletions(-) create mode 100644 examples/sample_model_example_10.yml diff --git a/examples/resolve_and_plot_model.py b/examples/resolve_and_plot_model.py index 11b62c4..1815b88 100644 --- a/examples/resolve_and_plot_model.py +++ b/examples/resolve_and_plot_model.py @@ -32,7 +32,7 @@ def main(txt=None): dtxt = dtxt["data_source"] if "sample" in dtxt: dtxt = dtxt["sample"]["model"] - sample = model_language.SampleModel(**dtxt) + sample = model_language.SampleModel.from_dict(dtxt) txt += f"\n{sample.stack}" else: sample = model_language.SampleModel(stack=txt) diff --git a/examples/sample_model_example_10.yml b/examples/sample_model_example_10.yml new file mode 100644 index 0000000..04fb329 --- /dev/null +++ b/examples/sample_model_example_10.yml @@ -0,0 +1,33 @@ +sample: + model: + stack: "Si | Bilayer{lipid: DMPC} | lbl2 | (lbl3) in H2O | H2O" + sub_stacks: + lbl2: + sub_stack_class: Bilayer + lipid: DMPC + apm: {magnitude: 0.8, unit: nm^2} + outer_hydration: 0.3 + inner_hydration: 0.3 + inner_hydration_2: 0.4 + outer_hydration_2: 0.5 + lbl3: + sub_stack_class: Bilayer + outer: + formula: C10H18O8NP + number_density: 0.3135 # 1/nm^2 + inner: + formula: C26H54 + number_density: 0.128 # 1/nm^2 + apm: 0.9 # nm^2 + outer_hydration: 0.3 + inner_hydration: 0.3 + inner_hydration_2: 0.4 + outer_hydration_2: 0.5 + globals: + default_solvent: D2O + + # below is just used for orsopy testing and not part of this model + test: + stacks: 5 + blocks: 5 + layers: 14 diff --git a/orsopy/fileio/model_complex.py b/orsopy/fileio/model_complex.py index d4d7906..cb19766 100644 --- a/orsopy/fileio/model_complex.py +++ b/orsopy/fileio/model_complex.py @@ -124,8 +124,8 @@ class Bilayer(Header, SubStackType): _known_lipids = { "DMPC": ( - Material(formula="C10H18O8NP", number_density=Value(1.0 / 3.19, "1/nm^2")), - Material(formula="C26H54", number_density=Value(1.0 / 7.82, "1/nm^2")), + Material(formula="C10H18O8NP", number_density=Value(1.0 / 3.19, "1/nm^3")), + Material(formula="C26H54", number_density=Value(1.0 / 7.82, "1/nm^3")), ), } @@ -137,6 +137,7 @@ class Bilayer(Header, SubStackType): inner_hydration: Optional[float] = 0.3 outer_hydration_2: Optional[float] = None inner_hydration_2: Optional[float] = None + roughness: Optional[Union[float, Value]] = None sub_stack_class: Literal["Bilayer"] = "Bilayer" _environment = None @@ -177,6 +178,22 @@ def resolve_defaults(self, defaults: ModelParameters) -> None: elif not isinstance(self.apm, Value): self.apm = Value(self.apm, unit=defaults.length_unit + "^2") + for mi in self._materials: + mi.resolve_defaults(defaults) + + if self.roughness is None: + self.roughness = defaults.roughness + elif not isinstance(self.roughness, Value): + self.roughness = Value(self.roughness, unit=defaults.length_unit) + elif self.roughness.unit is None: + self.roughness.unit = defaults.length_unit + + @staticmethod + def mixed_material(material: Material, solvent: Material, fraction: float): + srdens = solvent.number_density.magnitude / material.number_density.as_unit(solvent.number_density.unit) + FU = f"({material.formula}){1-fraction}({solvent.formula}){srdens*fraction}" + return Material(formula=FU, number_density=material.number_density) + def resolve_to_layers(self) -> List[Layer]: for material in [self._environment] + self._materials: material.generate_density() @@ -192,4 +209,20 @@ def resolve_to_layers(self) -> List[Layer]: material.mass_density.as_unit("g/cm^3") / fu_mass / u2g * 1e21, unit="1/nm^3" ) - return [] + d_head = Value(1.0 / (self._materials[0].number_density.as_unit("1/nm^3") * self.apm.as_unit("nm^2")), "nm") + d_tail = Value(1.0 / (self._materials[1].number_density.as_unit("1/nm^3") * self.apm.as_unit("nm^2")), "nm") + + m_head_1 = self.mixed_material(self._materials[0], self._environment, self.outer_hydration) + m_head_2 = self.mixed_material( + self._materials[0], self._environment, self.outer_hydration_2 or self.outer_hydration + ) + m_tail_1 = self.mixed_material(self._materials[1], self._environment, self.inner_hydration) + m_tail_2 = self.mixed_material( + self._materials[1], self._environment, self.inner_hydration_2 or self.inner_hydration + ) + + head = Layer(thickness=d_head, material=m_head_1, roughness=self.roughness) + tail = Layer(thickness=d_tail, material=m_tail_1, roughness=self.roughness) + head_2 = Layer(thickness=d_tail, material=m_tail_2, roughness=self.roughness) + tail_2 = Layer(thickness=d_head, material=m_head_2, roughness=self.roughness) + return [head, tail, tail_2, head_2] From 6639e24af2a0aa742b7a34c355116970d755a58d Mon Sep 17 00:00:00 2001 From: Artur Glavic Date: Tue, 19 May 2026 15:03:12 +0200 Subject: [PATCH 06/22] Add volume property to Material if number_density is already defined --- orsopy/fileio/model_building_blocks.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/orsopy/fileio/model_building_blocks.py b/orsopy/fileio/model_building_blocks.py index 9ad00a0..c6f7634 100644 --- a/orsopy/fileio/model_building_blocks.py +++ b/orsopy/fileio/model_building_blocks.py @@ -114,6 +114,19 @@ def generate_density(self): self.number_density = Value(magnitude=0.0, unit="1/nm^3") self.comment = "could not locate density information for material" + @property + def volume(self): + # returns the FU volume from the density information + self.generate_density() + if self.number_density is None: + return None + else: + from .base import get_unit_registry + + unit_registry = get_unit_registry() + val = 1.0 / (self.number_density.magnitude * unit_registry(self.number_density.unit)) + return Value(magnitude=val.magnitude, unit=f"{val.units:~}") + def get_sld(self, xray_energy=None) -> complex: if self.relative_density is None: rel = 1.0 From 90dd8d4494a5dd56073626eaaf8c8f443ad0cb5d Mon Sep 17 00:00:00 2001 From: Artur Glavic Date: Tue, 19 May 2026 15:12:32 +0200 Subject: [PATCH 07/22] Add coverage parameter --- examples/sample_model_example_10.yml | 1 + orsopy/fileio/model_complex.py | 7 ++++--- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/examples/sample_model_example_10.yml b/examples/sample_model_example_10.yml index 04fb329..82c39a6 100644 --- a/examples/sample_model_example_10.yml +++ b/examples/sample_model_example_10.yml @@ -23,6 +23,7 @@ sample: inner_hydration: 0.3 inner_hydration_2: 0.4 outer_hydration_2: 0.5 + coverage: 0.75 globals: default_solvent: D2O diff --git a/orsopy/fileio/model_complex.py b/orsopy/fileio/model_complex.py index cb19766..309ea66 100644 --- a/orsopy/fileio/model_complex.py +++ b/orsopy/fileio/model_complex.py @@ -137,6 +137,7 @@ class Bilayer(Header, SubStackType): inner_hydration: Optional[float] = 0.3 outer_hydration_2: Optional[float] = None inner_hydration_2: Optional[float] = None + coverage: Optional[float] = 1.0 roughness: Optional[Union[float, Value]] = None sub_stack_class: Literal["Bilayer"] = "Bilayer" @@ -188,10 +189,10 @@ def resolve_defaults(self, defaults: ModelParameters) -> None: elif self.roughness.unit is None: self.roughness.unit = defaults.length_unit - @staticmethod - def mixed_material(material: Material, solvent: Material, fraction: float): + def mixed_material(self, material: Material, solvent: Material, hydtration: float): + mat_frac = (1.0 - hydtration) * self.coverage srdens = solvent.number_density.magnitude / material.number_density.as_unit(solvent.number_density.unit) - FU = f"({material.formula}){1-fraction}({solvent.formula}){srdens*fraction}" + FU = f"({material.formula}){mat_frac}({solvent.formula}){srdens*(1.-mat_frac)}" return Material(formula=FU, number_density=material.number_density) def resolve_to_layers(self) -> List[Layer]: From 7ab6f3ed4dfddd5a2f87d36f8a618c492cf19ab0 Mon Sep 17 00:00:00 2001 From: Artur Glavic Date: Tue, 19 May 2026 15:33:50 +0200 Subject: [PATCH 08/22] minor fix --- examples/resolve_and_plot_model.py | 1 + orsopy/fileio/model_complex.py | 4 ++-- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/examples/resolve_and_plot_model.py b/examples/resolve_and_plot_model.py index 1815b88..b3ce5db 100644 --- a/examples/resolve_and_plot_model.py +++ b/examples/resolve_and_plot_model.py @@ -39,6 +39,7 @@ def main(txt=None): # initial model before resolving any names print(repr(sample), "\n") print("\n".join([repr(ss) for ss in sample.resolve_stack()]), "\n") + print("\n".join([repr(ss) for ss in sample.resolve_to_blocks()]), "\n") layers = sample.resolve_to_layers() print("\n".join([repr(li) for li in layers])) diff --git a/orsopy/fileio/model_complex.py b/orsopy/fileio/model_complex.py index 309ea66..566b311 100644 --- a/orsopy/fileio/model_complex.py +++ b/orsopy/fileio/model_complex.py @@ -224,6 +224,6 @@ def resolve_to_layers(self) -> List[Layer]: head = Layer(thickness=d_head, material=m_head_1, roughness=self.roughness) tail = Layer(thickness=d_tail, material=m_tail_1, roughness=self.roughness) - head_2 = Layer(thickness=d_tail, material=m_tail_2, roughness=self.roughness) - tail_2 = Layer(thickness=d_head, material=m_head_2, roughness=self.roughness) + tail_2 = Layer(thickness=d_tail, material=m_tail_2, roughness=self.roughness) + head_2 = Layer(thickness=d_head, material=m_head_2, roughness=self.roughness) return [head, tail, tail_2, head_2] From e7931223405b4f361baeb99d3f0f02bae2f40a17 Mon Sep 17 00:00:00 2001 From: Artur Glavic Date: Wed, 3 Jun 2026 11:09:07 +0200 Subject: [PATCH 09/22] Remove option for lipid to be explicit, add Leaflet class --- examples/sample_model_example_10.yml | 19 +++- orsopy/fileio/model_complex.py | 133 +++++++++++++++++++++++---- 2 files changed, 130 insertions(+), 22 deletions(-) diff --git a/examples/sample_model_example_10.yml b/examples/sample_model_example_10.yml index 82c39a6..8a08e0e 100644 --- a/examples/sample_model_example_10.yml +++ b/examples/sample_model_example_10.yml @@ -1,10 +1,16 @@ sample: model: - stack: "Si | Bilayer{lipid: DMPC} | lbl2 | (lbl3) in H2O | H2O" + stack: "Si | DMPC | Bilayer{outer: PC, inner: DM} | biDMPC | (lbl3) in H2O | H2O" sub_stacks: - lbl2: + DMPC: + sub_stack_class: Leaflet + heads: PC + tails: DM + heads_first: false + biDMPC: sub_stack_class: Bilayer - lipid: DMPC + outer: PC + inner: DM apm: {magnitude: 0.8, unit: nm^2} outer_hydration: 0.3 inner_hydration: 0.3 @@ -24,6 +30,13 @@ sample: inner_hydration_2: 0.4 outer_hydration_2: 0.5 coverage: 0.75 + materials: + PC: + formula: C10H18O8NP + number_density: 0.3135 + DM: + formula: C26H54 + number_density: 0.128 globals: default_solvent: D2O diff --git a/orsopy/fileio/model_complex.py b/orsopy/fileio/model_complex.py index 566b311..989a87a 100644 --- a/orsopy/fileio/model_complex.py +++ b/orsopy/fileio/model_complex.py @@ -112,6 +112,113 @@ def resolve_to_layers(self) -> List[Layer]: return output +def mix_hydrate_material(material: Material, solvent: Material, hydtration: float, coverage: float = 1.0): + mat_frac = (1.0 - hydtration) * coverage + srdens = solvent.number_density.magnitude / material.number_density.as_unit(solvent.number_density.unit) + FU = f"({material.formula}){mat_frac}({solvent.formula}){srdens*(1.-mat_frac)}" + return Material(formula=FU, number_density=material.number_density) + + +@dataclass +class Leaflet(Header, SubStackType): + """ + Building block corresponding to a single layer of lipids with heads and tails + molecule definition. The layer is calculated from molecular volume, area per molecule (apm) and + the level of hydration. + The solvent used is either the environment set by a SubStack above or the default_solvent attribute + of the ModelParameters. + The order of heads/tails with respect to the beam side can be flipped using heads_first=False. + """ + + heads: Union[Material, str] + tails: Union[Material, str] + heads_first: Optional[bool] = True + apm: Optional[Union[float, Value]] = field(default_factory=lambda: Value(0.7, unit="nm^2")) + heads_hydration: Optional[float] = 0.3 + tails_hydration: Optional[float] = 0.3 + coverage: Optional[float] = 1.0 + roughness: Optional[Union[float, Value]] = None + sub_stack_class: Literal["Leaflet"] = "Leaflet" + + _environment = None + + def resolve_names(self, resolvable_items): + self._materials = [] + for i, mi in enumerate((self.heads, self.tails)): + if isinstance(mi, Material): + material = mi + elif mi in resolvable_items: + material = resolvable_items[mi] + material.original_name = mi + elif mi in SPECIAL_MATERIALS: + material = SPECIAL_MATERIALS[mi] + material.original_name = mi + else: + material = Material(formula=mi) + material.original_name = mi + self._materials.append(material) + if "environment" in resolvable_items: + self._environment = resolvable_items["environment"] + if not isinstance(self._environment, Material): + if self._environment in resolvable_items: + self._environment = resolvable_items[self._environment] + elif self._environment in SPECIAL_MATERIALS: + self._environment = SPECIAL_MATERIALS[self._environment] + elif isinstance(self._environment, str): + self._environment = Material(formula=self._environment) + + def resolve_defaults(self, defaults: ModelParameters) -> None: + if not self._environment: + self._environment = defaults.default_solvent + else: + self._environment.resolve_defaults(defaults) + if isinstance(self.apm, Value) and self.apm.unit is None: + self.apm.unit = defaults.length_unit + "^2" + elif not isinstance(self.apm, Value): + self.apm = Value(self.apm, unit=defaults.length_unit + "^2") + + for mi in self._materials: + mi.resolve_defaults(defaults) + + if self.roughness is None: + self.roughness = defaults.roughness + elif not isinstance(self.roughness, Value): + self.roughness = Value(self.roughness, unit=defaults.length_unit) + elif self.roughness.unit is None: + self.roughness.unit = defaults.length_unit + + def mixed_material(self, material: Material, solvent: Material, hydtration: float): + return mix_hydrate_material(material, solvent, hydtration, self.coverage) + + def resolve_to_layers(self) -> List[Layer]: + for material in [self._environment] + self._materials: + material.generate_density() + if material.number_density is None: + # need to generate number density from mass density and formula + formula = Formula(material.formula, strict=True) + fu_mass = 0.0 + for element, number in formula.elements: + if element.mass is None: + raise ValueError(f"No mass known for element {element}") + fu_mass += number * element.mass + material.number_density = Value( + material.mass_density.as_unit("g/cm^3") / fu_mass / u2g * 1e21, unit="1/nm^3" + ) + + d_head = Value(1.0 / (self._materials[0].number_density.as_unit("1/nm^3") * self.apm.as_unit("nm^2")), "nm") + d_tail = Value(1.0 / (self._materials[1].number_density.as_unit("1/nm^3") * self.apm.as_unit("nm^2")), "nm") + + m_head = self.mixed_material(self._materials[0], self._environment, self.heads_hydration) + m_tail = self.mixed_material(self._materials[1], self._environment, self.tails_hydration) + + head = Layer(thickness=d_head, material=m_head, roughness=self.roughness) + tail = Layer(thickness=d_tail, material=m_tail, roughness=self.roughness) + if self.heads_first: + return [head, tail] + else: + return [tail, head] + + @dataclass class Bilayer(Header, SubStackType): """ @@ -122,16 +229,8 @@ class Bilayer(Header, SubStackType): of the ModelParameters. """ - _known_lipids = { - "DMPC": ( - Material(formula="C10H18O8NP", number_density=Value(1.0 / 3.19, "1/nm^3")), - Material(formula="C26H54", number_density=Value(1.0 / 7.82, "1/nm^3")), - ), - } - - lipid: Optional[Literal["DMPC"]] = None - outer: Optional[Union[Material, str]] = None - inner: Optional[Union[Material, str]] = None + outer: Union[Material, str] + inner: Union[Material, str] apm: Optional[Union[float, Value]] = field(default_factory=lambda: Value(0.7, unit="nm^2")) outer_hydration: Optional[float] = 0.3 inner_hydration: Optional[float] = 0.3 @@ -144,20 +243,19 @@ class Bilayer(Header, SubStackType): _environment = None def resolve_names(self, resolvable_items): - if self.lipid is None: - oi_mats = (self.outer, self.inner) - else: - oi_mats = self._known_lipids[self.lipid] self._materials = [] - for i, mi in enumerate(oi_mats): + for i, mi in enumerate((self.outer, self.inner)): if isinstance(mi, Material): material = mi elif mi in resolvable_items: material = resolvable_items[mi] + material.original_name = mi elif mi in SPECIAL_MATERIALS: material = SPECIAL_MATERIALS[mi] + material.original_name = mi else: material = Material(formula=mi) + material.original_name = mi self._materials.append(material) if "environment" in resolvable_items: self._environment = resolvable_items["environment"] @@ -190,10 +288,7 @@ def resolve_defaults(self, defaults: ModelParameters) -> None: self.roughness.unit = defaults.length_unit def mixed_material(self, material: Material, solvent: Material, hydtration: float): - mat_frac = (1.0 - hydtration) * self.coverage - srdens = solvent.number_density.magnitude / material.number_density.as_unit(solvent.number_density.unit) - FU = f"({material.formula}){mat_frac}({solvent.formula}){srdens*(1.-mat_frac)}" - return Material(formula=FU, number_density=material.number_density) + return mix_hydrate_material(material, solvent, hydtration, self.coverage) def resolve_to_layers(self) -> List[Layer]: for material in [self._environment] + self._materials: From 439022d4558a0d48bde9a95001f617e86147c617 Mon Sep 17 00:00:00 2001 From: Artur Glavic Date: Wed, 3 Jun 2026 11:24:00 +0200 Subject: [PATCH 10/22] Make sure that resolve_to_blocks includes all density information needed --- orsopy/fileio/model_complex.py | 72 ++++++++++++++++++++++------------ 1 file changed, 46 insertions(+), 26 deletions(-) diff --git a/orsopy/fileio/model_complex.py b/orsopy/fileio/model_complex.py index 989a87a..4ea453c 100644 --- a/orsopy/fileio/model_complex.py +++ b/orsopy/fileio/model_complex.py @@ -190,20 +190,15 @@ def resolve_defaults(self, defaults: ModelParameters) -> None: def mixed_material(self, material: Material, solvent: Material, hydtration: float): return mix_hydrate_material(material, solvent, hydtration, self.coverage) + def resolve_to_blocks(self) -> List[Union["Layer", "SubStackType"]]: + # Make sure the block includes full material data + self.ensure_densities() + self.heads = self._materials[0] + self.tails = self._materials[1] + return [self] + def resolve_to_layers(self) -> List[Layer]: - for material in [self._environment] + self._materials: - material.generate_density() - if material.number_density is None: - # need to generate number density from mass density and formula - formula = Formula(material.formula, strict=True) - fu_mass = 0.0 - for element, number in formula.elements: - if element.mass is None: - raise ValueError(f"No mass known for element {element}") - fu_mass += number * element.mass - material.number_density = Value( - material.mass_density.as_unit("g/cm^3") / fu_mass / u2g * 1e21, unit="1/nm^3" - ) + self.ensure_densities() d_head = Value(1.0 / (self._materials[0].number_density.as_unit("1/nm^3") * self.apm.as_unit("nm^2")), "nm") d_tail = Value(1.0 / (self._materials[1].number_density.as_unit("1/nm^3") * self.apm.as_unit("nm^2")), "nm") @@ -218,6 +213,21 @@ def resolve_to_layers(self) -> List[Layer]: else: return [tail, head] + def ensure_densities(self): + for material in [self._environment] + self._materials: + material.generate_density() + if material.number_density is None: + # need to generate number density from mass density and formula + formula = Formula(material.formula, strict=True) + fu_mass = 0.0 + for element, number in formula.elements: + if element.mass is None: + raise ValueError(f"No mass known for element {element}") + fu_mass += number * element.mass + material.number_density = Value( + material.mass_density.as_unit("g/cm^3") / fu_mass / u2g * 1e21, unit="1/nm^3" + ) + @dataclass class Bilayer(Header, SubStackType): @@ -290,20 +300,15 @@ def resolve_defaults(self, defaults: ModelParameters) -> None: def mixed_material(self, material: Material, solvent: Material, hydtration: float): return mix_hydrate_material(material, solvent, hydtration, self.coverage) + def resolve_to_blocks(self) -> List[Union["Layer", "SubStackType"]]: + self.ensure_densities() + # Make sure the block includes full material data + self.inner = self._materials[1] + self.outer = self._materials[0] + return [self] + def resolve_to_layers(self) -> List[Layer]: - for material in [self._environment] + self._materials: - material.generate_density() - if material.number_density is None: - # need to generate number density from mass density and formula - formula = Formula(material.formula, strict=True) - fu_mass = 0.0 - for element, number in formula.elements: - if element.mass is None: - raise ValueError(f"No mass known for element {element}") - fu_mass += number * element.mass - material.number_density = Value( - material.mass_density.as_unit("g/cm^3") / fu_mass / u2g * 1e21, unit="1/nm^3" - ) + self.ensure_densities() d_head = Value(1.0 / (self._materials[0].number_density.as_unit("1/nm^3") * self.apm.as_unit("nm^2")), "nm") d_tail = Value(1.0 / (self._materials[1].number_density.as_unit("1/nm^3") * self.apm.as_unit("nm^2")), "nm") @@ -322,3 +327,18 @@ def resolve_to_layers(self) -> List[Layer]: tail_2 = Layer(thickness=d_tail, material=m_tail_2, roughness=self.roughness) head_2 = Layer(thickness=d_head, material=m_head_2, roughness=self.roughness) return [head, tail, tail_2, head_2] + + def ensure_densities(self): + for material in [self._environment] + self._materials: + material.generate_density() + if material.number_density is None: + # need to generate number density from mass density and formula + formula = Formula(material.formula, strict=True) + fu_mass = 0.0 + for element, number in formula.elements: + if element.mass is None: + raise ValueError(f"No mass known for element {element}") + fu_mass += number * element.mass + material.number_density = Value( + material.mass_density.as_unit("g/cm^3") / fu_mass / u2g * 1e21, unit="1/nm^3" + ) From ed37922d333ae92d005e82cc16883754d1347105 Mon Sep 17 00:00:00 2001 From: Artur Glavic Date: Wed, 3 Jun 2026 11:25:32 +0200 Subject: [PATCH 11/22] Fix model 10 test values --- examples/sample_model_example_10.yml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/sample_model_example_10.yml b/examples/sample_model_example_10.yml index 8a08e0e..c3c57de 100644 --- a/examples/sample_model_example_10.yml +++ b/examples/sample_model_example_10.yml @@ -42,6 +42,6 @@ sample: # below is just used for orsopy testing and not part of this model test: - stacks: 5 - blocks: 5 - layers: 14 + stacks: 6 + blocks: 6 + layers: 16 From b0584c3db816a4857822efb3c18eca993f329495 Mon Sep 17 00:00:00 2001 From: Artur Glavic Date: Wed, 3 Jun 2026 13:28:55 +0200 Subject: [PATCH 12/22] Calculate number density in Material and don't overwrite Header __repr__ with dataclass default --- examples/sample_model_example_10.yml | 2 +- orsopy/fileio/base.py | 31 +++++++------ orsopy/fileio/data_source.py | 10 ++--- orsopy/fileio/model_building_blocks.py | 62 ++++++++++++++++++++------ orsopy/fileio/model_complex.py | 32 +++---------- orsopy/fileio/model_language.py | 6 +-- orsopy/fileio/orso.py | 4 +- orsopy/fileio/reduction.py | 4 +- orsopy/fileio/tests/test_base.py | 28 ++++++------ orsopy/fileio/tests/test_orso.py | 2 +- 10 files changed, 96 insertions(+), 85 deletions(-) diff --git a/examples/sample_model_example_10.yml b/examples/sample_model_example_10.yml index c3c57de..b41ab64 100644 --- a/examples/sample_model_example_10.yml +++ b/examples/sample_model_example_10.yml @@ -36,7 +36,7 @@ sample: number_density: 0.3135 DM: formula: C26H54 - number_density: 0.128 + mass_density: 0.0779 globals: default_solvent: D2O diff --git a/orsopy/fileio/base.py b/orsopy/fileio/base.py index 5b2ac12..6fee576 100644 --- a/orsopy/fileio/base.py +++ b/orsopy/fileio/base.py @@ -39,9 +39,9 @@ def _noop(self, *args, **kw): yaml.emitter.Emitter.process_tag = _noop # make sure that datetime strings get loaded as str not datetime instances -yaml.constructor.SafeConstructor.yaml_constructors["tag:yaml.org,2002:timestamp"] = ( - yaml.constructor.SafeConstructor.yaml_constructors["tag:yaml.org,2002:str"] -) +yaml.constructor.SafeConstructor.yaml_constructors[ + "tag:yaml.org,2002:timestamp" +] = yaml.constructor.SafeConstructor.yaml_constructors["tag:yaml.org,2002:str"] class ORSOResolveError(ValueError): @@ -296,8 +296,7 @@ def _resolve_type(hint: type, item: Any) -> Any: return item else: warnings.warn( - f"Has to be one of {get_args(hint)} got {item}", - ORSOSchemaWarning, + f"Has to be one of {get_args(hint)} got {item}", ORSOSchemaWarning, ) return str(item) return None @@ -562,7 +561,7 @@ def get_unit_registry(): return unit_registry -@dataclass +@dataclass(repr=False) class ErrorValue(Header): """ Information about errors on a value. @@ -612,7 +611,7 @@ def sigma(self): return self.error_value -@dataclass +@dataclass(repr=False) class Value(Header): """ A value or list of values with an optional unit. @@ -652,7 +651,7 @@ def as_unit(self, output_unit): return val.to(output_unit).magnitude -@dataclass +@dataclass(repr=False) class ComplexValue(Header): """ A value or list of values with an optional unit. @@ -692,7 +691,7 @@ def as_unit(self, output_unit): return val.to(output_unit).magnitude -@dataclass +@dataclass(repr=False) class ValueRange(Header): """ A range or list of ranges with mins, maxs, and an optional unit. @@ -731,7 +730,7 @@ def as_unit(self, output_unit): return (vmin.to(output_unit).magnitude, vmax.to(output_unit).magnitude) -@dataclass +@dataclass(repr=False) class ValueVector(Header): """ A vector or list of vectors with an optional unit. @@ -769,7 +768,7 @@ def as_unit(self, output_unit): return (vx.to(output_unit).magnitude, vy.to(output_unit).magnitude, vz.to(output_unit).magnitude) -@dataclass +@dataclass(repr=False) class AlternatingField(Header): """ A physical field with regular variations as AC magnetic field. @@ -780,7 +779,7 @@ class AlternatingField(Header): phase: Optional[Value] = None -@dataclass +@dataclass(repr=False) class Person(Header): """ Information about a person, including name, affiliation(s), and contact @@ -792,7 +791,7 @@ class Person(Header): contact: Optional[str] = field(default=None, metadata={"description": "Contact (email) address"}) -@dataclass +@dataclass(repr=False) class Column(Header): """ Information about a data column. @@ -817,7 +816,7 @@ class Column(Header): yaml_representer = Header.yaml_representer_compact -@dataclass +@dataclass(repr=False) class ErrorColumn(Header): """ Information about a data column. @@ -872,7 +871,7 @@ def to_sigma(self): return 1.0 -@dataclass +@dataclass(repr=False) class ContentHash(Header): """ A hash of some content, using standard algorithms @@ -922,7 +921,7 @@ def check_file(self, fname: Union[TextIO, str]): return self.digest == digest.hexdigest() -@dataclass +@dataclass(repr=False) class File(Header): """ A file with file path and a last modified timestamp. diff --git a/orsopy/fileio/data_source.py b/orsopy/fileio/data_source.py index bb9b2b3..8eec5bb 100644 --- a/orsopy/fileio/data_source.py +++ b/orsopy/fileio/data_source.py @@ -18,7 +18,7 @@ from .typing_backport import Literal -@dataclass +@dataclass(repr=False) class Experiment(Header): """ A definition of the experiment performed. @@ -43,7 +43,7 @@ class Experiment(Header): doi: Optional[str] = None -@dataclass +@dataclass(repr=False) class Sample(Header): """ A description of the sample measured. @@ -128,7 +128,7 @@ def yaml_representer(self, dumper: yaml.Dumper): return dumper.represent_str(output) -@dataclass +@dataclass(repr=False) class InstrumentSettings(Header): """ Settings associated with the instrumentation. @@ -158,7 +158,7 @@ class InstrumentSettings(Header): __repr__ = Header._staggered_repr -@dataclass +@dataclass(repr=False) class Measurement(Header): """ The measurement elements for the header. @@ -178,7 +178,7 @@ class Measurement(Header): __repr__ = Header._staggered_repr -@dataclass +@dataclass(repr=False) class DataSource(Header): """ The data_source object definition. diff --git a/orsopy/fileio/model_building_blocks.py b/orsopy/fileio/model_building_blocks.py index c6f7634..9976c5a 100644 --- a/orsopy/fileio/model_building_blocks.py +++ b/orsopy/fileio/model_building_blocks.py @@ -15,22 +15,26 @@ class SubStackType(ABC): @property @abstractmethod - def sub_stack_class(self) -> str: ... + def sub_stack_class(self) -> str: + ... @abstractmethod - def resolve_names(self, resolvable_items): ... + def resolve_names(self, resolvable_items): + ... @abstractmethod - def resolve_defaults(self, defaults: "ModelParameters"): ... + def resolve_defaults(self, defaults: "ModelParameters"): + ... @abstractmethod - def resolve_to_layers(self) -> List["Layer"]: ... + def resolve_to_layers(self) -> List["Layer"]: + ... def resolve_to_blocks(self) -> List[Union["Layer", "SubStackType"]]: return [self] -@dataclass +@dataclass(repr=False) class Material(Header): formula: Optional[str] = None mass_density: Optional[Union[float, Value]] = None @@ -68,9 +72,21 @@ def resolve_defaults(self, defaults: "ModelParameters"): elif not isinstance(self.magnetic_moment, Value): self.magnetic_moment = Value(self.magnetic_moment, unit=defaults.magnetic_moment_unit) + def _number_from_mass_density(self): + """ + Use chamical formula and mass density to define number density. + """ + mass_density = self.mass_density.as_unit("g/cm^3") + from ..slddb.material import Material + + material = Material(self.formula, dens=mass_density) + self.number_density = Value(magnitude=material.fu_dens * 1000.0, unit="1/nm^3") + def generate_density(self): if self.sld is not None or self.mass_density is not None or self.number_density is not None: # this material already contains density information + if self.number_density is None and self.mass_density is not None: + self._number_from_mass_density() return if len(DENSITY_RESOLVERS) == 0: from ..utils.resolver_slddb import ResolverSLDDB @@ -160,7 +176,7 @@ def get_sld(self, xray_energy=None) -> complex: return 0.0j -@dataclass +@dataclass(repr=False) class Composit(Header): composition: Dict[str, float] @@ -190,25 +206,43 @@ def resolve_defaults(self, defaults: "ModelParameters"): def generate_density(self, xray_energy=None): """ Create a material based on the composition attribute. + + If all items define a formula, return radiation agnostic Material. """ sld = 0.0 + use_formula = True for key, value in self.composition.items(): mi = self._composition_materials[key] mi.generate_density() - sldi = mi.get_sld(xray_energy=xray_energy) - sld += value * sldi + use_formula &= mi.formula is not None + mix_str = ";".join([f"{value}x{key}" for key, value in self.composition.items()]) - return Material( - sld=ComplexValue(real=sld.real, imag=sld.imag, unit="1/angstrom^2"), - comment=f"composition material: {mix_str}", - ) + + if use_formula: + reference_density = self._composition_materials[list(self.composition.keys())[0]].number_density + rd_unit = reference_density.unit + rd = reference_density.magnitude + FU = "" + for key, value in self.composition.items(): + mi = self._composition_materials[key] + FU += f"({mi.formula}){mi.number_density.as_unit(rd_unit)/rd*value}" + return Material(formula=FU, number_density=reference_density, comment=f"composition material: {mix_str}") + else: + for key, value in self.composition.items(): + mi = self._composition_materials[key] + sldi = mi.get_sld(xray_energy=xray_energy) + sld += value * sldi + return Material( + sld=ComplexValue(real=sld.real, imag=sld.imag, unit="1/angstrom^2"), + comment=f"composition material: {mix_str}", + ) def get_sld(self, xray_energy=None): material = self.generate_density(xray_energy=xray_energy) return material.get_sld(xray_energy=xray_energy) -@dataclass +@dataclass(repr=False) class ModelParameters(Header): roughness: Optional[Value] = field(default_factory=lambda: Value(0.3, "nm")) length_unit: Optional[str] = "nm" @@ -235,7 +269,7 @@ def __post_init__(self): CACHED_MATERIALS = {} -@dataclass +@dataclass(repr=False) class Layer(Header): thickness: Optional[Union[float, Value]] = None roughness: Optional[Union[float, Value]] = None diff --git a/orsopy/fileio/model_complex.py b/orsopy/fileio/model_complex.py index 4ea453c..115a8bc 100644 --- a/orsopy/fileio/model_complex.py +++ b/orsopy/fileio/model_complex.py @@ -14,7 +14,7 @@ from .model_building_blocks import SPECIAL_MATERIALS, Composit, Layer, Material, ModelParameters, SubStackType -@dataclass +@dataclass(repr=False) class FunctionTwoElements(Header, SubStackType): """ Models a continuous variation between two materials/SLDs according to an analytical function. @@ -119,7 +119,7 @@ def mix_hydrate_material(material: Material, solvent: Material, hydtration: floa return Material(formula=FU, number_density=material.number_density) -@dataclass +@dataclass(repr=False) class Leaflet(Header, SubStackType): """ Building block corresponding to a single layer of lipids with heads and tails @@ -216,20 +216,9 @@ def resolve_to_layers(self) -> List[Layer]: def ensure_densities(self): for material in [self._environment] + self._materials: material.generate_density() - if material.number_density is None: - # need to generate number density from mass density and formula - formula = Formula(material.formula, strict=True) - fu_mass = 0.0 - for element, number in formula.elements: - if element.mass is None: - raise ValueError(f"No mass known for element {element}") - fu_mass += number * element.mass - material.number_density = Value( - material.mass_density.as_unit("g/cm^3") / fu_mass / u2g * 1e21, unit="1/nm^3" - ) - - -@dataclass + + +@dataclass(repr=False) class Bilayer(Header, SubStackType): """ Building block corresponding to a bilayer of lipids with outer (heads) and inner (tails) @@ -331,14 +320,3 @@ def resolve_to_layers(self) -> List[Layer]: def ensure_densities(self): for material in [self._environment] + self._materials: material.generate_density() - if material.number_density is None: - # need to generate number density from mass density and formula - formula = Formula(material.formula, strict=True) - fu_mass = 0.0 - for element, number in formula.elements: - if element.mass is None: - raise ValueError(f"No mass known for element {element}") - fu_mass += number * element.mass - material.number_density = Value( - material.mass_density.as_unit("g/cm^3") / fu_mass / u2g * 1e21, unit="1/nm^3" - ) diff --git a/orsopy/fileio/model_language.py b/orsopy/fileio/model_language.py index 4880a88..b1efb29 100644 --- a/orsopy/fileio/model_language.py +++ b/orsopy/fileio/model_language.py @@ -42,7 +42,7 @@ def find_closing(string, start, brackets="()"): return -1 -@dataclass +@dataclass(repr=False) class SubStack(Header, SubStackType): repetitions: int = 1 stack: Optional[str] = None @@ -212,7 +212,7 @@ def resolve_to_layers(self) -> List[Layer]: SUBSTACK_TYPE = Union[SUBSTACK_TYPE, T] -@dataclass +@dataclass(repr=False) class ItemChanger(Header): """ Allows to define a simple change in SubStackType item by @@ -224,7 +224,7 @@ class ItemChanger(Header): original_name = None -@dataclass +@dataclass(repr=False) class SampleModel(Header): stack: str origin: Optional[str] = None diff --git a/orsopy/fileio/orso.py b/orsopy/fileio/orso.py index 0199bd5..55dca55 100644 --- a/orsopy/fileio/orso.py +++ b/orsopy/fileio/orso.py @@ -19,7 +19,7 @@ ) -@dataclass +@dataclass(repr=False) class Orso(Header): """ The Orso object collects the necessary metadata. @@ -140,7 +140,7 @@ def _to_object_dict(self): return out -@dataclass +@dataclass(repr=False) class OrsoDataset: """ :param info: The header information for the reflectivity measurement diff --git a/orsopy/fileio/reduction.py b/orsopy/fileio/reduction.py index 70f43d4..3529326 100644 --- a/orsopy/fileio/reduction.py +++ b/orsopy/fileio/reduction.py @@ -10,7 +10,7 @@ from .base import Header, Person -@dataclass +@dataclass(repr=False) class Software(Header): """ Software description. @@ -27,7 +27,7 @@ class Software(Header): yaml_representer = Header.yaml_representer_compact -@dataclass +@dataclass(repr=False) class Reduction(Header): """ A description of the reduction that has been performed. diff --git a/orsopy/fileio/tests/test_base.py b/orsopy/fileio/tests/test_base.py index e5fd265..2f2e344 100644 --- a/orsopy/fileio/tests/test_base.py +++ b/orsopy/fileio/tests/test_base.py @@ -28,7 +28,7 @@ class TestHeaderClass(unittest.TestCase): """ def test_resolve_any(self): - @dataclass + @dataclass(repr=False) class TestAny(base.Header): test: Any @@ -37,7 +37,7 @@ class TestAny(base.Header): assert res.test is test_object def test_resolve_datetime(self): - @dataclass + @dataclass(repr=False) class TestDatetime(base.Header): test: datetime @@ -55,7 +55,7 @@ class MockDatetime: def strptime(item, format): return datetime.strptime(item, format) - @dataclass + @dataclass(repr=False) class TestDatetime(base.Header): test: MockDatetime @@ -80,7 +80,7 @@ def test_resolve_dictof(self): # dict type annotation changed in 3.8 return - @dataclass + @dataclass(repr=False) class TestDictof(base.Header): test: Dict[str, datetime] test2: Dict[float, int] @@ -92,7 +92,7 @@ class TestDictof(base.Header): assert [(type(ki), type(vi)) for ki, vi in res.test2.items()] == [(float, int), (float, int)] # Dict should be used to define key/value types - @dataclass + @dataclass(repr=False) class TestDictNodef(base.Header): test: Dict @@ -101,7 +101,7 @@ class TestDictNodef(base.Header): assert res.test == {"ab": "cd"} def test_resolve_list(self): - @dataclass + @dataclass(repr=False) class TestList(base.Header): test: List[int] @@ -115,7 +115,7 @@ class TestList(base.Header): assert res.test == [134] def test_resolve_tuple(self): - @dataclass + @dataclass(repr=False) class TestTuple(base.Header): test: Tuple[int, int, int] @@ -134,13 +134,13 @@ class TestTuple(base.Header): assert res.test == (1, 2, 3) def test_subsubclass(self): - @dataclass + @dataclass(repr=False) class Test1(base.Header): pass with self.assertRaises(NotImplementedError): - @dataclass + @dataclass(repr=False) class Test2(Test1): pass @@ -150,7 +150,7 @@ def test_unexpected_error(self): reported as ORSOResolveError. """ - @dataclass + @dataclass(repr=False) class TestDatetime(base.Header): test: datetime @@ -166,7 +166,7 @@ def test_empty_creation(self): Make a class that tests all cases of the empty method. """ - @dataclass + @dataclass(repr=False) class TestEmpty(base.Header): test1: float test2: base.Value @@ -197,7 +197,7 @@ def _ast(self): res = base._todict(Test1()) assert res == {"a": "b"} - @dataclass + @dataclass(repr=False) class Test2: test: int test2: float @@ -207,7 +207,7 @@ class Test2: assert res == {"test": 13, "test2": 12.4, "test3": "1234", "TestClassKey": "Test2"} def test_empty_optional(self): - @dataclass + @dataclass(repr=False) class TestOptional(base.Header): test1: str test2: Optional[float] = None @@ -221,7 +221,7 @@ class TestOptional(base.Header): self.assertEqual(t2.to_yaml().strip(), "test1: abc") def test_default_optional(self): - @dataclass + @dataclass(repr=False) class TestOptional(base.Header): test1: str test2: Optional[float] = 14.5 diff --git a/orsopy/fileio/tests/test_orso.py b/orsopy/fileio/tests/test_orso.py index 29f3331..31b3090 100644 --- a/orsopy/fileio/tests/test_orso.py +++ b/orsopy/fileio/tests/test_orso.py @@ -265,7 +265,7 @@ def test_save_numpy_scalar_dtypes(self): assert i_n.incident_angle.magnitude == 2 def test_nxs_special_cases(self): - @dataclass + @dataclass(repr=False) class TestNxs(Header): test: list test2: tuple From fe1185d43ad3a6701a91e142e633cc42a55adf0b Mon Sep 17 00:00:00 2001 From: Artur Glavic Date: Wed, 3 Jun 2026 17:07:10 +0200 Subject: [PATCH 13/22] Fix mixing calculation, add hydration to thickness calculation, create neutorn refnx model with blocks --- examples/resolve_and_plot_model.py | 112 +++++++++++++- examples/sample_model_example_10.yml | 18 ++- orsopy/fileio/model_building_blocks.py | 9 +- orsopy/fileio/model_complex.py | 198 +++++++++++++------------ 4 files changed, 230 insertions(+), 107 deletions(-) diff --git a/examples/resolve_and_plot_model.py b/examples/resolve_and_plot_model.py index b3ce5db..b7854a2 100644 --- a/examples/resolve_and_plot_model.py +++ b/examples/resolve_and_plot_model.py @@ -16,13 +16,107 @@ from matplotlib import pyplot from numpy import linspace -from refnx.reflect import SLD, ReflectModel, Structure +from refnx.reflect import SLD, LipidLeaflet, ReflectModel, Slab, Structure -from orsopy.fileio import model_language +from orsopy.fileio import model_complex, model_language q = linspace(0.001, 0.2, 200) +def res_layer(layer: model_language.Layer): + m = SLD(layer.material.get_sld() * 1e6) + return Slab( + layer.thickness.as_unit("angstrom"), + m, + layer.roughness.as_unit("angstrom"), + name=getattr(layer, "original_name", ""), + ) + + +def res_leaflet(leaflet: model_complex.Leaflet): + apm = leaflet.apm.as_unit("angstrom^2") + vm_heads = leaflet.heads.volume.as_unit("angstrom^3") + vfrac_head = 1.0 - leaflet.heads_hydration + vfrac_tail = 1.0 - leaflet.tails_hydration + b_heads = leaflet.heads.get_sld() * 1e6 / vm_heads + d_heads = vm_heads / apm / vfrac_head + vm_tails = leaflet.tails.volume.as_unit("angstrom^3") + b_tails = leaflet.tails.get_sld() * 1e6 / vm_tails + d_tails = vm_tails / apm / vfrac_tail + solvent = leaflet.solvent.get_sld() * 1e6 + ll = LipidLeaflet( + apm=apm, + b_heads=b_heads, + vm_heads=vm_heads, + thickness_heads=d_heads, + b_tails=b_tails, + vm_tails=vm_tails, + thickness_tails=d_tails, + rough_head_tail=leaflet.roughness.as_unit("angstrom"), + rough_preceding_mono=leaflet.roughness.as_unit("angstrom"), + reverse_monolayer=not leaflet.heads_first, + name="LL " + getattr(leaflet, "original_name", ""), + head_solvent=solvent, + tail_solvent=solvent, + ) + return ll + + +def res_bilayer(bilayer: model_complex.Bilayer): + apm = bilayer.apm.as_unit("angstrom^2") + vm_heads = bilayer.outer.volume.as_unit("angstrom^3") + vfrac_head = 1.0 - bilayer.outer_hydration + vfrac_tail = 1.0 - bilayer.inner_hydration + b_heads = bilayer.outer.get_sld() * 1e6 / vm_heads + d_heads = vm_heads / apm / vfrac_head + vm_tails = bilayer.inner.volume.as_unit("angstrom^3") + b_tails = bilayer.inner.get_sld() * 1e6 / vm_tails + d_tails = vm_tails / apm / vfrac_tail + solvent = bilayer.solvent.get_sld() * 1e6 + ll = LipidLeaflet( + apm=apm, + b_heads=b_heads, + vm_heads=vm_heads, + thickness_heads=d_heads, + b_tails=b_tails, + vm_tails=vm_tails, + thickness_tails=d_tails, + rough_head_tail=bilayer.roughness.as_unit("angstrom"), + rough_preceding_mono=bilayer.roughness.as_unit("angstrom"), + reverse_monolayer=False, + name="LL " + getattr(bilayer, "original_name", ""), + head_solvent=solvent, + tail_solvent=solvent, + ) + vfrac_head = 1.0 - (bilayer.outer_hydration_2 or bilayer.outer_hydration) + vfrac_tail = 1.0 - (bilayer.inner_hydration_2 or bilayer.inner_hydration) + d_heads = vm_heads / apm / vfrac_head + d_tails = vm_tails / apm / vfrac_tail + rll = LipidLeaflet( + apm=apm, + b_heads=b_heads, + vm_heads=vm_heads, + thickness_heads=d_heads, + b_tails=b_tails, + vm_tails=vm_tails, + thickness_tails=d_tails, + rough_head_tail=bilayer.roughness.as_unit("angstrom"), + rough_preceding_mono=bilayer.roughness.as_unit("angstrom"), + reverse_monolayer=True, + name="rLL " + getattr(bilayer, "original_name", ""), + head_solvent=solvent, + tail_solvent=solvent, + ) + return ll | rll + + +refnx_resolvers = { + model_language.Layer: res_layer, + model_complex.Leaflet: res_leaflet, + model_complex.Bilayer: res_bilayer, +} + + def main(txt=None): if txt is None: txt = sys.argv[1] @@ -39,15 +133,21 @@ def main(txt=None): # initial model before resolving any names print(repr(sample), "\n") print("\n".join([repr(ss) for ss in sample.resolve_stack()]), "\n") - print("\n".join([repr(ss) for ss in sample.resolve_to_blocks()]), "\n") + + blocks = sample.resolve_to_blocks() + print("\n".join([repr(ss) for ss in blocks]), "\n") layers = sample.resolve_to_layers() print("\n".join([repr(li) for li in layers])) structure = Structure() - for lj in layers: - m = SLD(lj.material.get_sld() * 1e6) - structure |= m(lj.thickness.as_unit("angstrom"), lj.roughness.as_unit("angstrom")) + for block in blocks: + if type(block) in refnx_resolvers: + structure |= refnx_resolvers[type(block)](block) + else: + for li in block.resolve_to_layers(): + structure |= res_layer(li) + print("\n", structure, "\n") model = ReflectModel(structure, bkg=0.0) structurex = Structure() for lj in layers: diff --git a/examples/sample_model_example_10.yml b/examples/sample_model_example_10.yml index b41ab64..6c0a3eb 100644 --- a/examples/sample_model_example_10.yml +++ b/examples/sample_model_example_10.yml @@ -5,12 +5,18 @@ sample: DMPC: sub_stack_class: Leaflet heads: PC - tails: DM + tails: + composition: + DM: 0.9 + PC: 0.1 heads_first: false biDMPC: sub_stack_class: Bilayer outer: PC - inner: DM + inner: + composition: + DM: 0.9 + PC: 0.1 apm: {magnitude: 0.8, unit: nm^2} outer_hydration: 0.3 inner_hydration: 0.3 @@ -20,10 +26,10 @@ sample: sub_stack_class: Bilayer outer: formula: C10H18O8NP - number_density: 0.3135 # 1/nm^2 + number_density: 3.135 # 1/nm^3 inner: formula: C26H54 - number_density: 0.128 # 1/nm^2 + number_density: 1.28 # 1/nm^3 apm: 0.9 # nm^2 outer_hydration: 0.3 inner_hydration: 0.3 @@ -33,10 +39,10 @@ sample: materials: PC: formula: C10H18O8NP - number_density: 0.3135 + number_density: 3.135 DM: formula: C26H54 - mass_density: 0.0779 + mass_density: 0.779 globals: default_solvent: D2O diff --git a/orsopy/fileio/model_building_blocks.py b/orsopy/fileio/model_building_blocks.py index 9976c5a..56cff87 100644 --- a/orsopy/fileio/model_building_blocks.py +++ b/orsopy/fileio/model_building_blocks.py @@ -222,11 +222,14 @@ def generate_density(self, xray_energy=None): reference_density = self._composition_materials[list(self.composition.keys())[0]].number_density rd_unit = reference_density.unit rd = reference_density.magnitude - FU = "" + FU = Formula([]) for key, value in self.composition.items(): mi = self._composition_materials[key] - FU += f"({mi.formula}){mi.number_density.as_unit(rd_unit)/rd*value}" - return Material(formula=FU, number_density=reference_density, comment=f"composition material: {mix_str}") + fi = Formula(mi.formula, strict=False) + FU += fi * (mi.number_density.as_unit(rd_unit) / rd * value) + return Material( + formula=str(FU), number_density=reference_density, comment=f"composition material: {mix_str}" + ) else: for key, value in self.composition.items(): mi = self._composition_materials[key] diff --git a/orsopy/fileio/model_complex.py b/orsopy/fileio/model_complex.py index 115a8bc..75fa7e4 100644 --- a/orsopy/fileio/model_complex.py +++ b/orsopy/fileio/model_complex.py @@ -119,8 +119,51 @@ def mix_hydrate_material(material: Material, solvent: Material, hydtration: floa return Material(formula=FU, number_density=material.number_density) +class LipidBase: + """ + Common mathods used for resolving materials and defaults in Lipid based blocks. + """ + + _materials: List[Union[Composit, Material]] + apm: Union[float, Value] + coverage: float + roughness: Optional[Union[float, Value]] + + def resolve_defaults(self, defaults: ModelParameters) -> None: + if len(self._materials) != 3: + self._materials.append(defaults.default_solvent) + if isinstance(self.apm, Value) and self.apm.unit is None: + self.apm.unit = defaults.length_unit + "^2" + elif not isinstance(self.apm, Value): + self.apm = Value(self.apm, unit=defaults.length_unit + "^2") + + for mi in self._materials: + mi.resolve_defaults(defaults) + + if self.roughness is None: + self.roughness = defaults.roughness + elif not isinstance(self.roughness, Value): + self.roughness = Value(self.roughness, unit=defaults.length_unit) + elif self.roughness.unit is None: + self.roughness.unit = defaults.length_unit + + def mixed_material(self, material: Material, solvent: Material, hydtration: float): + return mix_hydrate_material(material, solvent, hydtration, self.coverage) + + def ensure_densities(self): + for i, mi in enumerate(self._materials): + res = mi.generate_density() + if res is not None: + # replace Composit by resulting Material + self._materials[i] = res + + @property + def solvent(self): + return self._materials[2] + + @dataclass(repr=False) -class Leaflet(Header, SubStackType): +class Leaflet(Header, LipidBase, SubStackType): """ Building block corresponding to a single layer of lipids with heads and tails molecule definition. The layer is calculated from molecular volume, area per molecule (apm) and @@ -130,8 +173,8 @@ class Leaflet(Header, SubStackType): The order of heads/tails with respect to the beam side can be flipped using heads_first=False. """ - heads: Union[Material, str] - tails: Union[Material, str] + heads: Union[Composit, Material, str] + tails: Union[Composit, Material, str] heads_first: Optional[bool] = True apm: Optional[Union[float, Value]] = field(default_factory=lambda: Value(0.7, unit="nm^2")) heads_hydration: Optional[float] = 0.3 @@ -140,13 +183,14 @@ class Leaflet(Header, SubStackType): roughness: Optional[Union[float, Value]] = None sub_stack_class: Literal["Leaflet"] = "Leaflet" - _environment = None - def resolve_names(self, resolvable_items): self._materials = [] for i, mi in enumerate((self.heads, self.tails)): if isinstance(mi, Material): material = mi + elif isinstance(mi, Composit): + mi.resolve_names(resolvable_items) + material = mi elif mi in resolvable_items: material = resolvable_items[mi] material.original_name = mi @@ -158,37 +202,15 @@ def resolve_names(self, resolvable_items): material.original_name = mi self._materials.append(material) if "environment" in resolvable_items: - self._environment = resolvable_items["environment"] - if not isinstance(self._environment, Material): - if self._environment in resolvable_items: - self._environment = resolvable_items[self._environment] - elif self._environment in SPECIAL_MATERIALS: - self._environment = SPECIAL_MATERIALS[self._environment] + environment = resolvable_items["environment"] + if not isinstance(self._environment, (Composit, Material)): + if environment in resolvable_items: + environment = resolvable_items[environment] + elif environment in SPECIAL_MATERIALS: + environment = SPECIAL_MATERIALS[environment] elif isinstance(self._environment, str): - self._environment = Material(formula=self._environment) - - def resolve_defaults(self, defaults: ModelParameters) -> None: - if not self._environment: - self._environment = defaults.default_solvent - else: - self._environment.resolve_defaults(defaults) - if isinstance(self.apm, Value) and self.apm.unit is None: - self.apm.unit = defaults.length_unit + "^2" - elif not isinstance(self.apm, Value): - self.apm = Value(self.apm, unit=defaults.length_unit + "^2") - - for mi in self._materials: - mi.resolve_defaults(defaults) - - if self.roughness is None: - self.roughness = defaults.roughness - elif not isinstance(self.roughness, Value): - self.roughness = Value(self.roughness, unit=defaults.length_unit) - elif self.roughness.unit is None: - self.roughness.unit = defaults.length_unit - - def mixed_material(self, material: Material, solvent: Material, hydtration: float): - return mix_hydrate_material(material, solvent, hydtration, self.coverage) + environment = Material(formula=environment) + self._materials.append(environment) def resolve_to_blocks(self) -> List[Union["Layer", "SubStackType"]]: # Make sure the block includes full material data @@ -200,11 +222,18 @@ def resolve_to_blocks(self) -> List[Union["Layer", "SubStackType"]]: def resolve_to_layers(self) -> List[Layer]: self.ensure_densities() - d_head = Value(1.0 / (self._materials[0].number_density.as_unit("1/nm^3") * self.apm.as_unit("nm^2")), "nm") - d_tail = Value(1.0 / (self._materials[1].number_density.as_unit("1/nm^3") * self.apm.as_unit("nm^2")), "nm") + Vf_head = 1.0 - self.heads_hydration + Vf_tail = 1.0 - self.tails_hydration + + d_head = Value( + 1.0 / (Vf_head * self._materials[0].number_density.as_unit("1/nm^3") * self.apm.as_unit("nm^2")), "nm" + ) + d_tail = Value( + 1.0 / (Vf_tail * self._materials[1].number_density.as_unit("1/nm^3") * self.apm.as_unit("nm^2")), "nm" + ) - m_head = self.mixed_material(self._materials[0], self._environment, self.heads_hydration) - m_tail = self.mixed_material(self._materials[1], self._environment, self.tails_hydration) + m_head = self.mixed_material(self._materials[0], self._materials[2], self.heads_hydration) + m_tail = self.mixed_material(self._materials[1], self._materials[2], self.tails_hydration) head = Layer(thickness=d_head, material=m_head, roughness=self.roughness) tail = Layer(thickness=d_tail, material=m_tail, roughness=self.roughness) @@ -213,13 +242,9 @@ def resolve_to_layers(self) -> List[Layer]: else: return [tail, head] - def ensure_densities(self): - for material in [self._environment] + self._materials: - material.generate_density() - @dataclass(repr=False) -class Bilayer(Header, SubStackType): +class Bilayer(Header, LipidBase, SubStackType): """ Building block corresponding to a bilayer of lipids with outer (heads) and inner (tails) molecule definition. The bilayer is calculated from molecular volume, area per molecule (apm) and @@ -228,8 +253,8 @@ class Bilayer(Header, SubStackType): of the ModelParameters. """ - outer: Union[Material, str] - inner: Union[Material, str] + outer: Union[Composit, Material, str] + inner: Union[Composit, Material, str] apm: Optional[Union[float, Value]] = field(default_factory=lambda: Value(0.7, unit="nm^2")) outer_hydration: Optional[float] = 0.3 inner_hydration: Optional[float] = 0.3 @@ -239,13 +264,14 @@ class Bilayer(Header, SubStackType): roughness: Optional[Union[float, Value]] = None sub_stack_class: Literal["Bilayer"] = "Bilayer" - _environment = None - def resolve_names(self, resolvable_items): self._materials = [] for i, mi in enumerate((self.outer, self.inner)): if isinstance(mi, Material): material = mi + elif isinstance(mi, Composit): + mi.resolve_names(resolvable_items) + material = mi elif mi in resolvable_items: material = resolvable_items[mi] material.original_name = mi @@ -257,37 +283,15 @@ def resolve_names(self, resolvable_items): material.original_name = mi self._materials.append(material) if "environment" in resolvable_items: - self._environment = resolvable_items["environment"] - if not isinstance(self._environment, Material): - if self._environment in resolvable_items: - self._environment = resolvable_items[self._environment] - elif self._environment in SPECIAL_MATERIALS: - self._environment = SPECIAL_MATERIALS[self._environment] + environment = resolvable_items["environment"] + if not isinstance(environment, Material): + if environment in resolvable_items: + environment = resolvable_items[self._environment] + elif environment in SPECIAL_MATERIALS: + environment = SPECIAL_MATERIALS[self._environment] elif isinstance(self._environment, str): - self._environment = Material(formula=self._environment) - - def resolve_defaults(self, defaults: ModelParameters) -> None: - if not self._environment: - self._environment = defaults.default_solvent - else: - self._environment.resolve_defaults(defaults) - if isinstance(self.apm, Value) and self.apm.unit is None: - self.apm.unit = defaults.length_unit + "^2" - elif not isinstance(self.apm, Value): - self.apm = Value(self.apm, unit=defaults.length_unit + "^2") - - for mi in self._materials: - mi.resolve_defaults(defaults) - - if self.roughness is None: - self.roughness = defaults.roughness - elif not isinstance(self.roughness, Value): - self.roughness = Value(self.roughness, unit=defaults.length_unit) - elif self.roughness.unit is None: - self.roughness.unit = defaults.length_unit - - def mixed_material(self, material: Material, solvent: Material, hydtration: float): - return mix_hydrate_material(material, solvent, hydtration, self.coverage) + environment = Material(formula=self._environment) + self._materials.append(environment) def resolve_to_blocks(self) -> List[Union["Layer", "SubStackType"]]: self.ensure_densities() @@ -299,24 +303,34 @@ def resolve_to_blocks(self) -> List[Union["Layer", "SubStackType"]]: def resolve_to_layers(self) -> List[Layer]: self.ensure_densities() - d_head = Value(1.0 / (self._materials[0].number_density.as_unit("1/nm^3") * self.apm.as_unit("nm^2")), "nm") - d_tail = Value(1.0 / (self._materials[1].number_density.as_unit("1/nm^3") * self.apm.as_unit("nm^2")), "nm") + Vf_head_1 = 1.0 - self.outer_hydration + Vf_tail_1 = 1.0 - self.inner_hydration + Vf_head_2 = 1.0 - (self.outer_hydration_2 or self.outer_hydration) + Vf_tail_2 = 1.0 - (self.inner_hydration_2 or self.inner_hydration) + d_head_1 = Value( + 1.0 / (Vf_head_1 * self._materials[0].number_density.as_unit("1/nm^3") * self.apm.as_unit("nm^2")), "nm" + ) + d_head_2 = Value( + 1.0 / (Vf_head_2 * self._materials[0].number_density.as_unit("1/nm^3") * self.apm.as_unit("nm^2")), "nm" + ) + d_tail_1 = Value( + 1.0 / (Vf_tail_1 * self._materials[1].number_density.as_unit("1/nm^3") * self.apm.as_unit("nm^2")), "nm" + ) + d_tail_2 = Value( + 1.0 / (Vf_tail_2 * self._materials[1].number_density.as_unit("1/nm^3") * self.apm.as_unit("nm^2")), "nm" + ) - m_head_1 = self.mixed_material(self._materials[0], self._environment, self.outer_hydration) + m_head_1 = self.mixed_material(self._materials[0], self._materials[2], self.outer_hydration) m_head_2 = self.mixed_material( - self._materials[0], self._environment, self.outer_hydration_2 or self.outer_hydration + self._materials[0], self._materials[2], self.outer_hydration_2 or self.outer_hydration ) - m_tail_1 = self.mixed_material(self._materials[1], self._environment, self.inner_hydration) + m_tail_1 = self.mixed_material(self._materials[1], self._materials[2], self.inner_hydration) m_tail_2 = self.mixed_material( - self._materials[1], self._environment, self.inner_hydration_2 or self.inner_hydration + self._materials[1], self._materials[2], self.inner_hydration_2 or self.inner_hydration ) - head = Layer(thickness=d_head, material=m_head_1, roughness=self.roughness) - tail = Layer(thickness=d_tail, material=m_tail_1, roughness=self.roughness) - tail_2 = Layer(thickness=d_tail, material=m_tail_2, roughness=self.roughness) - head_2 = Layer(thickness=d_head, material=m_head_2, roughness=self.roughness) + head = Layer(thickness=d_head_1, material=m_head_1, roughness=self.roughness) + tail = Layer(thickness=d_tail_1, material=m_tail_1, roughness=self.roughness) + tail_2 = Layer(thickness=d_tail_2, material=m_tail_2, roughness=self.roughness) + head_2 = Layer(thickness=d_head_2, material=m_head_2, roughness=self.roughness) return [head, tail, tail_2, head_2] - - def ensure_densities(self): - for material in [self._environment] + self._materials: - material.generate_density() From 4fe6ee6be016580f9cb72e501413945051750639 Mon Sep 17 00:00:00 2001 From: Artur Glavic Date: Thu, 4 Jun 2026 09:14:02 +0200 Subject: [PATCH 14/22] Allow named component resolution with LateResolver, currently implemented for hardcoded Lipids --- examples/sample_model_example_10.yml | 2 +- orsopy/fileio/model_building_blocks.py | 62 +++++++++++++++++++++++++- orsopy/fileio/model_complex.py | 39 ++++++++++++++++ orsopy/fileio/model_language.py | 16 ++++++- 4 files changed, 115 insertions(+), 4 deletions(-) diff --git a/examples/sample_model_example_10.yml b/examples/sample_model_example_10.yml index 6c0a3eb..3f303a8 100644 --- a/examples/sample_model_example_10.yml +++ b/examples/sample_model_example_10.yml @@ -1,6 +1,6 @@ sample: model: - stack: "Si | DMPC | Bilayer{outer: PC, inner: DM} | biDMPC | (lbl3) in H2O | H2O" + stack: "Si | DMPC | Bilayer{outer: PC, inner: DM} | biDMPC | (lbl3 | Bilayer{DMPC}) in H2O | Leaflet{rDMPC} | H2O" sub_stacks: DMPC: sub_stack_class: Leaflet diff --git a/orsopy/fileio/model_building_blocks.py b/orsopy/fileio/model_building_blocks.py index 56cff87..29b9df4 100644 --- a/orsopy/fileio/model_building_blocks.py +++ b/orsopy/fileio/model_building_blocks.py @@ -1,6 +1,6 @@ from abc import ABC, abstractmethod from dataclasses import dataclass, field -from typing import Dict, List, Optional, Union +from typing import Dict, List, Optional, Type, Union from ..utils.chemical_formula import Formula from ..utils.density_resolver import MaterialResolver @@ -9,6 +9,50 @@ DENSITY_RESOLVERS: List[MaterialResolver] = [] +class LateResolver: + """ + Placeholder object for later resolved named SubStackType classes. + Stores the name and potential resolution information for later + evaluation, so no resolution is needed for SampleModel.resolve_stack call. + """ + + name: str + sub_stack_class: Type["SubStackType"] + _resolved_object: Optional["SubStackType"] = None + + def __init__(self, name, sub_stack_class: Type["SubStackType"]): + self.name = name + self.sub_stack_class = sub_stack_class + + def resolve_names(self, resolvable_items): + self.resolvable_items = resolvable_items + + def resolve_defaults(self, defaults: "ModelParameters"): + self.defaults = defaults + + def resolve_object(self): + self._resolved_object = self.sub_stack_class.resolve_name(self.name) + self._resolved_object.original_name = self.name + self._resolved_object.resolve_names(self.resolvable_items) + self._resolved_object.resolve_defaults(self.defaults) + + def resolve_to_blocks(self) -> List[Union["Layer", "SubStackType"]]: + if self._resolved_object is None: + self.resolve_object() + return self._resolved_object.resolve_to_blocks() + + def resolve_to_layers(self) -> List["Layer"]: + if self._resolved_object is None: + self.resolve_object() + return self._resolved_object.resolve_to_layers() + + def __repr__(self): + if self._resolved_object is None: + return "LateResolver('" + self.sub_stack_class.__name__ + "{" + self.name + "}'" + ")" + else: + return "LateResolver(" + repr(self._resolved_object) + ")" + + class SubStackType(ABC): # Protocol for all items that can be placed in sub_stack _orso_name_export_priority = ["sub_stack_class"] @@ -33,6 +77,22 @@ def resolve_to_layers(self) -> List["Layer"]: def resolve_to_blocks(self) -> List[Union["Layer", "SubStackType"]]: return [self] + @classmethod + def resolve_name(cls, name): + """ + Actually perform resolution of object by name. + """ + raise NotImplementedError(f"{cls.__name__} does not implement name based resolution") + + @classmethod + def from_name(cls, name): + """ + This is called if a SubStackType class shall be created from a name at a later stage. + Stores the name and potential parameters in a LateResolver object that is replaced + with the actual instance on resolve_to_blocks or resolve_to_layers calls. + """ + return LateResolver(name, cls) + @dataclass(repr=False) class Material(Header): diff --git a/orsopy/fileio/model_complex.py b/orsopy/fileio/model_complex.py index 75fa7e4..e811dfe 100644 --- a/orsopy/fileio/model_complex.py +++ b/orsopy/fileio/model_complex.py @@ -161,6 +161,13 @@ def ensure_densities(self): def solvent(self): return self._materials[2] + known_lipids = { + "DMPC": { + "heads": Material(formula="C10H18O8NP", number_density=Value(3.135, unit="1/nm^3")), + "tails": Material(formula="C26H54", number_density=Value(1.2, unit="1/nm^3")), + } + } + @dataclass(repr=False) class Leaflet(Header, LipidBase, SubStackType): @@ -183,6 +190,19 @@ class Leaflet(Header, LipidBase, SubStackType): roughness: Optional[Union[float, Value]] = None sub_stack_class: Literal["Leaflet"] = "Leaflet" + @classmethod + def resolve_name(cls, name): + kwds = {} + if name.startswith("r"): + name = name[1:] + kwds["heads_first"] = False + if name in cls.known_lipids: + data = cls.known_lipids[name] + kwds.update(data) + return cls.from_dict(kwds) + else: + raise ValueError(f"Unknown lipid name: {name}") + def resolve_names(self, resolvable_items): self._materials = [] for i, mi in enumerate((self.heads, self.tails)): @@ -264,6 +284,25 @@ class Bilayer(Header, LipidBase, SubStackType): roughness: Optional[Union[float, Value]] = None sub_stack_class: Literal["Bilayer"] = "Bilayer" + @classmethod + def resolve_name(cls, name): + kwds = {} + if name in cls.known_lipids: + data = cls.known_lipids[name] + for src, dest in [ + ("heads", "outer"), + ("tails", "inner"), + ("apm", "apm"), + ("outer_hydration", "heads_hydration"), + ("inner_hydration", "tails_hydration"), + ("roughness", "roughness"), + ]: + if src in data: + kwds[dest] = data[src] + return cls.from_dict(kwds) + else: + raise ValueError(f"Unknown lipid name: {name}") + def resolve_names(self, resolvable_items): self._materials = [] for i, mi in enumerate((self.outer, self.inner)): diff --git a/orsopy/fileio/model_language.py b/orsopy/fileio/model_language.py index b1efb29..ae44ba2 100644 --- a/orsopy/fileio/model_language.py +++ b/orsopy/fileio/model_language.py @@ -99,7 +99,13 @@ def resolve_names(self, resolvable_items): if T.sub_stack_class == class_name: ssclass = T break - obj = ssclass.from_dict(yaml.load(stack[open_idx : close_idx + 1], Loader=yaml.FullLoader)) + class_data = stack[open_idx : close_idx + 1] + if ":" in class_data: + # data describes yaml dictionary of keywords for this class + obj = ssclass.from_dict(yaml.load(class_data, Loader=yaml.FullLoader)) + else: + # data describes the name of an object to be resolved later + obj = ssclass.from_name(class_data[1:-1].strip()) try: thickness = float(stack[close_idx + 1 : next_idx].strip()) except ValueError: @@ -309,7 +315,13 @@ def resolve_stack(self): if T.sub_stack_class == class_name: ssclass = T break - obj = ssclass.from_dict(yaml.load(stack[open_idx : close_idx + 1], Loader=yaml.FullLoader)) + class_data = stack[open_idx : close_idx + 1] + if ":" in class_data: + # data describes yaml dictionary of keywords for this class + obj = ssclass.from_dict(yaml.load(class_data, Loader=yaml.FullLoader)) + else: + # data describes the name of an object to be resolved later + obj = ssclass.from_name(class_data[1:-1].strip()) try: thickness = float(stack[close_idx + 1 : next_idx].strip()) except ValueError: From f5b4655ade0db1759550842df85d88389ac8bcf8 Mon Sep 17 00:00:00 2001 From: Artur Glavic Date: Thu, 4 Jun 2026 09:26:37 +0200 Subject: [PATCH 15/22] Add simple example showing multiple lipids with just using stack string --- examples/sample_model_example_10.yml | 48 ++--------------------- examples/sample_model_example_11.yml | 53 ++++++++++++++++++++++++++ orsopy/fileio/model_building_blocks.py | 4 ++ orsopy/fileio/model_complex.py | 34 +++++++++++++++-- 4 files changed, 91 insertions(+), 48 deletions(-) create mode 100644 examples/sample_model_example_11.yml diff --git a/examples/sample_model_example_10.yml b/examples/sample_model_example_10.yml index 3f303a8..124b1a7 100644 --- a/examples/sample_model_example_10.yml +++ b/examples/sample_model_example_10.yml @@ -1,50 +1,8 @@ sample: model: - stack: "Si | DMPC | Bilayer{outer: PC, inner: DM} | biDMPC | (lbl3 | Bilayer{DMPC}) in H2O | Leaflet{rDMPC} | H2O" - sub_stacks: - DMPC: - sub_stack_class: Leaflet - heads: PC - tails: - composition: - DM: 0.9 - PC: 0.1 - heads_first: false - biDMPC: - sub_stack_class: Bilayer - outer: PC - inner: - composition: - DM: 0.9 - PC: 0.1 - apm: {magnitude: 0.8, unit: nm^2} - outer_hydration: 0.3 - inner_hydration: 0.3 - inner_hydration_2: 0.4 - outer_hydration_2: 0.5 - lbl3: - sub_stack_class: Bilayer - outer: - formula: C10H18O8NP - number_density: 3.135 # 1/nm^3 - inner: - formula: C26H54 - number_density: 1.28 # 1/nm^3 - apm: 0.9 # nm^2 - outer_hydration: 0.3 - inner_hydration: 0.3 - inner_hydration_2: 0.4 - outer_hydration_2: 0.5 - coverage: 0.75 - materials: - PC: - formula: C10H18O8NP - number_density: 3.135 - DM: - formula: C26H54 - mass_density: 0.779 - globals: - default_solvent: D2O + stack: > + Si | Leaflet{DMPC} | Bilayer{d54-DMPC} | Bilayer{DMPG} | Bilayer{d54-DMPG} | Bilayer{POPC} + | Bilayer{POPG} | Bilayer{POPS} | Bilayer{POPE} | Leaflet{rDMPC} | H2O # below is just used for orsopy testing and not part of this model test: diff --git a/examples/sample_model_example_11.yml b/examples/sample_model_example_11.yml new file mode 100644 index 0000000..3f303a8 --- /dev/null +++ b/examples/sample_model_example_11.yml @@ -0,0 +1,53 @@ +sample: + model: + stack: "Si | DMPC | Bilayer{outer: PC, inner: DM} | biDMPC | (lbl3 | Bilayer{DMPC}) in H2O | Leaflet{rDMPC} | H2O" + sub_stacks: + DMPC: + sub_stack_class: Leaflet + heads: PC + tails: + composition: + DM: 0.9 + PC: 0.1 + heads_first: false + biDMPC: + sub_stack_class: Bilayer + outer: PC + inner: + composition: + DM: 0.9 + PC: 0.1 + apm: {magnitude: 0.8, unit: nm^2} + outer_hydration: 0.3 + inner_hydration: 0.3 + inner_hydration_2: 0.4 + outer_hydration_2: 0.5 + lbl3: + sub_stack_class: Bilayer + outer: + formula: C10H18O8NP + number_density: 3.135 # 1/nm^3 + inner: + formula: C26H54 + number_density: 1.28 # 1/nm^3 + apm: 0.9 # nm^2 + outer_hydration: 0.3 + inner_hydration: 0.3 + inner_hydration_2: 0.4 + outer_hydration_2: 0.5 + coverage: 0.75 + materials: + PC: + formula: C10H18O8NP + number_density: 3.135 + DM: + formula: C26H54 + mass_density: 0.779 + globals: + default_solvent: D2O + + # below is just used for orsopy testing and not part of this model + test: + stacks: 6 + blocks: 6 + layers: 16 diff --git a/orsopy/fileio/model_building_blocks.py b/orsopy/fileio/model_building_blocks.py index 29b9df4..ce9d9a6 100644 --- a/orsopy/fileio/model_building_blocks.py +++ b/orsopy/fileio/model_building_blocks.py @@ -23,6 +23,8 @@ class LateResolver: def __init__(self, name, sub_stack_class: Type["SubStackType"]): self.name = name self.sub_stack_class = sub_stack_class + if hasattr(sub_stack_class, "thickness"): + self.thickness = None def resolve_names(self, resolvable_items): self.resolvable_items = resolvable_items @@ -33,6 +35,8 @@ def resolve_defaults(self, defaults: "ModelParameters"): def resolve_object(self): self._resolved_object = self.sub_stack_class.resolve_name(self.name) self._resolved_object.original_name = self.name + if getattr(self._resolved_object, "thickness", "ignore") is None: + self._resolved_object.thickness = self.thickness self._resolved_object.resolve_names(self.resolvable_items) self._resolved_object.resolve_defaults(self.defaults) diff --git a/orsopy/fileio/model_complex.py b/orsopy/fileio/model_complex.py index e811dfe..2e7c808 100644 --- a/orsopy/fileio/model_complex.py +++ b/orsopy/fileio/model_complex.py @@ -163,9 +163,37 @@ def solvent(self): known_lipids = { "DMPC": { - "heads": Material(formula="C10H18O8NP", number_density=Value(3.135, unit="1/nm^3")), - "tails": Material(formula="C26H54", number_density=Value(1.2, unit="1/nm^3")), - } + "heads": Material(formula="C10H18O8NP", number_density=Value(1.0 / 0.319, unit="1/nm^3")), + "tails": Material(formula="C26H54", number_density=Value(1.0 / 0.782, unit="1/nm^3")), + }, + "d54-DMPC": { + "heads": Material(formula="C10H18O8NP", number_density=Value(1.0 / 0.319, unit="1/nm^3")), + "tails": Material(formula="C26D54", number_density=Value(1.0 / 0.782, unit="1/nm^3")), + }, + "DMPG": { + "heads": Material(formula="C8H12O10P1", number_density=Value(1.0 / 0.257, unit="1/nm^3")), + "tails": Material(formula="C26H54", number_density=Value(1.0 / 0.782, unit="1/nm^3")), + }, + "d54-DMPG": { + "heads": Material(formula="C8H12O10P1", number_density=Value(1.0 / 0.257, unit="1/nm^3")), + "tails": Material(formula="C26D54", number_density=Value(1.0 / 0.782, unit="1/nm^3")), + }, + "POPC": { + "heads": Material(formula="C10H18O8NP", number_density=Value(1.0 / 0.319, unit="1/nm^3")), + "tails": Material(formula="C32H64", number_density=Value(1.0 / 0.944, unit="1/nm^3")), + }, + "POPG": { + "heads": Material(formula="C8H12O10P", number_density=Value(1.0 / 0.257, unit="1/nm^3")), + "tails": Material(formula="C32H64", number_density=Value(1.0 / 0.944, unit="1/nm^3")), + }, + "POPS": { + "heads": Material(formula="C8H11O10NP", number_density=Value(1.0 / 0.244, unit="1/nm^3")), + "tails": Material(formula="C32H64", number_density=Value(1.0 / 0.944, unit="1/nm^3")), + }, + "POPE": { + "heads": Material(formula="C7H12O8PN", number_density=Value(1.0 / 0.252, unit="1/nm^3")), + "tails": Material(formula="C32H64", number_density=Value(1.0 / 0.944, unit="1/nm^3")), + }, } From 5d885854185f4b72df41035815af684140116054 Mon Sep 17 00:00:00 2001 From: Artur Glavic Date: Thu, 4 Jun 2026 09:31:06 +0200 Subject: [PATCH 16/22] Fix test values for changed models --- examples/sample_model_example_10.yml | 6 +++--- examples/sample_model_example_11.yml | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/examples/sample_model_example_10.yml b/examples/sample_model_example_10.yml index 124b1a7..0986e0b 100644 --- a/examples/sample_model_example_10.yml +++ b/examples/sample_model_example_10.yml @@ -6,6 +6,6 @@ sample: # below is just used for orsopy testing and not part of this model test: - stacks: 6 - blocks: 6 - layers: 16 + stacks: 11 + blocks: 11 + layers: 34 diff --git a/examples/sample_model_example_11.yml b/examples/sample_model_example_11.yml index 3f303a8..84b66f7 100644 --- a/examples/sample_model_example_11.yml +++ b/examples/sample_model_example_11.yml @@ -48,6 +48,6 @@ sample: # below is just used for orsopy testing and not part of this model test: - stacks: 6 - blocks: 6 - layers: 16 + stacks: 7 + blocks: 8 + layers: 22 From faab03dc4495e4c28496e0a0817a83fb9d153f0f Mon Sep 17 00:00:00 2001 From: Artur Glavic Date: Thu, 4 Jun 2026 11:48:26 +0200 Subject: [PATCH 17/22] Still needs to validate refnx implementation of leaflet --- examples/resolve_and_plot_model.py | 44 ++++++++++++++++++++---------- 1 file changed, 29 insertions(+), 15 deletions(-) diff --git a/examples/resolve_and_plot_model.py b/examples/resolve_and_plot_model.py index b7854a2..e35aac1 100644 --- a/examples/resolve_and_plot_model.py +++ b/examples/resolve_and_plot_model.py @@ -23,8 +23,8 @@ q = linspace(0.001, 0.2, 200) -def res_layer(layer: model_language.Layer): - m = SLD(layer.material.get_sld() * 1e6) +def res_layer(layer: model_language.Layer, xray_energy=None): + m = SLD(layer.material.get_sld(xray_energy=xray_energy) * 1e6) return Slab( layer.thickness.as_unit("angstrom"), m, @@ -33,17 +33,17 @@ def res_layer(layer: model_language.Layer): ) -def res_leaflet(leaflet: model_complex.Leaflet): +def res_leaflet(leaflet: model_complex.Leaflet, xray_energy=None): apm = leaflet.apm.as_unit("angstrom^2") vm_heads = leaflet.heads.volume.as_unit("angstrom^3") vfrac_head = 1.0 - leaflet.heads_hydration vfrac_tail = 1.0 - leaflet.tails_hydration - b_heads = leaflet.heads.get_sld() * 1e6 / vm_heads + b_heads = leaflet.heads.get_sld(xray_energy=xray_energy) * 1e5 / vm_heads d_heads = vm_heads / apm / vfrac_head vm_tails = leaflet.tails.volume.as_unit("angstrom^3") - b_tails = leaflet.tails.get_sld() * 1e6 / vm_tails + b_tails = leaflet.tails.get_sld(xray_energy=xray_energy) * 1e5 / vm_tails d_tails = vm_tails / apm / vfrac_tail - solvent = leaflet.solvent.get_sld() * 1e6 + solvent = leaflet.solvent.get_sld(xray_energy=xray_energy) * 1e6 ll = LipidLeaflet( apm=apm, b_heads=b_heads, @@ -62,17 +62,17 @@ def res_leaflet(leaflet: model_complex.Leaflet): return ll -def res_bilayer(bilayer: model_complex.Bilayer): +def res_bilayer(bilayer: model_complex.Bilayer, xray_energy=None): apm = bilayer.apm.as_unit("angstrom^2") vm_heads = bilayer.outer.volume.as_unit("angstrom^3") vfrac_head = 1.0 - bilayer.outer_hydration vfrac_tail = 1.0 - bilayer.inner_hydration - b_heads = bilayer.outer.get_sld() * 1e6 / vm_heads + b_heads = bilayer.outer.get_sld(xray_energy=xray_energy) * 1e5 / vm_heads d_heads = vm_heads / apm / vfrac_head vm_tails = bilayer.inner.volume.as_unit("angstrom^3") - b_tails = bilayer.inner.get_sld() * 1e6 / vm_tails + b_tails = bilayer.inner.get_sld(xray_energy=xray_energy) * 1e5 / vm_tails d_tails = vm_tails / apm / vfrac_tail - solvent = bilayer.solvent.get_sld() * 1e6 + solvent = bilayer.solvent.get_sld(xray_energy=xray_energy) * 1e6 ll = LipidLeaflet( apm=apm, b_heads=b_heads, @@ -140,6 +140,7 @@ def main(txt=None): layers = sample.resolve_to_layers() print("\n".join([repr(li) for li in layers])) + # high-level resolution based on blocks structure = Structure() for block in blocks: if type(block) in refnx_resolvers: @@ -150,11 +151,22 @@ def main(txt=None): print("\n", structure, "\n") model = ReflectModel(structure, bkg=0.0) structurex = Structure() - for lj in layers: - m = SLD(lj.material.get_sld(xray_energy="Cu") * 1e6) - structurex |= m(lj.thickness.as_unit("angstrom"), lj.roughness.as_unit("angstrom")) + for block in blocks: + if type(block) in refnx_resolvers: + structurex |= refnx_resolvers[type(block)](block, xray_energy="Cu") + else: + for li in block.resolve_to_layers(): + structurex |= res_layer(li, xray_energy="Cu") modelx = ReflectModel(structurex, bkg=0.0) + # ORSOpy slab resolution + orsopy_neutron = Structure() + for li in layers: + orsopy_neutron |= res_layer(li) + orsopy_xray = Structure() + for li in layers: + orsopy_xray |= res_layer(li, xray_energy="Cu") + pyplot.figure(figsize=(12, 5)) pyplot.subplot(121) pyplot.semilogy(q, model(q), label="neutron") @@ -164,8 +176,10 @@ def main(txt=None): pyplot.xlabel("q [Å$^{-1}$]") pyplot.ylabel("Neutron-reflectivity") pyplot.subplot(122) - pyplot.plot(*structure.sld_profile(), label="neutron") - pyplot.plot(*structurex.sld_profile(), label="x-ray (Cu)") + pyplot.plot(*structure.sld_profile(), label="neutron", color="C0") + pyplot.plot(*orsopy_neutron.sld_profile(), "--", label="", color="C0") + pyplot.plot(*structurex.sld_profile(), label="x-ray (Cu)", color="C1") + pyplot.plot(*orsopy_xray.sld_profile(), "--", label="", color="C1") pyplot.legend() pyplot.title(txt) pyplot.ylabel("SLD / $10^{-6} \\AA^{-2}$") From 955424eac175a0aead7865cb73dad7c00afa7588 Mon Sep 17 00:00:00 2001 From: Artur Glavic Date: Fri, 5 Jun 2026 16:37:11 +0200 Subject: [PATCH 18/22] Fix resolver for refnx LipidLeaflet to calculate b correctly and include coverage parameter --- examples/resolve_and_plot_model.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/examples/resolve_and_plot_model.py b/examples/resolve_and_plot_model.py index e35aac1..ca2b690 100644 --- a/examples/resolve_and_plot_model.py +++ b/examples/resolve_and_plot_model.py @@ -38,14 +38,14 @@ def res_leaflet(leaflet: model_complex.Leaflet, xray_energy=None): vm_heads = leaflet.heads.volume.as_unit("angstrom^3") vfrac_head = 1.0 - leaflet.heads_hydration vfrac_tail = 1.0 - leaflet.tails_hydration - b_heads = leaflet.heads.get_sld(xray_energy=xray_energy) * 1e5 / vm_heads + b_heads = leaflet.heads.get_sld(xray_energy=xray_energy) * vm_heads d_heads = vm_heads / apm / vfrac_head vm_tails = leaflet.tails.volume.as_unit("angstrom^3") - b_tails = leaflet.tails.get_sld(xray_energy=xray_energy) * 1e5 / vm_tails + b_tails = leaflet.tails.get_sld(xray_energy=xray_energy) * vm_tails d_tails = vm_tails / apm / vfrac_tail solvent = leaflet.solvent.get_sld(xray_energy=xray_energy) * 1e6 ll = LipidLeaflet( - apm=apm, + apm=apm / leaflet.coverage, b_heads=b_heads, vm_heads=vm_heads, thickness_heads=d_heads, @@ -67,14 +67,14 @@ def res_bilayer(bilayer: model_complex.Bilayer, xray_energy=None): vm_heads = bilayer.outer.volume.as_unit("angstrom^3") vfrac_head = 1.0 - bilayer.outer_hydration vfrac_tail = 1.0 - bilayer.inner_hydration - b_heads = bilayer.outer.get_sld(xray_energy=xray_energy) * 1e5 / vm_heads + b_heads = bilayer.outer.get_sld(xray_energy=xray_energy) * vm_heads d_heads = vm_heads / apm / vfrac_head vm_tails = bilayer.inner.volume.as_unit("angstrom^3") - b_tails = bilayer.inner.get_sld(xray_energy=xray_energy) * 1e5 / vm_tails + b_tails = bilayer.inner.get_sld(xray_energy=xray_energy) * vm_tails d_tails = vm_tails / apm / vfrac_tail solvent = bilayer.solvent.get_sld(xray_energy=xray_energy) * 1e6 ll = LipidLeaflet( - apm=apm, + apm=apm / bilayer.coverage, b_heads=b_heads, vm_heads=vm_heads, thickness_heads=d_heads, @@ -93,7 +93,7 @@ def res_bilayer(bilayer: model_complex.Bilayer, xray_energy=None): d_heads = vm_heads / apm / vfrac_head d_tails = vm_tails / apm / vfrac_tail rll = LipidLeaflet( - apm=apm, + apm=apm / bilayer.coverage, b_heads=b_heads, vm_heads=vm_heads, thickness_heads=d_heads, @@ -177,9 +177,9 @@ def main(txt=None): pyplot.ylabel("Neutron-reflectivity") pyplot.subplot(122) pyplot.plot(*structure.sld_profile(), label="neutron", color="C0") - pyplot.plot(*orsopy_neutron.sld_profile(), "--", label="", color="C0") + pyplot.plot(*orsopy_neutron.sld_profile(), "--", lw=2, label="", color="C0") pyplot.plot(*structurex.sld_profile(), label="x-ray (Cu)", color="C1") - pyplot.plot(*orsopy_xray.sld_profile(), "--", label="", color="C1") + pyplot.plot(*orsopy_xray.sld_profile(), "--", lw=2, label="", color="C1") pyplot.legend() pyplot.title(txt) pyplot.ylabel("SLD / $10^{-6} \\AA^{-2}$") From 4be199322f1cebe947cd25bde447db24907743dd Mon Sep 17 00:00:00 2001 From: Artur Glavic Date: Fri, 5 Jun 2026 17:16:30 +0200 Subject: [PATCH 19/22] Change resolve_and_plot_model to first generate refnx model code and print the results --- examples/resolve_and_plot_model.py | 261 +++++++++++++++-------------- 1 file changed, 139 insertions(+), 122 deletions(-) diff --git a/examples/resolve_and_plot_model.py b/examples/resolve_and_plot_model.py index ca2b690..d4878de 100644 --- a/examples/resolve_and_plot_model.py +++ b/examples/resolve_and_plot_model.py @@ -1,6 +1,8 @@ """ Python script to showcase the simple model language by resolving a stack string -and plotting SLD and reflectivity for neutrons. +and plotting SLD and reflectivity for neutrons and x-rays. +The refnx model is first generate as a string and then evaluated to be able to export +the model directly. Script requires matplotlib and refnx to run: pip install matplotlib refnx @@ -16,105 +18,121 @@ from matplotlib import pyplot from numpy import linspace -from refnx.reflect import SLD, LipidLeaflet, ReflectModel, Slab, Structure from orsopy.fileio import model_complex, model_language -q = linspace(0.001, 0.2, 200) + +class RefNxResolver: + def __init__(self, sample: model_language.SampleModel): + self.sample = sample + + def get_model(self, xray_energy=None, layers_only=False): + self.xray_energy = xray_energy + + self.model_header() + if layers_only: + self.resolve_to_layers(self.sample) + else: + for block in self.sample.resolve_to_blocks(): + block_resolver = getattr(self, "res" + type(block).__name__, self.resolve_to_layers) + block_resolver(block) + self.model_end() + return self.model + + def model_header(self): + self.model = "from refnx.reflect import SLD, LipidLeaflet, ReflectModel, Slab, Structure\n\n" + self.model += "structure = Structure()\n\n" + + def model_end(self): + self.model += "model = ReflectModel(structure, bkg=0.0)\n" + + def resolve_to_layers(self, block: model_language.SubStackType): + for li in block.resolve_to_layers(): + self.resLayer(li) + + def resLayer(self, layer: model_language.Layer): + self.model += f"m = SLD({layer.material.get_sld(xray_energy=self.xray_energy) * 1e6})\n" + self.model += ( + f"structure |= Slab({layer.thickness.as_unit('angstrom')}, m, " + f"{layer.roughness.as_unit('angstrom')}, " + f"name='{getattr(layer, 'original_name', '')}')\n\n" + ) + + def resLeaflet(self, leaflet: model_complex.Leaflet): + apm = leaflet.apm.as_unit("angstrom^2") + vm_heads = leaflet.heads.volume.as_unit("angstrom^3") + vfrac_head = 1.0 - leaflet.heads_hydration + vfrac_tail = 1.0 - leaflet.tails_hydration + b_heads = leaflet.heads.get_sld(xray_energy=self.xray_energy) * vm_heads + d_heads = vm_heads / apm / vfrac_head + vm_tails = leaflet.tails.volume.as_unit("angstrom^3") + b_tails = leaflet.tails.get_sld(xray_energy=self.xray_energy) * vm_tails + d_tails = vm_tails / apm / vfrac_tail + solvent = leaflet.solvent.get_sld(xray_energy=self.xray_energy) * 1e6 + self.model += f"""structure |= LipidLeaflet( + apm={apm / leaflet.coverage}, + b_heads={b_heads}, + vm_heads={vm_heads}, + thickness_heads={d_heads}, + b_tails={b_tails}, + vm_tails={vm_tails}, + thickness_tails={d_tails}, + rough_head_tail={leaflet.roughness.as_unit("angstrom")}, + rough_preceding_mono={leaflet.roughness.as_unit("angstrom")}, + reverse_monolayer={not leaflet.heads_first}, + name='{'LL ' + getattr(leaflet, 'original_name', '')}', + head_solvent={solvent}, + tail_solvent={solvent}, + )\n\n""" + + def resBilayer(self, bilayer: model_complex.Bilayer): + apm = bilayer.apm.as_unit("angstrom^2") + vm_heads = bilayer.outer.volume.as_unit("angstrom^3") + vfrac_head = 1.0 - bilayer.outer_hydration + vfrac_tail = 1.0 - bilayer.inner_hydration + b_heads = bilayer.outer.get_sld(xray_energy=self.xray_energy) * vm_heads + d_heads = vm_heads / apm / vfrac_head + vm_tails = bilayer.inner.volume.as_unit("angstrom^3") + b_tails = bilayer.inner.get_sld(xray_energy=self.xray_energy) * vm_tails + d_tails = vm_tails / apm / vfrac_tail + solvent = bilayer.solvent.get_sld(xray_energy=self.xray_energy) * 1e6 + self.model += f"""structure |= LipidLeaflet( + apm={apm / bilayer.coverage}, + b_heads={b_heads}, + vm_heads={vm_heads}, + thickness_heads={d_heads}, + b_tails={b_tails}, + vm_tails={vm_tails}, + thickness_tails={d_tails}, + rough_head_tail={bilayer.roughness.as_unit("angstrom")}, + rough_preceding_mono={bilayer.roughness.as_unit("angstrom")}, + reverse_monolayer=False, + name='{'LL ' + getattr(bilayer, 'original_name', '')}', + head_solvent={solvent}, + tail_solvent={solvent}, + )\n""" + vfrac_head = 1.0 - (bilayer.outer_hydration_2 or bilayer.outer_hydration) + vfrac_tail = 1.0 - (bilayer.inner_hydration_2 or bilayer.inner_hydration) + d_heads = vm_heads / apm / vfrac_head + d_tails = vm_tails / apm / vfrac_tail + self.model += f"""structure |= LipidLeaflet( + apm={apm / bilayer.coverage}, + b_heads={b_heads}, + vm_heads={vm_heads}, + thickness_heads={d_heads}, + b_tails={b_tails}, + vm_tails={vm_tails}, + thickness_tails={d_tails}, + rough_head_tail={bilayer.roughness.as_unit("angstrom")}, + rough_preceding_mono={bilayer.roughness.as_unit("angstrom")}, + reverse_monolayer=True, + name='{'fLL ' + getattr(bilayer, 'original_name', '')}', + head_solvent={solvent}, + tail_solvent={solvent}, + )\n\n""" -def res_layer(layer: model_language.Layer, xray_energy=None): - m = SLD(layer.material.get_sld(xray_energy=xray_energy) * 1e6) - return Slab( - layer.thickness.as_unit("angstrom"), - m, - layer.roughness.as_unit("angstrom"), - name=getattr(layer, "original_name", ""), - ) - - -def res_leaflet(leaflet: model_complex.Leaflet, xray_energy=None): - apm = leaflet.apm.as_unit("angstrom^2") - vm_heads = leaflet.heads.volume.as_unit("angstrom^3") - vfrac_head = 1.0 - leaflet.heads_hydration - vfrac_tail = 1.0 - leaflet.tails_hydration - b_heads = leaflet.heads.get_sld(xray_energy=xray_energy) * vm_heads - d_heads = vm_heads / apm / vfrac_head - vm_tails = leaflet.tails.volume.as_unit("angstrom^3") - b_tails = leaflet.tails.get_sld(xray_energy=xray_energy) * vm_tails - d_tails = vm_tails / apm / vfrac_tail - solvent = leaflet.solvent.get_sld(xray_energy=xray_energy) * 1e6 - ll = LipidLeaflet( - apm=apm / leaflet.coverage, - b_heads=b_heads, - vm_heads=vm_heads, - thickness_heads=d_heads, - b_tails=b_tails, - vm_tails=vm_tails, - thickness_tails=d_tails, - rough_head_tail=leaflet.roughness.as_unit("angstrom"), - rough_preceding_mono=leaflet.roughness.as_unit("angstrom"), - reverse_monolayer=not leaflet.heads_first, - name="LL " + getattr(leaflet, "original_name", ""), - head_solvent=solvent, - tail_solvent=solvent, - ) - return ll - - -def res_bilayer(bilayer: model_complex.Bilayer, xray_energy=None): - apm = bilayer.apm.as_unit("angstrom^2") - vm_heads = bilayer.outer.volume.as_unit("angstrom^3") - vfrac_head = 1.0 - bilayer.outer_hydration - vfrac_tail = 1.0 - bilayer.inner_hydration - b_heads = bilayer.outer.get_sld(xray_energy=xray_energy) * vm_heads - d_heads = vm_heads / apm / vfrac_head - vm_tails = bilayer.inner.volume.as_unit("angstrom^3") - b_tails = bilayer.inner.get_sld(xray_energy=xray_energy) * vm_tails - d_tails = vm_tails / apm / vfrac_tail - solvent = bilayer.solvent.get_sld(xray_energy=xray_energy) * 1e6 - ll = LipidLeaflet( - apm=apm / bilayer.coverage, - b_heads=b_heads, - vm_heads=vm_heads, - thickness_heads=d_heads, - b_tails=b_tails, - vm_tails=vm_tails, - thickness_tails=d_tails, - rough_head_tail=bilayer.roughness.as_unit("angstrom"), - rough_preceding_mono=bilayer.roughness.as_unit("angstrom"), - reverse_monolayer=False, - name="LL " + getattr(bilayer, "original_name", ""), - head_solvent=solvent, - tail_solvent=solvent, - ) - vfrac_head = 1.0 - (bilayer.outer_hydration_2 or bilayer.outer_hydration) - vfrac_tail = 1.0 - (bilayer.inner_hydration_2 or bilayer.inner_hydration) - d_heads = vm_heads / apm / vfrac_head - d_tails = vm_tails / apm / vfrac_tail - rll = LipidLeaflet( - apm=apm / bilayer.coverage, - b_heads=b_heads, - vm_heads=vm_heads, - thickness_heads=d_heads, - b_tails=b_tails, - vm_tails=vm_tails, - thickness_tails=d_tails, - rough_head_tail=bilayer.roughness.as_unit("angstrom"), - rough_preceding_mono=bilayer.roughness.as_unit("angstrom"), - reverse_monolayer=True, - name="rLL " + getattr(bilayer, "original_name", ""), - head_solvent=solvent, - tail_solvent=solvent, - ) - return ll | rll - - -refnx_resolvers = { - model_language.Layer: res_layer, - model_complex.Leaflet: res_leaflet, - model_complex.Bilayer: res_bilayer, -} +q = linspace(0.001, 0.2, 200) def main(txt=None): @@ -137,48 +155,47 @@ def main(txt=None): blocks = sample.resolve_to_blocks() print("\n".join([repr(ss) for ss in blocks]), "\n") + resolver = RefNxResolver(sample) + model_n = resolver.get_model() + model_x = resolver.get_model(xray_energy="Cu") + layers = sample.resolve_to_layers() print("\n".join([repr(li) for li in layers])) # high-level resolution based on blocks - structure = Structure() - for block in blocks: - if type(block) in refnx_resolvers: - structure |= refnx_resolvers[type(block)](block) - else: - for li in block.resolve_to_layers(): - structure |= res_layer(li) - print("\n", structure, "\n") - model = ReflectModel(structure, bkg=0.0) - structurex = Structure() - for block in blocks: - if type(block) in refnx_resolvers: - structurex |= refnx_resolvers[type(block)](block, xray_energy="Cu") - else: - for li in block.resolve_to_layers(): - structurex |= res_layer(li, xray_energy="Cu") - modelx = ReflectModel(structurex, bkg=0.0) + resn = {} + exec(model_n, resn) + resx = {} + exec(model_x, resx) + + print("####################### Neutron Model ##################") + print(model_n) + + print("\n\n####################### X-ray Model ##################") + print(model_n) # ORSOpy slab resolution - orsopy_neutron = Structure() - for li in layers: - orsopy_neutron |= res_layer(li) - orsopy_xray = Structure() - for li in layers: - orsopy_xray |= res_layer(li, xray_energy="Cu") + model = resolver.get_model(layers_only=True) + res = {} + exec(model, res) + orsopy_neutron = res["structure"] + model = resolver.get_model(xray_energy="Cu", layers_only=True) + res = {} + exec(model, res) + orsopy_xray = res["structure"] pyplot.figure(figsize=(12, 5)) pyplot.subplot(121) - pyplot.semilogy(q, model(q), label="neutron") - pyplot.semilogy(q, modelx(q), label="x-ray (Cu)") + pyplot.semilogy(q, resn["model"](q), label="neutron") + pyplot.semilogy(q, resx["model"](q), label="x-ray (Cu)") pyplot.legend() pyplot.title(txt) pyplot.xlabel("q [Å$^{-1}$]") pyplot.ylabel("Neutron-reflectivity") pyplot.subplot(122) - pyplot.plot(*structure.sld_profile(), label="neutron", color="C0") + pyplot.plot(*resn["structure"].sld_profile(), label="neutron", color="C0") pyplot.plot(*orsopy_neutron.sld_profile(), "--", lw=2, label="", color="C0") - pyplot.plot(*structurex.sld_profile(), label="x-ray (Cu)", color="C1") + pyplot.plot(*resx["structure"].sld_profile(), label="x-ray (Cu)", color="C1") pyplot.plot(*orsopy_xray.sld_profile(), "--", lw=2, label="", color="C1") pyplot.legend() pyplot.title(txt) From 659feb062c4beb05c086ab64c221b3f0eebe0ef1 Mon Sep 17 00:00:00 2001 From: Artur Glavic Date: Fri, 5 Jun 2026 17:24:46 +0200 Subject: [PATCH 20/22] Minor fix --- orsopy/fileio/model_complex.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/orsopy/fileio/model_complex.py b/orsopy/fileio/model_complex.py index 2e7c808..c0d709a 100644 --- a/orsopy/fileio/model_complex.py +++ b/orsopy/fileio/model_complex.py @@ -251,12 +251,12 @@ def resolve_names(self, resolvable_items): self._materials.append(material) if "environment" in resolvable_items: environment = resolvable_items["environment"] - if not isinstance(self._environment, (Composit, Material)): + if not isinstance(environment, (Composit, Material)): if environment in resolvable_items: environment = resolvable_items[environment] elif environment in SPECIAL_MATERIALS: environment = SPECIAL_MATERIALS[environment] - elif isinstance(self._environment, str): + elif isinstance(environment, str): environment = Material(formula=environment) self._materials.append(environment) @@ -353,11 +353,11 @@ def resolve_names(self, resolvable_items): environment = resolvable_items["environment"] if not isinstance(environment, Material): if environment in resolvable_items: - environment = resolvable_items[self._environment] + environment = resolvable_items[environment] elif environment in SPECIAL_MATERIALS: - environment = SPECIAL_MATERIALS[self._environment] - elif isinstance(self._environment, str): - environment = Material(formula=self._environment) + environment = SPECIAL_MATERIALS[environment] + elif isinstance(environment, str): + environment = Material(formula=environment) self._materials.append(environment) def resolve_to_blocks(self) -> List[Union["Layer", "SubStackType"]]: From fc5959a6e9e46f2cd585099b45261cdc4f967f89 Mon Sep 17 00:00:00 2001 From: Artur Glavic Date: Mon, 15 Jun 2026 13:47:44 +0200 Subject: [PATCH 21/22] Modify definition for heads/tails instead of inner/outer and only define one hydration parameter per leaflet --- examples/resolve_and_plot_model.py | 29 +++++----- examples/sample_model_example_11.yml | 22 ++++---- orsopy/fileio/model_complex.py | 80 ++++++++++++++-------------- 3 files changed, 62 insertions(+), 69 deletions(-) diff --git a/examples/resolve_and_plot_model.py b/examples/resolve_and_plot_model.py index d4878de..f4a3a30 100644 --- a/examples/resolve_and_plot_model.py +++ b/examples/resolve_and_plot_model.py @@ -61,13 +61,12 @@ def resLayer(self, layer: model_language.Layer): def resLeaflet(self, leaflet: model_complex.Leaflet): apm = leaflet.apm.as_unit("angstrom^2") vm_heads = leaflet.heads.volume.as_unit("angstrom^3") - vfrac_head = 1.0 - leaflet.heads_hydration - vfrac_tail = 1.0 - leaflet.tails_hydration + vfrac = 1.0 - leaflet.hydration b_heads = leaflet.heads.get_sld(xray_energy=self.xray_energy) * vm_heads - d_heads = vm_heads / apm / vfrac_head + d_heads = vm_heads / apm / vfrac vm_tails = leaflet.tails.volume.as_unit("angstrom^3") b_tails = leaflet.tails.get_sld(xray_energy=self.xray_energy) * vm_tails - d_tails = vm_tails / apm / vfrac_tail + d_tails = vm_tails / apm / vfrac solvent = leaflet.solvent.get_sld(xray_energy=self.xray_energy) * 1e6 self.model += f"""structure |= LipidLeaflet( apm={apm / leaflet.coverage}, @@ -87,14 +86,13 @@ def resLeaflet(self, leaflet: model_complex.Leaflet): def resBilayer(self, bilayer: model_complex.Bilayer): apm = bilayer.apm.as_unit("angstrom^2") - vm_heads = bilayer.outer.volume.as_unit("angstrom^3") - vfrac_head = 1.0 - bilayer.outer_hydration - vfrac_tail = 1.0 - bilayer.inner_hydration - b_heads = bilayer.outer.get_sld(xray_energy=self.xray_energy) * vm_heads - d_heads = vm_heads / apm / vfrac_head - vm_tails = bilayer.inner.volume.as_unit("angstrom^3") - b_tails = bilayer.inner.get_sld(xray_energy=self.xray_energy) * vm_tails - d_tails = vm_tails / apm / vfrac_tail + vm_heads = bilayer.heads.volume.as_unit("angstrom^3") + vfrac = 1.0 - bilayer.hydration + b_heads = bilayer.heads.get_sld(xray_energy=self.xray_energy) * vm_heads + d_heads = vm_heads / apm / vfrac + vm_tails = bilayer.tails.volume.as_unit("angstrom^3") + b_tails = bilayer.tails.get_sld(xray_energy=self.xray_energy) * vm_tails + d_tails = vm_tails / apm / vfrac solvent = bilayer.solvent.get_sld(xray_energy=self.xray_energy) * 1e6 self.model += f"""structure |= LipidLeaflet( apm={apm / bilayer.coverage}, @@ -111,10 +109,9 @@ def resBilayer(self, bilayer: model_complex.Bilayer): head_solvent={solvent}, tail_solvent={solvent}, )\n""" - vfrac_head = 1.0 - (bilayer.outer_hydration_2 or bilayer.outer_hydration) - vfrac_tail = 1.0 - (bilayer.inner_hydration_2 or bilayer.inner_hydration) - d_heads = vm_heads / apm / vfrac_head - d_tails = vm_tails / apm / vfrac_tail + vfrac = 1.0 - (bilayer.hydration_2 or bilayer.hydration) + d_heads = vm_heads / apm / vfrac + d_tails = vm_tails / apm / vfrac self.model += f"""structure |= LipidLeaflet( apm={apm / bilayer.coverage}, b_heads={b_heads}, diff --git a/examples/sample_model_example_11.yml b/examples/sample_model_example_11.yml index 84b66f7..4a2ce6c 100644 --- a/examples/sample_model_example_11.yml +++ b/examples/sample_model_example_11.yml @@ -1,6 +1,6 @@ sample: model: - stack: "Si | DMPC | Bilayer{outer: PC, inner: DM} | biDMPC | (lbl3 | Bilayer{DMPC}) in H2O | Leaflet{rDMPC} | H2O" + stack: "Si | DMPC | Bilayer{heads: PC, tails: DM} | biDMPC | (lbl3 | Bilayer{DMPC}) in H2O | Leaflet{rDMPC} | H2O" sub_stacks: DMPC: sub_stack_class: Leaflet @@ -12,29 +12,25 @@ sample: heads_first: false biDMPC: sub_stack_class: Bilayer - outer: PC - inner: + heads: PC + tails: composition: DM: 0.9 PC: 0.1 apm: {magnitude: 0.8, unit: nm^2} - outer_hydration: 0.3 - inner_hydration: 0.3 - inner_hydration_2: 0.4 - outer_hydration_2: 0.5 + hydration: 0.3 + hydration_2: 0.4 lbl3: sub_stack_class: Bilayer - outer: + heads: formula: C10H18O8NP number_density: 3.135 # 1/nm^3 - inner: + tails: formula: C26H54 number_density: 1.28 # 1/nm^3 apm: 0.9 # nm^2 - outer_hydration: 0.3 - inner_hydration: 0.3 - inner_hydration_2: 0.4 - outer_hydration_2: 0.5 + hydration: 0.3 + hydration_2: 0.4 coverage: 0.75 materials: PC: diff --git a/orsopy/fileio/model_complex.py b/orsopy/fileio/model_complex.py index c0d709a..9b384d2 100644 --- a/orsopy/fileio/model_complex.py +++ b/orsopy/fileio/model_complex.py @@ -212,8 +212,7 @@ class Leaflet(Header, LipidBase, SubStackType): tails: Union[Composit, Material, str] heads_first: Optional[bool] = True apm: Optional[Union[float, Value]] = field(default_factory=lambda: Value(0.7, unit="nm^2")) - heads_hydration: Optional[float] = 0.3 - tails_hydration: Optional[float] = 0.3 + hydration: Optional[float] = 0.3 coverage: Optional[float] = 1.0 roughness: Optional[Union[float, Value]] = None sub_stack_class: Literal["Leaflet"] = "Leaflet" @@ -270,18 +269,17 @@ def resolve_to_blocks(self) -> List[Union["Layer", "SubStackType"]]: def resolve_to_layers(self) -> List[Layer]: self.ensure_densities() - Vf_head = 1.0 - self.heads_hydration - Vf_tail = 1.0 - self.tails_hydration + Vf = 1.0 - self.hydration d_head = Value( - 1.0 / (Vf_head * self._materials[0].number_density.as_unit("1/nm^3") * self.apm.as_unit("nm^2")), "nm" + 1.0 / (Vf * self._materials[0].number_density.as_unit("1/nm^3") * self.apm.as_unit("nm^2")), "nm" ) d_tail = Value( - 1.0 / (Vf_tail * self._materials[1].number_density.as_unit("1/nm^3") * self.apm.as_unit("nm^2")), "nm" + 1.0 / (Vf * self._materials[1].number_density.as_unit("1/nm^3") * self.apm.as_unit("nm^2")), "nm" ) - m_head = self.mixed_material(self._materials[0], self._materials[2], self.heads_hydration) - m_tail = self.mixed_material(self._materials[1], self._materials[2], self.tails_hydration) + m_head = self.mixed_material(self._materials[0], self._materials[2], self.hydration) + m_tail = self.mixed_material(self._materials[1], self._materials[2], self.hydration) head = Layer(thickness=d_head, material=m_head, roughness=self.roughness) tail = Layer(thickness=d_tail, material=m_tail, roughness=self.roughness) @@ -294,20 +292,21 @@ def resolve_to_layers(self) -> List[Layer]: @dataclass(repr=False) class Bilayer(Header, LipidBase, SubStackType): """ - Building block corresponding to a bilayer of lipids with outer (heads) and inner (tails) - molecule definition. The bilayer is calculated from molecular volume, area per molecule (apm) and + Building block corresponding to a bilayer of lipids with heads and tails molecule definition, + optinallly with different lipids for the second Leaflet. + The bilayer is calculated from molecular volume, area per molecule (apm) and the level of hydration. The solvent used is either the environment set by a SubStack above or the default_solvent attribute of the ModelParameters. """ - outer: Union[Composit, Material, str] - inner: Union[Composit, Material, str] + heads: Union[Composit, Material, str] + tails: Union[Composit, Material, str] apm: Optional[Union[float, Value]] = field(default_factory=lambda: Value(0.7, unit="nm^2")) - outer_hydration: Optional[float] = 0.3 - inner_hydration: Optional[float] = 0.3 - outer_hydration_2: Optional[float] = None - inner_hydration_2: Optional[float] = None + hydration: Optional[float] = 0.3 + heads_2: Optional[Union[Composit, Material, str]] = None + tails_2: Optional[Union[Composit, Material, str]] = None + hydration_2: Optional[float] = None coverage: Optional[float] = 1.0 roughness: Optional[Union[float, Value]] = None sub_stack_class: Literal["Bilayer"] = "Bilayer" @@ -318,11 +317,10 @@ def resolve_name(cls, name): if name in cls.known_lipids: data = cls.known_lipids[name] for src, dest in [ - ("heads", "outer"), - ("tails", "inner"), + ("heads", "heads"), + ("tails", "tails"), ("apm", "apm"), - ("outer_hydration", "heads_hydration"), - ("inner_hydration", "tails_hydration"), + ("hydration", "hydration"), ("roughness", "roughness"), ]: if src in data: @@ -331,10 +329,16 @@ def resolve_name(cls, name): else: raise ValueError(f"Unknown lipid name: {name}") + @property + def solvent(self): + return self._materials[4] + def resolve_names(self, resolvable_items): self._materials = [] - for i, mi in enumerate((self.outer, self.inner)): - if isinstance(mi, Material): + for i, mi in enumerate((self.heads, self.tails, self.tails_2, self.heads_2)): + if i > 1 and mi is None: + material = self._materials[[0, 0, 1, 0][i]] # map self.heads/self.tails to self.tails_2/self.heads_2 + elif isinstance(mi, Material): material = mi elif isinstance(mi, Composit): mi.resolve_names(resolvable_items) @@ -363,38 +367,34 @@ def resolve_names(self, resolvable_items): def resolve_to_blocks(self) -> List[Union["Layer", "SubStackType"]]: self.ensure_densities() # Make sure the block includes full material data - self.inner = self._materials[1] - self.outer = self._materials[0] + self.heads = self._materials[0] + self.tails = self._materials[1] + self.tails_2 = self._materials[2] + self.heads_2 = self._materials[3] return [self] def resolve_to_layers(self) -> List[Layer]: self.ensure_densities() - Vf_head_1 = 1.0 - self.outer_hydration - Vf_tail_1 = 1.0 - self.inner_hydration - Vf_head_2 = 1.0 - (self.outer_hydration_2 or self.outer_hydration) - Vf_tail_2 = 1.0 - (self.inner_hydration_2 or self.inner_hydration) + Vf_1 = 1.0 - self.hydration + Vf_2 = 1.0 - (self.hydration_2 or self.hydration) d_head_1 = Value( - 1.0 / (Vf_head_1 * self._materials[0].number_density.as_unit("1/nm^3") * self.apm.as_unit("nm^2")), "nm" + 1.0 / (Vf_1 * self._materials[0].number_density.as_unit("1/nm^3") * self.apm.as_unit("nm^2")), "nm" ) d_head_2 = Value( - 1.0 / (Vf_head_2 * self._materials[0].number_density.as_unit("1/nm^3") * self.apm.as_unit("nm^2")), "nm" + 1.0 / (Vf_2 * self._materials[3].number_density.as_unit("1/nm^3") * self.apm.as_unit("nm^2")), "nm" ) d_tail_1 = Value( - 1.0 / (Vf_tail_1 * self._materials[1].number_density.as_unit("1/nm^3") * self.apm.as_unit("nm^2")), "nm" + 1.0 / (Vf_1 * self._materials[1].number_density.as_unit("1/nm^3") * self.apm.as_unit("nm^2")), "nm" ) d_tail_2 = Value( - 1.0 / (Vf_tail_2 * self._materials[1].number_density.as_unit("1/nm^3") * self.apm.as_unit("nm^2")), "nm" + 1.0 / (Vf_2 * self._materials[2].number_density.as_unit("1/nm^3") * self.apm.as_unit("nm^2")), "nm" ) - m_head_1 = self.mixed_material(self._materials[0], self._materials[2], self.outer_hydration) - m_head_2 = self.mixed_material( - self._materials[0], self._materials[2], self.outer_hydration_2 or self.outer_hydration - ) - m_tail_1 = self.mixed_material(self._materials[1], self._materials[2], self.inner_hydration) - m_tail_2 = self.mixed_material( - self._materials[1], self._materials[2], self.inner_hydration_2 or self.inner_hydration - ) + m_head_1 = self.mixed_material(self._materials[0], self._materials[4], self.hydration) + m_head_2 = self.mixed_material(self._materials[3], self._materials[4], self.hydration_2 or self.hydration) + m_tail_1 = self.mixed_material(self._materials[1], self._materials[4], self.hydration) + m_tail_2 = self.mixed_material(self._materials[2], self._materials[4], self.hydration_2 or self.hydration) head = Layer(thickness=d_head_1, material=m_head_1, roughness=self.roughness) tail = Layer(thickness=d_tail_1, material=m_tail_1, roughness=self.roughness) From 62d0a6292dc1d9a931e52866dbed83633e2da2c3 Mon Sep 17 00:00:00 2001 From: Artur Glavic Date: Mon, 15 Jun 2026 13:57:14 +0200 Subject: [PATCH 22/22] minor fix to keep tails_2/heads_2 None for resolve_to_block to signify symmetric bi-layer. --- examples/resolve_and_plot_model.py | 4 ++++ orsopy/fileio/model_complex.py | 10 ++++++++-- 2 files changed, 12 insertions(+), 2 deletions(-) diff --git a/examples/resolve_and_plot_model.py b/examples/resolve_and_plot_model.py index f4a3a30..25a8b31 100644 --- a/examples/resolve_and_plot_model.py +++ b/examples/resolve_and_plot_model.py @@ -110,7 +110,11 @@ def resBilayer(self, bilayer: model_complex.Bilayer): tail_solvent={solvent}, )\n""" vfrac = 1.0 - (bilayer.hydration_2 or bilayer.hydration) + if bilayer.heads_2 is not None: + b_heads = bilayer.heads_2.get_sld(xray_energy=self.xray_energy) * vm_heads d_heads = vm_heads / apm / vfrac + if bilayer.tails_2 is not None: + b_tails = bilayer.tails_2.get_sld(xray_energy=self.xray_energy) * vm_tails d_tails = vm_tails / apm / vfrac self.model += f"""structure |= LipidLeaflet( apm={apm / bilayer.coverage}, diff --git a/orsopy/fileio/model_complex.py b/orsopy/fileio/model_complex.py index 9b384d2..131e69e 100644 --- a/orsopy/fileio/model_complex.py +++ b/orsopy/fileio/model_complex.py @@ -333,6 +333,10 @@ def resolve_name(cls, name): def solvent(self): return self._materials[4] + @property + def is_symmetric(self): + return self.hydration_2 is None and self.heads_2 is None and self.tails_2 is None + def resolve_names(self, resolvable_items): self._materials = [] for i, mi in enumerate((self.heads, self.tails, self.tails_2, self.heads_2)): @@ -369,8 +373,10 @@ def resolve_to_blocks(self) -> List[Union["Layer", "SubStackType"]]: # Make sure the block includes full material data self.heads = self._materials[0] self.tails = self._materials[1] - self.tails_2 = self._materials[2] - self.heads_2 = self._materials[3] + if self.tails_2 is not None: + self.tails_2 = self._materials[2] + if self.heads_2 is not None: + self.heads_2 = self._materials[3] return [self] def resolve_to_layers(self) -> List[Layer]: