diff --git a/examples/resolve_and_plot_model.py b/examples/resolve_and_plot_model.py index 11b62c4..25a8b31 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,9 +18,120 @@ from matplotlib import pyplot from numpy import linspace -from refnx.reflect import SLD, ReflectModel, Structure -from orsopy.fileio import model_language +from orsopy.fileio import model_complex, model_language + + +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 = 1.0 - leaflet.hydration + b_heads = leaflet.heads.get_sld(xray_energy=self.xray_energy) * vm_heads + 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 + 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.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}, + 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 = 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}, + 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""" + q = linspace(0.001, 0.2, 200) @@ -32,7 +145,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) @@ -40,31 +153,51 @@ def main(txt=None): print(repr(sample), "\n") print("\n".join([repr(ss) for ss in sample.resolve_stack()]), "\n") + 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])) - 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")) - 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")) - modelx = ReflectModel(structurex, bkg=0.0) + # high-level resolution based on blocks + 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 + 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") - pyplot.plot(*structurex.sld_profile(), label="x-ray (Cu)") + pyplot.plot(*resn["structure"].sld_profile(), label="neutron", color="C0") + pyplot.plot(*orsopy_neutron.sld_profile(), "--", lw=2, label="", color="C0") + 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) pyplot.ylabel("SLD / $10^{-6} \\AA^{-2}$") diff --git a/examples/sample_model_example_10.yml b/examples/sample_model_example_10.yml new file mode 100644 index 0000000..0986e0b --- /dev/null +++ b/examples/sample_model_example_10.yml @@ -0,0 +1,11 @@ +sample: + model: + 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: + stacks: 11 + blocks: 11 + layers: 34 diff --git a/examples/sample_model_example_11.yml b/examples/sample_model_example_11.yml new file mode 100644 index 0000000..4a2ce6c --- /dev/null +++ b/examples/sample_model_example_11.yml @@ -0,0 +1,49 @@ +sample: + model: + stack: "Si | DMPC | Bilayer{heads: PC, tails: 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 + heads: PC + tails: + composition: + DM: 0.9 + PC: 0.1 + apm: {magnitude: 0.8, unit: nm^2} + hydration: 0.3 + hydration_2: 0.4 + lbl3: + sub_stack_class: Bilayer + heads: + formula: C10H18O8NP + number_density: 3.135 # 1/nm^3 + tails: + formula: C26H54 + number_density: 1.28 # 1/nm^3 + apm: 0.9 # nm^2 + hydration: 0.3 + hydration_2: 0.4 + 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: 7 + blocks: 8 + layers: 22 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 f4e431b..ce9d9a6 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,28 +9,96 @@ 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 + if hasattr(sub_stack_class, "thickness"): + self.thickness = None + + 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 + 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) + + 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"] @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] + @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 +@dataclass(repr=False) class Material(Header): formula: Optional[str] = None mass_density: Optional[Union[float, Value]] = None @@ -68,9 +136,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 @@ -114,6 +194,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 @@ -147,21 +240,7 @@ 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 +@dataclass(repr=False) class Composit(Header): composition: Dict[str, float] @@ -184,31 +263,71 @@ 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) 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 = Formula([]) + for key, value in self.composition.items(): + mi = self._composition_materials[key] + 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] + 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(repr=False) +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")), @@ -217,7 +336,7 @@ def get_sld(self, xray_energy=None): 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 b82cccc..131e69e 100644 --- a/orsopy/fileio/model_complex.py +++ b/orsopy/fileio/model_complex.py @@ -5,14 +5,16 @@ 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 -@dataclass +@dataclass(repr=False) class FunctionTwoElements(Header, SubStackType): """ Models a continuous variation between two materials/SLDs according to an analytical function. @@ -108,3 +110,300 @@ def resolve_to_layers(self) -> List[Layer]: output.append(Layer(material=composition, thickness=thickness, roughness=roughness)) output[0].roughness = self.roughness 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) + + +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] + + known_lipids = { + "DMPC": { + "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")), + }, + } + + +@dataclass(repr=False) +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 + 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[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")) + hydration: Optional[float] = 0.3 + coverage: Optional[float] = 1.0 + 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)): + 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 + 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: + environment = resolvable_items["environment"] + 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(environment, str): + 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 + self.ensure_densities() + self.heads = self._materials[0] + self.tails = self._materials[1] + return [self] + + def resolve_to_layers(self) -> List[Layer]: + self.ensure_densities() + + Vf = 1.0 - self.hydration + + d_head = Value( + 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 * 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.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) + if self.heads_first: + return [head, tail] + else: + return [tail, head] + + +@dataclass(repr=False) +class Bilayer(Header, LipidBase, SubStackType): + """ + 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. + """ + + 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")) + 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" + + @classmethod + def resolve_name(cls, name): + kwds = {} + if name in cls.known_lipids: + data = cls.known_lipids[name] + for src, dest in [ + ("heads", "heads"), + ("tails", "tails"), + ("apm", "apm"), + ("hydration", "hydration"), + ("roughness", "roughness"), + ]: + if src in data: + kwds[dest] = data[src] + return cls.from_dict(kwds) + else: + raise ValueError(f"Unknown lipid name: {name}") + + @property + 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)): + 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) + 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: + environment = resolvable_items["environment"] + if not isinstance(environment, Material): + if environment in resolvable_items: + environment = resolvable_items[environment] + elif environment in SPECIAL_MATERIALS: + 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"]]: + self.ensure_densities() + # Make sure the block includes full material data + self.heads = self._materials[0] + self.tails = self._materials[1] + 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]: + self.ensure_densities() + + Vf_1 = 1.0 - self.hydration + Vf_2 = 1.0 - (self.hydration_2 or self.hydration) + d_head_1 = Value( + 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_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_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_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[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) + 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] diff --git a/orsopy/fileio/model_language.py b/orsopy/fileio/model_language.py index e859d62..ae44ba2 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 @@ -40,7 +42,7 @@ def find_closing(string, start): return -1 -@dataclass +@dataclass(repr=False) class SubStack(Header, SubStackType): repetitions: int = 1 stack: Optional[str] = None @@ -87,6 +89,30 @@ 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 + 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: + 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() @@ -192,7 +218,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 @@ -204,7 +230,7 @@ class ItemChanger(Header): original_name = None -@dataclass +@dataclass(repr=False) class SampleModel(Header): stack: str origin: Optional[str] = None @@ -279,6 +305,30 @@ 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 + 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: + 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() 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