"""
Component for studying lipid membranes at an interface
"""
import numpy as np
from scipy.optimize import NonlinearConstraint
from refnx.reflect import Component, SLD, ReflectModel, Structure
from refnx.analysis import possibly_create_parameter, Parameters, Parameter
from refnx.reflect.structure import overall_sld
[docs]class LipidLeaflet(Component):
r"""
Describes a lipid leaflet Component at an interface
Parameters
----------
APM: float or refnx.analysis.Parameter
b_heads: float, refnx.analysis.Parameter, complex or SLD
Sum of coherent scattering lengths of head group (Angstrom).
When an SLD is provided it is simply an easy way to provide a complex
value. `LipidLeaflet.b_heads_real` is set to `SLD.real`, etc.
vm_heads: float or refnx.analysis.Parameter
Molecular volume of head group (Angstrom**2)
thickness_heads: float or refnx.analysis.Parameter
Thickness of head group region (Angstrom)
b_tails: float, refnx.analysis.Parameter, complex or SLD
Sum of coherent scattering lengths of tail group (Angstrom).
When an SLD is provided it is simply an easy way to provide a complex
value. `LipidLeaflet.b_tails_real` is set to `SLD.real`, etc.
vm_tails: float or refnx.analysis.Parameter
Molecular volume of tail group (Angstrom**2)
thickness_tails: float or refnx.analysis.Parameter
Thickness of head group region (Angstrom)
rough_head_tail: float or refnx.analysis.Parameter
Roughness of head-tail group (Angstrom)
rough_preceding_mono: float or refnx.analysis.Parameter
Roughness between preceding component (in the fronting direction) and
the monolayer (Angstrom). If `reverse_monolayer is False` then this is
the roughness between the preceding component and the heads, if
`reverse_monolayer is True` then this is the roughness between the
preceding component and the tails.
head_solvent: None, float, complex, refnx.reflect.SLD
Solvent for the head region. If `None`, then solvation will be
performed by the parent `Structure`, using the `Structure.solvent`
attribute. Other options are coerced to an `SLD` object using
`SLD(float | complex)`. A float/complex argument is the SLD of the
solvent (10**-6 Angstrom**-2).
tail_solvent: None, float, complex, refnx.reflect.SLD
Solvent for the tail region. If `None`, then solvation will be
performed by the parent `Structure`, using the `Structure.solvent`
attribute. Other options are coerced to an `SLD` object using
`SLD(float | complex)`. A float/complex argument is the SLD of the
solvent (10**-6 Angstrom**-2).
reverse_monolayer: bool, optional
The default is to have heads closer to the fronting medium and
tails closer to the backing medium. If `reverse_monolayer is True`
then the tails will be closer to the fronting medium and heads
closer to the backing medium.
name: str, optional
The name for the component
Notes
-----
The sum of coherent scattering lengths must be in Angstroms, the volume
must be in cubic Angstroms. This is because the SLD of a tail group is
calculated as `b_tails / vm_tails * 1e6` to achieve the units
10**6 Angstrom**-2.
"""
# TODO: use SLD of head instead of b_heads, vm_heads?
def __init__(
self,
apm,
b_heads,
vm_heads,
thickness_heads,
b_tails,
vm_tails,
thickness_tails,
rough_head_tail,
rough_preceding_mono,
head_solvent=None,
tail_solvent=None,
reverse_monolayer=False,
name="",
):
"""
Parameters
----------
apm: float or Parameter
Area per molecule
b_heads: float, Parameter or complex
Sum of coherent scattering lengths of head group (Angstrom)
vm_heads: float or Parameter
Molecular volume of head group (Angstrom**3)
thickness_heads: float or Parameter
Thickness of head group region (Angstrom)
b_tails: float, Parameter or complex
Sum of coherent scattering lengths of tail group (Angstrom)
vm_tails: float or Parameter
Molecular volume of tail group (Angstrom**3)
thickness_tails: float or Parameter
Thickness of head group region (Angstrom)
rough_head_tail: float or refnx.analysis.Parameter
Roughness of head-tail group (Angstrom)
rough_preceding_mono: float or Parameter
Roughness between preceding component (in the fronting direction)
and the monolayer (Angstrom). If `reverse_monolayer is False` then
this is the roughness between the preceding component and the
heads, if `reverse_monolayer is True` then this is the roughness
between the preceding component and the tails.
head_solvent: None, float, complex, SLD
Solvent for the head region. If `None`, then solvation will be
performed by the parent `Structure`, using the `Structure.solvent`
attribute. Other options are coerced to an `SLD` object using
`SLD(float | complex)`. A float/complex argument is the SLD of the
solvent (10**-6 Angstrom**-2).
tail_solvent: None, float, complex, SLD
Solvent for the tail region. If `None`, then solvation will be
performed by the parent `Structure`, using the `Structure.solvent`
attribute. Other options are coerced to an `SLD` object using
`SLD(float | complex)`. A float/complex argument is the SLD of the
solvent (10**-6 Angstrom**-2).
reverse_monolayer: bool, optional
The default is to have heads closer to the fronting medium and
tails closer to the backing medium. If `reverse_monolayer is True`
then the tails will be closer to the fronting medium and heads
closer to the backing medium.
name: str, optional
The name for the component
"""
super().__init__()
self.apm = possibly_create_parameter(
apm, "%s - area_per_molecule" % name, units="Å**2"
)
if isinstance(b_heads, complex):
self.b_heads_real = possibly_create_parameter(
b_heads.real, name="%s - b_heads_real" % name
)
self.b_heads_imag = possibly_create_parameter(
b_heads.imag, name="%s - b_heads_imag" % name
)
elif isinstance(b_heads, SLD):
self.b_heads_real = b_heads.real
self.b_heads_imag = b_heads.imag
else:
self.b_heads_real = possibly_create_parameter(
b_heads, name="%s - b_heads_real" % name
)
self.b_heads_imag = possibly_create_parameter(
0, name="%s - b_heads_imag" % name
)
self.b_heads_real.units = self.b_heads_imag.units = "Å"
self.vm_heads = possibly_create_parameter(
vm_heads, name="%s - vm_heads" % name, units="Å**3"
)
self.thickness_heads = possibly_create_parameter(
thickness_heads, name="%s - thickness_heads" % name, units="Å"
)
if isinstance(b_tails, complex):
self.b_tails_real = possibly_create_parameter(
b_tails.real, name="%s - b_tails_real" % name
)
self.b_tails_imag = possibly_create_parameter(
b_tails.imag, name="%s - b_tails_imag" % name
)
elif isinstance(b_tails, SLD):
self.b_tails_real = b_tails.real
self.b_tails_imag = b_tails.imag
else:
self.b_tails_real = possibly_create_parameter(
b_tails, name="%s - b_tails_real" % name
)
self.b_tails_imag = possibly_create_parameter(
0, name="%s - b_tails_imag" % name
)
self.b_tails_real.units = self.b_tails_imag.units = "Å"
self.vm_tails = possibly_create_parameter(
vm_tails, name="%s - vm_tails" % name, units="Å**3"
)
self.thickness_tails = possibly_create_parameter(
thickness_tails, name="%s - thickness_tails" % name, units="Å"
)
self.rough_head_tail = possibly_create_parameter(
rough_head_tail, name="%s - rough_head_tail" % name, units="Å"
)
self.rough_preceding_mono = possibly_create_parameter(
rough_preceding_mono,
name="%s - rough_fronting_mono" % name,
units="Å",
)
self.head_solvent = self.tail_solvent = None
if head_solvent is not None:
self.head_solvent = SLD(head_solvent)
if tail_solvent is not None:
self.tail_solvent = SLD(tail_solvent)
self.reverse_monolayer = reverse_monolayer
self.name = name
def __repr__(self):
sld_bh = SLD([self.b_heads_real, self.b_heads_imag])
sld_bt = SLD([self.b_tails_real, self.b_tails_imag])
s = (
f"LipidLeaflet("
f"{self.apm!r}, "
f"{sld_bh!r}, "
f"{self.vm_heads!r}, "
f"{self.thickness_heads!r}, "
f"{sld_bt!r}, "
f"{self.vm_tails!r}, "
f"{self.thickness_tails!r}, "
f"{self.rough_head_tail!r}, "
f"{self.rough_preceding_mono!r}, "
f"head_solvent={self.head_solvent!r}, "
f"tail_solvent={self.tail_solvent!r}, "
f"reverse_monolayer={self.reverse_monolayer}, "
f"name={self.name!r})"
)
return s
[docs] def slabs(self, structure=None):
"""
Slab representation of monolayer, as an array
Parameters
----------
structure : refnx.reflect.Structure
The Structure hosting this Component
"""
layers = np.zeros((2, 5))
# thicknesses
layers[0, 0] = float(self.thickness_heads)
layers[1, 0] = float(self.thickness_tails)
# real and imag SLD's
layers[0, 1] = float(self.b_heads_real) / float(self.vm_heads) * 1.0e6
layers[0, 2] = float(self.b_heads_imag) / float(self.vm_heads) * 1.0e6
layers[1, 1] = float(self.b_tails_real) / float(self.vm_tails) * 1.0e6
layers[1, 2] = float(self.b_tails_imag) / float(self.vm_tails) * 1.0e6
# roughnesses
layers[0, 3] = float(self.rough_preceding_mono)
layers[1, 3] = float(self.rough_head_tail)
# volume fractions
# head region
layers[0, 4] = 1 - self.volfrac_h
if self.head_solvent is not None:
_head_solvent = self.head_solvent.complex(
getattr(structure, "wavelength", None)
)
# we do the solvation here, not in Structure.slabs
layers[0] = overall_sld(layers[0], _head_solvent)
layers[0, 4] = 0
# tail region
layers[1, 4] = 1 - self.volfrac_t
if self.tail_solvent is not None:
_tail_solvent = self.tail_solvent.complex(
getattr(structure, "wavelength", None)
)
# we do the solvation here, not in Structure.slabs
layers[1] = overall_sld(layers[1], _tail_solvent)
layers[1, 4] = 0
if self.reverse_monolayer:
layers = np.flipud(layers)
layers[:, 3] = layers[::-1, 3]
return layers
@property
def parameters(self):
p = Parameters(name=self.name)
p.extend(
[
self.apm,
self.b_heads_real,
self.b_heads_imag,
self.vm_heads,
self.thickness_heads,
self.b_tails_real,
self.b_tails_imag,
self.vm_tails,
self.thickness_tails,
self.rough_head_tail,
self.rough_preceding_mono,
]
)
if self.head_solvent is not None:
p.append(self.head_solvent.parameters)
if self.tail_solvent is not None:
p.append(self.tail_solvent.parameters)
return p
[docs] def logp(self):
# penalise unphysical volume fractions.
if self.volfrac_h > 1 or self.volfrac_t > 1:
return -np.inf
return 0
@property
def volfrac_h(self):
# Volume fraction of head group in head group region
return self.vm_heads.value / (
self.apm.value * self.thickness_heads.value
)
@property
def volfrac_t(self):
# Volume fraction of tail group in tail group region
return self.vm_tails.value / (
self.apm.value * self.thickness_tails.value
)
[docs] def make_constraint(self, objective):
"""
Creates a NonlinearConstraint for a LipidLeaflet, ensuring that volume
fraction of lipid in the head+tail regions lies in [0, 1]. Suitable for
use by differential_evolution.
Parameters
----------
objective: refnx.analysis.Objective
Objective containing the LipidLeaflet. Must be the Objective that is
being minimised by differential_evolution.
Returns
-------
nlc: NonlinearConstraint
Notes
-----
You must create separate constraints for each LipidLeaflet object in your
system.
The Objective you supply must be for the overall curve fitting system.
i.e. possibly a GlobalObjective.
Examples
--------
>>> # leaflet is a LipidLeaflet, used in an Objective, obj
>>> con = leaflet.make_constraint(obj)
>>> fitter = CurveFitter(obj)
>>> fitter.fit("differential_evolution", constraints=(con,))
"""
def con(x):
objective.setp(x)
return self.volfrac_h, self.volfrac_t
return NonlinearConstraint(con, 0, 1)