Lipid bilayer model

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,39 +145,59 @@ 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)
     # initial model before resolving any names
     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}$")

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

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

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.