Source code for refnx.reflect.structure

"""
refnx is distributed under the following license:

Copyright (c) 2015 A. R. J. Nelson, ANSTO

Permission to use and redistribute the source code or binary forms of this
software and its documentation, with or without modification is hereby
granted provided that the above notice of copyright, these terms of use,
and the disclaimer of warranty below appear in the source code and
documentation, and that none of the names of above institutions or
authors appear in advertising or endorsement of works derived from this
software without specific prior written permission from all parties.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL
THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
DEALINGS IN THIS SOFTWARE.

"""
# -*- coding: utf-8 -*-

from collections import UserList
import numbers
import operator as op

import numpy as np
from scipy.stats import norm
from scipy.interpolate import interp1d

try:
    from refnx.reflect import _creflect as refcalc
except ImportError:
    from refnx.reflect import _reflect as refcalc

from refnx._lib import flatten
from refnx.analysis import Parameters, Parameter, possibly_create_parameter
from refnx.analysis.parameter import BaseParameter
from refnx.reflect.interface import Interface, Erf, Step
from refnx.reflect.reflect_model import get_reflect_backend, SpinChannel

# contracting the SLD profile can greatly speed a reflectivity calculation up.
contract_by_area = refcalc._contract_by_area


[docs]class Structure(UserList): """ Represents the interfacial Structure of a reflectometry sample. Successive Components are added to the Structure to construct the interface. Parameters ---------- components : sequence A sequence of Components to initialise the Structure. name : str Name of this structure solvent : refnx.reflect.Scatterer Specifies the scattering length density used for solvation. If no solvent is specified then the SLD of the solvent is assumed to be the SLD of `Structure[-1].slabs()[-1]` (after any possible slab order reversal). reverse_structure : bool If `Structure.reverse_structure` is `True` then the slab representation produced by `Structure.slabs` is reversed. The sld profile and calculated reflectivity will correspond to this reversed structure. contract : float If contract > 0 then an attempt to contract/shrink the slab representation is made. Use larger values for coarser profiles (and vice versa). A typical starting value to try might be 1.0. wavelength : float, None Wavelength the sample was measured at. Notes ----- If `Structure.reverse_structure is True` then the slab representation order is reversed. If no solvent is specified then the volume fraction of solvent in each of the Components is *assumed* to have the scattering length density of `Structure[-1].slabs()[-1]` after any possible slab order reversal. This slab corresponds to the scattering length density of the semi-infinite backing medium. Normally the reflectivity will be calculated using the Nevot-Croce approximation for Gaussian roughness between different layers. However, if individual components have non-Gaussian roughness (e.g. Tanh), then the overall reflectivity and SLD profile are calculated by micro-slicing. Micro-slicing involves calculating the specific SLD profile, dividing it up into small-slabs, and calculating the reflectivity from those. This normally takes much longer than the Nevot-Croce approximation. To speed the calculation up the `Structure.contract` property can be used. Contracting too far may mask the subtle differences between different roughness types. The profile contraction specified by this property can greatly improve calculation time for Structures created with micro-slicing. If you use this option it is recommended to check the reflectivity signal with and without contraction to ensure they are comparable. Example ------- >>> from refnx.reflect import SLD, Linear, Tanh, Interface >>> # make the materials >>> air = SLD(0, 0) >>> # overall SLD of polymer is (1.0 + 0.001j) x 10**-6 A**-2 >>> polymer = SLD(1.0 + 0.0001j) >>> si = SLD(2.07) >>> # Make the structure, s, from slabs. >>> # The polymer slab has a thickness of 200 A and a air/polymer roughness >>> # of 4 A. >>> s = air(0, 0) | polymer(200, 4) | si(0, 3) Use Linear roughness between air and polymer (rather than default Gaussian roughness). Use Tanh roughness between si and polymer. If non-default roughness is used then the reflectivity is calculated via micro-slicing - set the `contract` property to speed the calculation up. >>> s[1].interfaces = Linear() >>> s[2].interfaces = Tanh() >>> s.contract = 0.5 Create a user defined interfacial roughness based on the cumulative distribution function (CDF) of a Cauchy. >>> from scipy.stats import cauchy >>> class Cauchy(Interface): ... def __call__(self, x, loc=0, scale=1): ... return cauchy.cdf(x, loc=loc, scale=scale) >>> >>> c = Cauchy() >>> s[1].interfaces = c """ def __init__( self, components=(), name="", solvent=None, reverse_structure=False, contract=0, wavelength=None, ): super().__init__() self._name = name self._solvent = solvent self._reverse_structure = bool(reverse_structure) #: **float** if contract > 0 then an attempt to contract/shrink the #: slab representation is made. Use larger values for coarser profiles #: (and vice versa). A typical starting value to try might be 1.0. self.contract = contract # used for energy dispersive measurements. self.wavelength = wavelength # if you provide a list of components to start with, then initialise # the structure from that self.data = [c for c in components if isinstance(c, Component)] # private attribute for future work. Users should not count on this # staying around. self._spin = None def __copy__(self): s = Structure(name=self.name, solvent=self._solvent) s.data = self.data.copy() return s
[docs] def __setitem__(self, i, v): self.data[i] = v
[docs] def __str__(self): s = list() s.append(f"{'':_>80}") s.append(f"Structure: {self.name: ^15}") s.append(f"solvent: {(self._solvent)!r}") s.append(f"reverse structure: {self.reverse_structure}") s.append(f"contract: {self.contract}\n") for component in self: s.append(str(component)) return "\n".join(s)
def __repr__(self): return ( f"Structure(components={self.data!r}," f" name={self._name!r}," f" solvent={self._solvent!r}," f" reverse_structure={self._reverse_structure}," f" contract={self.contract})" )
[docs] def append(self, item): """ Append a :class:`Component` to the Structure. Parameters ---------- item: refnx.reflect.Component The component to be added. """ if isinstance(item, Scatterer): self.append(item()) return if not isinstance(item, Component): raise ValueError( "You can only add Component objects to a structure" ) super().append(item)
@property def name(self): return self._name @name.setter def name(self, name): self._name = name @property def solvent(self): if self._solvent is None: if not self.reverse_structure: solv_slab = self[-1].slabs(self) else: solv_slab = self[0].slabs(self) return SLD(complex(solv_slab[-1, 1], solv_slab[-1, 2])) else: return self._solvent @solvent.setter def solvent(self, sld): if sld is None: self._solvent = None elif isinstance(sld, Scatterer): # don't make a new SLD object, use its reference self._solvent = sld else: solv = SLD(sld) self._solvent = solv @property def reverse_structure(self): """ **bool** if `True` then the slab representation produced by :meth:`Structure.slabs` is reversed. The sld profile and calculated reflectivity will correspond to this reversed structure. """ return bool(self._reverse_structure) @reverse_structure.setter def reverse_structure(self, reverse_structure): self._reverse_structure = reverse_structure
[docs] @classmethod def from_slabs(cls, slabs): s = [] for slab in slabs: _slab = Slab(slab[0], complex(slab[1], slab[2]), slab[3]) s.append(_slab) return cls(s)
[docs] def slabs(self, **kwds): r""" Returns ------- slabs : :class:`np.ndarray` Slab representation of this structure. Has shape (N, 5). - slab[N, 0] thickness of layer N - slab[N, 1] *overall* SLD.real of layer N (material AND solvent) - slab[N, 2] *overall* SLD.imag of layer N (material AND solvent) - slab[N, 3] roughness between layer N and N-1 - slab[N, 4] volume fraction of solvent in layer N. Notes ----- If `Structure.reversed is True` then the slab representation order is reversed. The slab order is reversed before the solvation calculation is done. I.e. if `Structure.solvent == 'backing'` and `Structure.reversed is True` then the material that solvates the system is the component in `Structure[0]`, which corresponds to `Structure.slab[-1]`. """ if not len(self): return None if not ( isinstance(self.data[-1], (Slab, MixedSlab)) and isinstance(self.data[0], (Slab, MixedSlab)) ): raise ValueError( "The first and last Components in a Structure" " need to be Slabs" ) # over-ride the wavelength if "wavelength" in kwds: self.wavelength = float(kwds["wavelength"]) # Each layer can be given a different type of roughness profile # that defines transition between successive layers. # The default interface is specified by None (= Gaussian roughness) interfaces = flatten(self.interfaces) if all([i is None for i in interfaces]): # if all the interfaces are Gaussian, then simply concatenate # the default slabs property of each component. sl = [c.slabs(structure=self) for c in self.components] try: slabs = np.concatenate(sl) except ValueError: # some of slabs may be None. np can't concatenate arr and None slabs = np.concatenate([s for s in sl if s is not None]) else: # there is a non-default interfacial roughness, create a microslab # representation slabs = self._micro_slabs() # if the slab representation needs to be reversed. reverse = self.reverse_structure if reverse: roughnesses = slabs[1:, 3] slabs = np.flipud(slabs) slabs[1:, 3] = roughnesses[::-1] slabs[0, 3] = 0.0 if (slabs[:, 4] > 0).any(): # overall SLD is a weighted average of the vfs and slds # accessing self.solvent leads to overhead from object # creation. if self._solvent is not None: solv = self._solvent else: # we should always choose the solvating material to be the last # slab. If the structure is not reversed then you want the last # slab. If the structure is reversed then you should want to # use the first slab, but the code block above reverses the # slab order, so we still want the last one solv = complex(slabs[-1, 1], slabs[-1, 2]) slabs[1:-1] = self.overall_sld(slabs[1:-1], solv) if self.contract > 0: return contract_by_area(slabs, self.contract) else: return slabs
def _micro_slabs(self, slice_size=0.5): """ Creates a microslab representation of the Structure. Parameters ---------- slice_size : float Thickness of each slab in the micro-slab representation Returns ------- micro_slabs : np.ndarray The micro-slab representation of the model. See the `Structure.slabs` method for a description of the array. """ # solvate the slabs from each component sl = [c.slabs(structure=self) for c in self.components] total_slabs = np.concatenate(sl) total_slabs[1:-1] = self.overall_sld(total_slabs[1:-1], self.solvent) total_slabs[:, 0] = np.fabs(total_slabs[:, 0]) total_slabs[:, 3] = np.fabs(total_slabs[:, 3]) # interfaces between all the slabs _interfaces = self.interfaces erf_interface = Erf() i = 0 # the default Interface is None. # The Component.interfaces property may not have the same length as the # Component.slabs. Expand it so it matches the number of slabs, # otherwise the calculation of microslabs fails. for _interface, _slabs in zip(_interfaces, sl): if _interface is None or isinstance(_interface, Interface): f = _interface or erf_interface _interfaces[i] = [f] * len(_slabs) i += 1 _interfaces = list(flatten(_interfaces)) _interfaces = [erf_interface if i is None else i for i in _interfaces] # distance of each interface from the fronting interface dist = np.cumsum(total_slabs[:-1, 0]) # workout how much space the SLD profile should encompass zstart = -5.0 - 8 * total_slabs[1, 3] zend = 5.0 + dist[-1] + 8 * total_slabs[-1, 3] nsteps = int((zend - zstart) / slice_size + 1) zed = np.linspace(zstart, zend, num=nsteps) # the output arrays sld = np.ones_like(zed, dtype=float) * total_slabs[0, 1] isld = np.ones_like(zed, dtype=float) * total_slabs[0, 2] # work out the step in SLD at an interface delta_rho = total_slabs[1:, 1] - total_slabs[:-1, 1] delta_irho = total_slabs[1:, 2] - total_slabs[:-1, 2] # the RMS roughness of each step sigma = total_slabs[1:, 3] step = Step() # accumulate the SLD of each step. for i in range(len(total_slabs) - 1): f = _interfaces[i + 1] if sigma[i] == 0: f = step p = f(zed, scale=sigma[i], loc=dist[i]) sld += delta_rho[i] * p isld += delta_irho[i] * p sld[0] = total_slabs[0, 1] isld[0] = total_slabs[0, 2] sld[-1] = total_slabs[-1, 1] isld[-1] = total_slabs[-1, 2] micro_slabs = np.zeros((len(zed), 5), float) micro_slabs[:, 0] = zed[1] - zed[0] micro_slabs[:, 1] = sld micro_slabs[:, 2] = isld return micro_slabs @property def interfaces(self): """ A nested list containing the interfacial roughness types for each of the `Component`s. `len(Structure.interfaces) == len(Structure.components)` """ return [c.interfaces for c in self.components]
[docs] def overall_sld(self, slabs, solvent): """ Performs a volume fraction weighted average of the material SLD in a layer and the solvent in a layer. Parameters ---------- slabs : np.ndarray Slab representation of the layers to be averaged. solvent : complex or reflect.Scatterer SLD of solvating material. Returns ------- averaged_slabs : np.ndarray the averaged slabs. """ solv = solvent if isinstance(solvent, Scatterer): solv = solvent.complex(self.wavelength) return overall_sld(slabs, solv)
[docs] def reflectivity(self, q, threads=0): """ Calculate theoretical reflectivity of this structure Parameters ---------- q : array-like Q values (Angstrom**-1) for evaluation threads : int, optional Specifies the number of threads for parallel calculation. This option is only applicable if you are using the ``_creflect`` module. The option is ignored if using the pure python calculator, ``_reflect``. If `threads == 0` then all available processors are used. Notes ----- Normally the reflectivity will be calculated using the Nevot-Croce approximation for Gaussian roughness between different layers. However, if individual components have non-Gaussian roughness (e.g. Tanh), then the overall reflectivity and SLD profile are calculated by micro-slicing. Micro-slicing involves calculating the specific SLD profile, dividing it up into small-slabs, and calculating the reflectivity from those. This normally takes much longer than the Nevot-Croce approximation. To speed the calculation up the `Structure.contract` property can be used. """ abeles = get_reflect_backend() return abeles(q, self.slabs()[..., :4], threads=threads)
[docs] def sld_profile(self, z=None, align=0, max_delta_z=None): """ Calculates an SLD profile, as a function of distance through the interface. Parameters ---------- z : float Interfacial distance (Angstrom) measured from interface between the fronting medium and the first layer. align: int, optional Places a specified interface in the slab representation of a Structure at z = 0. Python indexing is allowed, e.g. supplying -1 will place the backing medium at z = 0. max_delta_z : {None, float}, optional If specified this will control the maximum spacing between SLD points. Only used if `z is None`. Returns ------- sld : float Scattering length density / 1e-6 Angstrom**-2 Notes ----- This can be called in vectorised fashion. """ slabs = self.slabs() if ( (slabs is None) or (len(slabs) < 2) or (not isinstance(self.data[0], Slab)) or (not isinstance(self.data[-1], Slab)) ): raise ValueError( "Structure requires fronting and backing" " Slabs in order to calculate." ) zed, sld = sld_profile(slabs, z=z, max_delta_z=max_delta_z) offset = 0 if align != 0: align = int(align) if align >= len(slabs) - 1 or align < -1 * len(slabs): raise RuntimeError( "abs(align) has to be less than len(slabs) - 1" ) # to figure out the offset you need to know the cumulative distance # to the interface slabs[0, 0] = slabs[-1, 0] = 0.0 if align >= 0: offset = np.sum(slabs[: align + 1, 0]) else: offset = np.sum(slabs[:align, 0]) return zed - offset, sld
[docs] def __ior__(self, other): """ Build a structure by `IOR`'ing Structures/Components/SLDs. Parameters ---------- other: :class:`Structure`, :class:`Component`, :class:`SLD` The object to add to the structure. Examples -------- >>> air = SLD(0, name='air') >>> sio2 = SLD(3.47, name='SiO2') >>> si = SLD(2.07, name='Si') >>> structure = air | sio2(20, 3) >>> structure |= si(0, 4) """ # self |= other if isinstance(other, Component): self.append(other) elif isinstance(other, Structure): self.extend(other.data) elif isinstance(other, Scatterer): slab = other(0, 0) self.append(slab) else: raise ValueError() return self
[docs] def __or__(self, other): """ Build a structure by `OR`'ing Structures/Components/SLDs. Parameters ---------- other: :class:`Structure`, :class:`Component`, :class:`SLD` The object to add to the structure. Examples -------- >>> air = SLD(0, name='air') >>> sio2 = SLD(3.47, name='SiO2') >>> si = SLD(2.07, name='Si') >>> structure = Structure() >>> structure = air | sio2(20, 3) | si(0, 3) """ # c = self | other p = Structure() p |= self p |= other return p
@property def components(self): """ The list of components in the sample. """ return self.data @property def parameters(self): r""" :class:`refnx.analysis.Parameters`, all the parameters associated with this structure. """ p = Parameters(name=f"Structure - {self.name}") p.extend([component.parameters for component in self.components]) if self._solvent is not None: p.append(self.solvent.parameters) return p
[docs] def logp(self): """ log-probability for the interfacial structure. Note that if a given component is present more than once in a Structure then it's log-prob will be counted twice. Returns ------- logp : float log-prior for the Structure. """ logp = 0 for component in self.components: logp += component.logp() return logp
def _generate_sld_profile_mcmc( self, samples=0, z=None, max_delta_z=None, align=0, random_state=None ): """ Yield SLD profiles from the MCMC samples contained within the varying parameters. The parameter state is altered during generation, but should be restored when the generator exits. If this generator is not exhausted then the parameter state will not be restored! Parameters ---------- samples: int How many SLD profiles to generate z : float Interfacial distance (Angstrom) measured from interface between the fronting medium and the first layer. align: int, optional Places a specified interface in the slab representation of a Structure at z = 0. Python indexing is allowed, e.g. supplying -1 will place the backing medium at z = 0. max_delta_z : {None, float}, optional If specified this will control the maximum spacing between SLD points. Only used if `z is None`. random_state : {int, np.random.Generator, None} random number generator that picks the samples Yields ------ sld_profile : np.ndarray `Structure.sld_profile` points for each of the samples. """ parameters = self.parameters.varying_parameters() saved_pars = np.array(parameters) _pgen = parameters.pgen(ngen=samples, random_state=random_state) try: for pars in _pgen: parameters.pvals = pars yield self.sld_profile( z=z, max_delta_z=max_delta_z, align=align ) finally: parameters.pvals = saved_pars
[docs] def plot(self, pvals=None, samples=0, fig=None, align=0): """ Plot the structure. Requires matplotlib be installed. Parameters ---------- pvals : np.ndarray, optional Numeric values for the Parameter's that are varying samples: number If this structures constituent parameters have been sampled, how many samples you wish to plot on the graph. fig: Figure instance, optional If `fig` is not supplied then a new figure is created. Otherwise the graph is created on the current axes on the supplied figure. align: int, optional Aligns the plotted structures around a specified interface in the slab representation of a Structure. This interface will appear at z = 0 in the sld plot. Note that Components can consist of more than a single slab, so some thought is required if the interface to be aligned around lies in the middle of a Component. Python indexing is allowed, e.g. supplying -1 will align at the backing medium. Returns ------- fig, ax : :class:`matplotlib.Figure`, :class:`matplotlib.Axes` `matplotlib` figure and axes objects. """ import matplotlib.pyplot as plt params = self.parameters if pvals is not None: params.pvals = pvals if fig is None: fig = plt.figure() ax = fig.add_subplot(111) else: ax = fig.gca() if samples > 0: # Get a number of chains, chosen randomly, and plot the model. for sld_profile in self._generate_sld_profile_mcmc( samples=samples, align=align ): ax.plot(*sld_profile, color="k", alpha=0.01) ax.plot(*self.sld_profile(align=align), color="red", zorder=20) ax.set_ylabel("SLD / 1e-6 $\\AA^{-2}$") ax.set_xlabel("z / $\\AA$") return fig, ax
def overall_sld(slabs, solvent): """ Performs a volume fraction weighted average of the material SLD in a layer and the solvent in a layer. Parameters ---------- slabs : np.ndarray Slab representation of the layers to be averaged. solvent : complex or reflect.Scatterer SLD of solvating material. Returns ------- averaged_slabs : np.ndarray the averaged slabs. """ slabs[..., 1:3] *= (1 - slabs[..., 4])[..., np.newaxis] slabs[..., 1] += solvent.real * slabs[..., 4] slabs[..., 2] += solvent.imag * slabs[..., 4] return slabs class Scatterer: """ Abstract base class for something that will have a scattering length density """ def __init__(self, name=""): self.name = name # by default energy dispersive scatterers are not energy dispersive self.dispersive = False def __str__(self): sld = complex(self) return f"SLD = {sld} x10**-6 Å**-2" def __complex__(self): raise NotImplementedError def complex(self, wavelength): """ Wavelength dispersive evaluation of a Scatterer SLD/RI Parameters ---------- wavelength : float Returns ------- sldc : complex """ return complex(self) @property def parameters(self): raise NotImplementedError def __call__(self, thick=0, rough=0, vfsolv=0): """ Create a :class:`Slab`. Parameters ---------- thick: refnx.analysis.Parameter or float Thickness of slab in Angstrom rough: refnx.analysis.Parameter or float Roughness of slab in Angstrom vfsolv: refnx.analysis.Parameter or float Volume fraction of water in slab Returns ------- slab : refnx.reflect.Slab The newly made Slab. Example -------- >>> # an SLD object representing Silicon Dioxide >>> sio2 = SLD(3.47, name='SiO2') >>> # create a Slab of SiO2 20 A in thickness, with a 3 A roughness >>> sio2_layer = sio2(20, 3) """ return Slab(thick, self, rough, name=self.name, vfsolv=vfsolv) def __or__(self, other): # c = self | other slab = self() return slab | other
[docs]class SLD(Scatterer): """ Object representing freely varying SLD of a material Parameters ---------- value : float, complex, Parameter, Parameters Scattering length density of a material. Units (10**-6 Angstrom**-2) name : str, optional Name of material. Notes ----- An SLD object can be used to create a Slab: >>> # an SLD object representing Silicon Dioxide >>> sio2 = SLD(3.47, name='SiO2') >>> # create a Slab of SiO2 20 A in thickness, with a 3 A roughness >>> sio2_layer = sio2(20, 3) The SLD object can also be made from a complex number, or from Parameters >>> sio2 = SLD(3.47+0.01j) >>> re = Parameter(3.47) >>> im = Parameter(0.01) >>> sio2 = SLD(re) >>> sio2 = SLD([re, im]) """ def __init__(self, value, name=""): super().__init__(name=name) self.imag = Parameter(0, name="%s - isld" % name) if isinstance(value, numbers.Real): self.real = Parameter(value.real, name="%s - sld" % name) elif isinstance(value, numbers.Complex): self.real = Parameter(value.real, name="%s - sld" % name) self.imag = Parameter(value.imag, name="%s - isld" % name) elif isinstance(value, SLD): self.real = value.real self.imag = value.imag elif isinstance(value, BaseParameter): self.real = value elif ( hasattr(value, "__len__") and isinstance(value[0], BaseParameter) and isinstance(value[1], BaseParameter) ): self.real = value[0] self.imag = value[1] self.real.units = self.imag.units = "10**-6 Å**-2" self._parameters = Parameters(name=name) def __repr__(self): return f"SLD([{self.real!r}, {self.imag!r}], name={self.name!r})"
[docs] def __complex__(self): sldc = complex(self.real.value, self.imag.value) return sldc
@property def parameters(self): """ :class:`refnx.analysis.Parameters` associated with this component """ self._parameters.name = self.name self._parameters.data = [self.real, self.imag] return self._parameters
[docs]class MaterialSLD(Scatterer): """ Object representing SLD of a chemical formula. You can fit the mass density of the material. Parameters ---------- formula : str Chemical formula density : float or Parameter mass density of compound in g / cm**3 probe : {'x-ray', 'neutron'}, optional Are you using neutrons or X-rays? wavelength : float, optional wavelength of radiation (Angstrom) name : str, optional Name of material Notes ----- You need to have the `periodictable` package installed to use this object. An SLD object can be used to create a Slab: >>> # a MaterialSLD object representing Silicon Dioxide >>> sio2 = MaterialSLD('SiO2', 2.2, name='SiO2') >>> # create a silica slab of SiO2 20 A in thickness, with a 3 A roughness >>> sio2_layer = sio2(20, 3) >>> # allow the mass density of the silica to vary between 2.1 and 2.3 >>> # g/cm**3 >>> sio2.density.setp(vary=True, bounds=(2.1, 2.3)) """ def __init__( self, formula, density, probe="neutron", wavelength=1.8, name="" ): import periodictable as pt super().__init__(name=name) self.__formula = pt.formula(formula) self._compound = formula self.density = possibly_create_parameter( density, name="density", units="g / cm**3" ) if probe.lower() not in ["x-ray", "neutron"]: raise RuntimeError("'probe' must be one of 'x-ray' or 'neutron'") self.probe = probe.lower() self.wavelength = wavelength self._parameters = Parameters(name=name) self.dispersive = True def __repr__(self): return ( f"MaterialSLD({self._compound!r}, " f"{self.density!r}, " f"probe={self.probe!r}, " f"wavelength={self.wavelength!r}, " f"name={self.name!r})" ) @property def formula(self): return self._compound @formula.setter def formula(self, formula): import periodictable as pt self.__formula = pt.formula(formula) self._compound = formula
[docs] def __complex__(self): import periodictable as pt if self.probe == "neutron": sldc = pt.neutron_sld( self.__formula, density=self.density.value, wavelength=self.wavelength, ) elif self.probe == "x-ray": sldc = pt.xray_sld( self.__formula, density=self.density.value, wavelength=self.wavelength, ) return complex(sldc[0], sldc[1])
[docs] def complex(self, wavelength): """ Wavelength dispersive evaluation of a Scatterer SLD/RI Parameters ---------- wavelength : float Returns ------- sldc : complex """ if wavelength is not None: bak = self.wavelength self.wavelength = wavelength sldc = complex(self) self.wavelength = bak return sldc return complex(self)
@property def parameters(self): self._parameters.data = [self.density] return self._parameters
[docs]class Component: """ A base class for describing the structure of a subset of an interface. Parameters ---------- name : str, optional The name associated with the Component Notes ----- By setting the `Component.interfaces` property one can control the type of interfacial roughness between all the layers of an interfacial profile. """ def __init__(self, name=""): self.name = name self._interfaces = None
[docs] def __or__(self, other): """ OR'ing components can create a :class:`Structure`. Parameters ---------- other: refnx.reflect.Structure, refnx.reflect.Component Combines with this component to make a Structure Returns ------- s: refnx.reflect.Structure The created Structure Examples -------- >>> air = SLD(0, name='air') >>> sio2 = SLD(3.47, name='SiO2') >>> si = SLD(2.07, name='Si') >>> structure = air | sio2(20, 3) | si(0, 3) """ # c = self | other p = Structure() p |= self p |= other return p
[docs] def __mul__(self, n): """ MUL'ing components makes them repeat. Parameters ---------- n: int How many times you want to repeat the Component Returns ------- s: refnx.reflect.Structure The created Structure """ # convert to integer, should raise an error if there's a problem n = op.index(n) if n < 1: return Structure() elif n == 1: return self else: s = Structure() s.extend([self] * n) return s
[docs] def __str__(self): return str(self.parameters)
@property def parameters(self): """ :class:`refnx.analysis.Parameters` associated with this component """ raise NotImplementedError( "A component should override the parameters property" ) @property def interfaces(self): """ The interfacial roughness type between each layer in `Component.slabs`. Should be one of {None, :class:`Interface`, or sequence of :class:`Interface`}. """ return self._interfaces @interfaces.setter def interfaces(self, interfaces): # Sentinel for default roughness. if interfaces is None: self._interfaces = None return if isinstance(interfaces, Interface): self._interfaces = interfaces return # this will raise TypeError is interfaces is not iterable _interfaces = [i for i in interfaces if isinstance(i, Interface)] if len(_interfaces) == 1: self._interfaces = _interfaces[0] return n_slabs = len(self.slabs()) if len(_interfaces) == n_slabs: self._interfaces = _interfaces else: raise ValueError( "Interface property must be set with one of:" " {None, Interface, sequence of Interface. If a" " sequence is provided it must have the same" " length as `Component.slabs`." )
[docs] def slabs(self, structure=None): """ The slab representation of this component Parameters ---------- structure : refnx.reflect.Structure The Structure hosting the Component. Returns ------- slabs : np.ndarray Slab representation of this Component. Has shape (N, 5). - slab[N, 0] thickness of layer N - slab[N, 1] SLD.real of layer N (not including solvent) - slab[N, 2] SLD.imag of layer N (not including solvent) - slab[N, 3] roughness between layer N and N-1 - slab[N, 4] volume fraction of solvent in layer N. If a Component returns None, then it doesn't have any slabs. """ raise NotImplementedError( "A component should override the slabs property" )
[docs] def logp(self): """ The log-probability that this Component adds to the total log-prior term. Do not include log-probability terms for the actual parameters, these are automatically included elsewhere. Returns ------- logp : float Log-probability """ return 0
[docs]class Slab(Component): """ A slab component has uniform SLD over its thickness. Parameters ---------- thick : refnx.analysis.Parameter or float thickness of slab (Angstrom) sld : :class:`refnx.reflect.Scatterer`, complex, or float (complex) SLD of film (/1e-6 Angstrom**-2) rough : refnx.analysis.Parameter or float roughness on top of this slab (Angstrom) name : str Name of this slab vfsolv : refnx.analysis.Parameter or float Volume fraction of solvent [0, 1] interface : {:class:`Interface`, None}, optional The type of interfacial roughness associated with the Slab. If `None`, then the default interfacial roughness is an Error function (also known as Gaussian roughness). """ def __init__(self, thick, sld, rough, name="", vfsolv=0, interface=None): super().__init__(name=name) self.thick = possibly_create_parameter( thick, name=f"{name} - thick", units="Å" ) if isinstance(sld, Scatterer): self.sld = sld else: self.sld = SLD(sld) self.rough = possibly_create_parameter( rough, name=f"{name} - rough", units="Å" ) self.vfsolv = possibly_create_parameter( vfsolv, name=f"{name} - volfrac solvent", bounds=(0.0, 1.0) ) self._parameters = Parameters(name=self.name) self.interfaces = interface def __repr__(self): return ( f"Slab({self.thick!r}, {self.sld!r}, {self.rough!r}," f" name={self.name!r}, vfsolv={self.vfsolv!r}," f" interface={self.interfaces!r})" )
[docs] def __str__(self): # sld = repr(self.sld) # # s = 'Slab: {0}\n thick = {1} Å, {2}, rough = {3} Å, # \u03D5_solv = {4}' # t = s.format(self.name, self.thick.value, sld, self.rough.value, # self.vfsolv.value) return str(self.parameters)
@property def parameters(self): """ :class:`refnx.analysis.Parameters` associated with this component """ self._parameters.name = self.name self._parameters.data = [ self.thick, self.sld.parameters, self.rough, self.vfsolv, ] return self._parameters
[docs] def slabs(self, structure=None): """ Slab representation of this component. See :class:`Component.slabs` """ # speculative shortcut to prevent a number of attribute retrievals if self.sld.dispersive: sldc = self.sld.complex(getattr(structure, "wavelength", None)) else: sldc = complex(self.sld) return np.array( [ [ self.thick.value, sldc.real, sldc.imag, self.rough.value, self.vfsolv.value, ] ], dtype=float, )
[docs]class MixedSlab(Component): """ A slab component made of several components Parameters ---------- thick : refnx.analysis.Parameter or float thickness of slab (Angstrom) sld_list : sequence of {refnx.reflect.Scatterer, complex, float} Sequence of (complex) SLDs that are contained in film (/1e-6 Angstrom**-2) vf_list : sequence of refnx.analysis.Parameter or float relative volume fractions of each of the materials contained in the film. rough : refnx.analysis.Parameter or float roughness on top of this slab (Angstrom) name : str Name of this slab vfsolv : refnx.analysis.Parameter or float Volume fraction of solvent [0, 1] interface : {:class:`Interface`, None}, optional The type of interfacial roughness associated with the Slab. If `None`, then the default interfacial roughness is an Error function (also known as Gaussian roughness). Notes ----- The SLD of this Slab is calculated using the normalised volume fractions of each of the constituent Scatterers: >>> np.sum([complex(sld) * vf / np.sum(vf_list) for sld, vf in ... zip(sld_list, vf_list)]). The overall SLD then takes into account the volume fraction of solvent, `vfsolv`. """ def __init__( self, thick, sld_list, vf_list, rough, name="", vfsolv=0, interface=None, ): super().__init__(name=name) self.thick = possibly_create_parameter( thick, name="%s - thick" % name, units="Å" ) self.sld = [] self.vf = [] self._sld_parameters = Parameters(name=f"{name} - slds") self._vf_parameters = Parameters(name=f"{name} - volfracs") i = 0 for s, v in zip(sld_list, vf_list): if isinstance(s, Scatterer): self.sld.append(s) else: self.sld.append(SLD(s)) vf = possibly_create_parameter( v, name=f"vf{i} - {name}", bounds=(0.0, 1.0) ) self.vf.append(vf) i += 1 self.vfsolv = possibly_create_parameter( vfsolv, name=f"{name} - volfrac solvent", bounds=(0.0, 1.0) ) self.rough = possibly_create_parameter( rough, name=f"{name} - rough", units="Å" ) p = Parameters(name=self.name) self._parameters = p self.interfaces = interface def __repr__(self): return ( f"MixedSlab({self.thick!r}, {self.sld!r}, {self.vf!r}," f" {self.rough!r}, vfsolv={self.vfsolv!r}, name={self.name!r}," f" interface={self.interfaces!r})" )
[docs] def __str__(self): return str(self.parameters)
@property def parameters(self): """ :class:`refnx.analysis.Parameters` associated with this component """ p = self._parameters p.name = self.name self._sld_parameters.data = [s.parameters for s in self.sld] self._vf_parameters.data = [vf for vf in self.vf] p.data = [self.thick] p.data.extend(self._sld_parameters) p.data.extend(self._vf_parameters) p.data.extend([self.vfsolv, self.rough]) return p
[docs] def slabs(self, structure=None): """ Slab representation of this component. See :class:`Component.slabs` """ vfs = np.array([vf.value for vf in self.vf]) sum_vfs = np.sum(vfs) sldc = np.sum( [ sld.complex(getattr(structure, "wavelength", None)) * vf / sum_vfs for sld, vf in zip(self.sld, vfs) ] ) return np.array( [ [ self.thick.value, sldc.real, sldc.imag, self.rough.value, self.vfsolv.value, ] ] )
[docs]class Stack(Component, UserList): r""" A series of Components to be considered as one. When part of a Structure the Stack can represent a multilayer by setting the `repeats` attribute. Parameters ---------- components : sequence A series of Components to initialise the stack with name : str Name of the Stack repeats : number, Parameter When viewed from a parent Structure the Components in this Stack will appear to be repeated `repeats` times. Internally `repeats` is rounded to the nearest integer before use, allowing it to be used as a fitting parameter. Notes ----- To add Components to the Stack you can: - initialise the object with a list of Components - utilise list methods (`extend`, `append`, `insert`, etc) - Add by `__ior__` (e.g. `stack |= component`) You can't use `__or__` to add Components to a stack (e.g. ``Stack() | component``) ORing a Stack with other Components will make a Structure. """ def __init__(self, components=(), name="", repeats=1): Component.__init__(self, name=name) UserList.__init__(self) # explicit calls without super self.repeats = possibly_create_parameter(repeats, "repeats") self.repeats.bounds.lb = 1 # if you provide a list of components to start with, then initialise # the Stack from that for c in components: if isinstance(c, Component): self.data.append(c) else: raise ValueError( "You can only initialise a Stack with Components" )
[docs] def __setitem__(self, i, v): self.data[i] = v
[docs] def __str__(self): s = list() s.append(f"{'' :=>80}") s.append(f"Stack start: {int(round(abs(self.repeats.value)))} repeats") for component in self: s.append(str(component)) s.append("Stack finish") s.append(f"{'' :=>80}") return "\n".join(s)
def __repr__(self): return ( f"Stack(name={self.name!r}," f" components={self.data!r}," f" repeats={self.repeats!r})" )
[docs] def append(self, item): """ Append a :class:`Component` to the Stack. Parameters ---------- item: refnx.reflect.Component The component to be added. """ if isinstance(item, Scatterer): self.append(item()) return if not isinstance(item, Component): raise ValueError( "You can only add Component objects to a structure" ) self.data.append(item)
[docs] def slabs(self, structure=None): """ Slab representation of this component. See :class:`Component.slabs` Notes ----- The overall set of slabs returned by this method consists of the concatenated constituent Component slabs repeated `Stack.repeats` times. """ if not len(self): return None # a sub stack member may want to know what the solvent is. if structure is not None: self.solvent = structure.solvent repeats = int(round(abs(self.repeats.value))) slabs = np.concatenate( [c.slabs(structure=self) for c in self.components] ) if repeats > 1: slabs = np.concatenate([slabs] * repeats) if hasattr(self, "solvent"): delattr(self, "solvent") return slabs
def _interfaces_get(self): repeats = int(round(abs(self.repeats.value))) interfaces = list(flatten([i.interfaces for i in self.data])) if repeats > 1: interfaces = interfaces * repeats return interfaces def _interfaces_set(self, interfaces): raise RuntimeError( "Cannot set interfaces property for a Stack" " Component. Please set the interfaces property" " for the constituent Components." ) # override the interfaces property for this subclass interfaces = property(_interfaces_get, _interfaces_set) @property def components(self): """ The list of components in the sample. """ return self.data @property def parameters(self): r""" :class:`refnx.analysis.Parameters`, all the parameters associated with this structure. """ p = Parameters(name=f"Stack - {self.name}") p.append(self.repeats) p.extend([component.parameters for component in self.components]) return p
[docs] def __ior__(self, other): """ Build a Stack by `IOR`'ing. Parameters ---------- other: :class:`Component`, :class:`Scatterer` The object to add to the structure. """ # self |= other if isinstance(other, Component): self.append(other) elif isinstance(other, Scatterer): slab = other(0, 0) self.append(slab) else: raise ValueError() return self
class _PolarisedSlab(Component): """ A slab component has uniform SLD over its thickness. Parameters ---------- thick : refnx.analysis.Parameter or float thickness of slab (Angstrom) sld : :class:`refnx.reflect.Scatterer`, complex, or float (complex) nuclear SLD of film (/1e-6 Angstrom**-2) rough : refnx.analysis.Parameter or float roughness on top of this slab (Angstrom) rhoM : refnx.analysis.Parameter or float magnetic SLD of film (/1e-6 Angstrom**-2) thetaM : refnx.analysis.Parameter or float Magnetic angle of the layer name : str Name of this slab vfsolv : refnx.analysis.Parameter or float Volume fraction of solvent [0, 1] interface : {:class:`Interface`, None}, optional The type of interfacial roughness associated with the Slab. If `None`, then the default interfacial roughness is an Error function (also known as Gaussian roughness). """ def __init__( self, thick, sld, rough, rhoM, thetaM, name="", vfsolv=0, interface=None, ): super().__init__(name=name) self.thick = possibly_create_parameter( thick, name=f"{name} - thick", units="Å" ) if isinstance(sld, Scatterer): self.sld = sld else: self.sld = SLD(sld) self.rough = possibly_create_parameter( rough, name=f"{name} - rough", units="Å" ) self.vfsolv = possibly_create_parameter( vfsolv, name=f"{name} - volfrac solvent", bounds=(0.0, 1.0) ) self.rhoM = possibly_create_parameter( rhoM, name=f"{name} - rhoM", ) self.thetaM = possibly_create_parameter( thetaM, name=f"{name} - thetaM", ) self._parameters = Parameters(name=self.name) self.interfaces = interface def __repr__(self): return ( f"Slab({self.thick!r}, {self.sld!r}, {self.rough!r}," f" name={self.name!r}, vfsolv={self.vfsolv!r}," f" rhoM={self.rhoM!r}, thetaM={self.thetaM!r}" f" interface={self.interfaces!r})" ) def __str__(self): return str(self.parameters) @property def parameters(self): """ :class:`refnx.analysis.Parameters` associated with this component """ self._parameters.name = self.name self._parameters.data = [ self.thick, self.sld.parameters, self.rough, self.vfsolv, self.rhoM, self.thetaM, ] return self._parameters def slabs(self, structure=None): """ Slab representation of this component. See :class:`Component.slabs` """ # speculative shortcut to prevent a number of attribute retrievals if self.sld.dispersive: sldc = self.sld.complex(getattr(structure, "wavelength", None)) else: sldc = complex(self.sld) mag_proj = self.rhoM.value * np.cos(np.degrees(self.thetaM.value)) if structure._spin in [SpinChannel.UP_UP, SpinChannel.UP_DOWN]: _op = op.add elif structure._spin in [SpinChannel.DOWN_UP, SpinChannel.DOWN_DOWN]: _op = op.sub else: raise ValueError( "Invalid spin state. Set Structure._spin to one of the" " refnx.reflect.SpinChannel enumerations." ) return np.array( [ [ self.thick.value, _op(sldc.real, mag_proj), sldc.imag, self.rough.value, self.vfsolv.value, ] ], dtype=float, ) def _profile_slicer(z, sld_profile, slice_size=None): """ Converts a scattering length density profile into a Structure by approximating with Slabs. Parameters ---------- z : array-like Distance (Angstrom) through the interface at which the SLD profile is given. sld_profile : array-like Scattering length density (10**-6 Angstrom**-2) at a given distance through the interface. Both the real and imaginary terms of the SLD can be provided - either by making `sld_profile` a complex array, by supplying an array with two columns (representing the real and imaginary parts). slice_size : None, float, optional if `slice_size is None` then `np.min(np.diff(z))/4` is used to determine the rough size of the created slabs. Otherwise `float(slice_size)` is used. Returns ------- structure : Structure A Structure representation of the sld profile Notes ----- `sld_profile` is quadratically interpolated to obtain equally spaced points. In testing the round trip structure->sld_profile->structure the maximum relative difference in reflectivity profiles from the original and final structures is on the order of fractions of a percent, with the largest difference around the critical edge. """ sld = np.asarray(sld_profile).astype(complex) if len(sld.shape) > 1 and sld.shape[1] == 2: sld[:, 0].imag = sld[:, 1].real sld = sld[:, 0] real_interp = interp1d(z, sld.real, kind="quadratic") imag_interp = interp1d(z, sld.imag, kind="quadratic") if slice_size is None: slice_size = np.min(np.diff(z)) / 4 else: slice_size = float(slice_size) # figure out the z values to calculate the slabs at z_min, z_max = np.min(z), np.max(z) n_steps = np.ceil((z_max - z_min) / slice_size) zeds = np.linspace(z_min, z_max, int(n_steps) + 1) # this is the true thickness of the slab slice_size = np.diff(zeds)[0] zeds -= slice_size / 2 zeds = zeds[1:] reals = real_interp(zeds) imags = imag_interp(zeds) slabs = [ Slab(slice_size, complex(real, imag), 0) for real, imag in zip(reals, imags) ] structure = Structure(name="sliced sld profile") structure.extend(slabs) return structure
[docs]def sld_profile(slabs, z=None, max_delta_z=None): """ Calculates an SLD profile, as a function of distance through the interface. Parameters ---------- slabs : :class:`np.ndarray` Slab representation of this structure. Has shape (N, 5). - slab[N, 0] thickness of layer N - slab[N, 1] *overall* SLD.real of layer N (material AND solvent) - slab[N, 2] *overall* SLD.imag of layer N (material AND solvent) - slab[N, 3] roughness between layer N and N-1 - slab[N, 4] volume fraction of solvent in layer N, *ignored* for this function z : {None, `np.ndarray`}, optional If `z is None` then the SLD profile will have 500 points (unless `max_delta_z` is specified). If `z` is an array, then the array specifies the interfacial locations at which the SLD will be evaluated. max_delta_z : {None, float}, optional If specified this will control the maximum spacing between SLD points. Only used if `z is None`. Returns ------- sld : float Scattering length density / 1e-6 Angstrom**-2 Notes ----- This can be called in vectorised fashion. """ nlayers = np.size(slabs, 0) - 2 # work on a copy of the input array layers = np.copy(slabs) layers[:, 0] = np.fabs(slabs[:, 0]) layers[:, 3] = np.fabs(slabs[:, 3]) # bounding layers should have zero thickness layers[0, 0] = layers[-1, 0] = 0 # distance of each interface from the fronting interface dist = np.cumsum(layers[:-1, 0]) zstart = -5 - 4 * np.fabs(slabs[1, 3]) zend = 5 + dist[-1] + 4 * layers[-1, 3] # workout how much space the SLD profile should encompass # (if z array not provided) if z is None: npnts = 500 if max_delta_z is not None: max_delta_z = float(max_delta_z) npnts_deltaz = int(np.ceil((zend - zstart) / max_delta_z)) + 1 npnts = max(500, npnts_deltaz) zed = np.linspace(zstart, zend, num=npnts) else: zed = np.asarray(z).astype(float, copy=False) # the output array sld = np.ones_like(zed, dtype=float) * layers[0, 1] # work out the step in SLD at an interface delta_rho = layers[1:, 1] - layers[:-1, 1] # use erf for roughness function, but step if the roughness is zero step_f = Step() erf_f = Erf() sigma = layers[1:, 3] # accumulate the SLD of each step. for i in range(nlayers + 1): f = erf_f if sigma[i] == 0: f = step_f sld += delta_rho[i] * f(zed, scale=sigma[i], loc=dist[i]) return zed, sld