# This program is in the public domain
# Author Paul Kienzle
"""
Reflectometry materials.
Materials (see :class:`Material`) have a composition and a density.
Density may not be known, either because it has not been measured or
because the measurement of the bulk value does not apply to thin films.
The density parameter can be fitted directly, or the bulk density can be
used, and a stretch parameter can be fitted.
Mixtures (see :class:`Mixture`) are a special kind of material which are
composed of individual parts in proportion. A mixture can be constructed
in a number of ways, such as by measuring proportional masses and mixing
or measuring proportional volumes and mixing. The parameter of interest
may also be the relative number of atoms of one material versus another.
The fractions of the different mixture components are fitted parameters,
with the remainder of the bulk filled by the final component.
SLDs (see :class:`SLD`) are raw scattering length density values. These
should be used if the material composition is not known. In that case,
you will need separate SLD materials for each wavelength and probe.
*air* (see :class:`Vacuum`) is a predefined scatterer transparent to
all probes.
Scatter (see :class:`Scatterer`) is the abstract base class from which
all scatterers are derived.
The probe cache (see :class:`ProbeCache`) stores the scattering factors
for the various materials and calls the material sld method on demand.
Because the same material can be used for multiple measurements, the
scattering factors cannot be stored with material itself, nor does it
make sense to store them with the probe. The scattering factor lookup
for the material is separate from the scattering length density
calculation so that you only need to look up the material once per fit.
The probe itself deals with all computations relating to the radiation
type and energy. Unlike the normally tabulated scattering factors f', f''
for X-ray, there is no need to scale by probe by electron radius. In
the end, sld is just the returned scattering factors times density.
"""
__all__ = ['Material', 'Mixture', 'SLD', 'Vacuum', 'Scatterer', 'ProbeCache']
import numpy as np
from numpy import inf, NaN
import periodictable
from periodictable.constants import avogadro_number
from bumps.parameter import Parameter, to_dict
[docs]
class Scatterer(object):
"""
A generic scatterer separates the lookup of the scattering factors
from the calculation of the scattering length density. This allows
programs to fit density and alloy composition more efficiently.
.. Note::
the Scatterer base class is extended by
:class:`_MaterialStacker <refl1d.model._MaterialStacker>` so that materials
can be implicitly converted to slabs when used in stack construction
expressions. It is not done directly to avoid circular dependencies
between :mod:`model <refl1d.model>` and :mod:`material <refl1d.material>`.
"""
name = None
[docs]
def sld(self, sf):
"""
Return the scattering length density expected for the given
scattering factors, as returned from a call to scattering_factors()
for a particular probe.
"""
raise NotImplementedError()
def __or__(self, other):
"""
Interface for a material stacker, to support e.g., *stack = Si | air*.
"""
raise NotImplementedError("failed monkey-patch: material stacker needs"
" to replace __or__ in Scatterer")
def __call__(self, *args, **kw):
"""
Interface for a material stacker, to support e.g., *stack = Si(thickness=)*.
"""
raise NotImplementedError("failed monkey-patch: material stacker needs"
" to replace __call__ in Scatterer")
def __str__(self):
return self.name
# ============================ No scatterer =============================
[docs]
class Vacuum(Scatterer):
"""
Empty layer
"""
name = 'air'
[docs]
def parameters(self):
return []
[docs]
def to_dict(self):
return {
'type': type(self).__name__,
}
[docs]
def sld(self, probe):
return 0, 0
def __repr__(self):
return "Vacuum()"
# ============================ Unknown scatterer ========================
[docs]
class SLD(Scatterer):
r"""
Unknown composition.
Use this when you don't know the composition of the sample. The
absorption and scattering length density are stored directly rather
than trying to guess at the composition from details about the sample.
The complex scattering potential is defined by $\rho + j \rho_i$.
Note that this differs from $\rho + j \mu/(2 \lambda)$ more
traditionally used in neutron reflectometry, and $N r_e (f_1 + j f_2)$
traditionally used in X-ray reflectometry.
Given that $f_1$ and $f_2$ are always wavelength dependent for X-ray
reflectometry, it will not make much sense to uses this for wavelength
varying X-ray measurements. Similarly, some isotopes, particularly
rare earths, show wavelength dependence for neutrons, and so
time-of-flight measurements should not be fit with a fixed SLD scatterer.
"""
def __init__(self, name="SLD", rho=0, irho=0):
self.name = name
self.rho = Parameter.default(rho, name=name+" rho")
self.irho = Parameter.default(irho, name=name+" irho")
[docs]
def parameters(self):
return {'rho':self.rho, 'irho':self.irho}
[docs]
def to_dict(self):
return to_dict({
'type': type(self).__name__,
'name': self.name,
'rho': self.rho,
'irho': self.irho,
})
[docs]
def sld(self, probe):
return self.rho.value, self.irho.value
# ============================ Substances =============================
[docs]
class Material(Scatterer):
"""
Description of a solid block of material.
:Parameters:
*formula* : Formula
Composition can be initialized from either a string or a chemical
formula. Valid values are defined in periodictable.formula.
*density* : float | |g/cm^3|
If specified, set the bulk density for the material.
*natural_density* : float | |g/cm^3|
If specified, set the natural bulk density for the material.
*use_incoherent* = False : boolean
True if incoherent scattering should be interpreted as absorption.
*fitby* = 'bulk_density' : string
Which density parameter is the fitting parameter. The choices
are *bulk_density*, *natural_density*, *relative_density* or
*cell_volume*. See :meth:`fitby` for details.
*value* : Parameter or float | units depends on fitby type
Initial value for the fitted density parameter. If None, the
value will be initialized from the material density.
For example, to fit Pd by cell volume use::
>>> m = Material('Pd', fitby='cell_volume')
>>> m.cell_volume.range(1, 10)
Parameter(Pd cell volume)
>>> print("%.2f %.2f"%(m.density.value, m.cell_volume.value))
12.02 14.70
You can change density representation by calling *material.fitby(type)*.
"""
def __init__(self, formula=None, name=None, use_incoherent=False,
density=None, natural_density=None,
fitby='bulk_density', value=None):
self.formula = periodictable.formula(formula, density=density,
natural_density=natural_density)
self.name = name if name is not None else str(self.formula)
self.use_incoherent = use_incoherent
self.fitby(type=fitby, value=value)
[docs]
def fitby(self, type='bulk_density', value=None):
"""
Specify the fitting parameter to use for material density.
:Parameters:
*type* : string
Density representation
*value* : Parameter
Initial value, or associated parameter.
Density type can be one of the following:
*bulk_density* : |g/cm^3| or kg/L
Density is *bulk_density*
*natural_density* : |g/cm^3| or kg/L
Density is *natural_density* / (natural mass/isotope mass)
*relative_density* : unitless
Density is *relative_density* * formula density
*cell_volume* : |Ang^3|
Density is mass / *cell_volume*
*number_density*: [atoms/cm^3]
Density is *number_density* * molar mass / avogadro constant
The resulting material will have a *density* attribute with the
computed material density in addition to the *fitby*
attribute specified.
.. Note::
Calling *fitby* replaces the *density* parameter in the
material, so be sure to do so before using *density* in a
parameter expression. Using *bumps.parameter.WrappedParameter*
for *density* is another alternative.
"""
# Clean out old parameter
for attr in ('bulk_density', 'natural_density', 'cell_volume',
'relative_density', 'number_density'):
try:
delattr(self, attr)
except Exception:
pass
# Put in new parameters
if type == 'bulk_density':
if value is None:
value = self.formula.density
self.bulk_density = Parameter.default(
value, name=self.name+" density", limits=(0, inf))
self.density = self.bulk_density
elif type == "number_density":
if value is None:
value = avogadro_number / self.formula.mass * self.formula.density
self.number_density = Parameter.default(
value, name=self.name+" number density", limits=(0, inf))
self.density = self.number_density / avogadro_number * self.formula.mass
elif type == 'natural_density':
if value is None:
value = self.formula.natural_density
self.natural_density = Parameter.default(
value, name=self.name+" nat. density", limits=(0, inf))
self.density = self.natural_density / self.formula.natural_mass_ratio()
elif type == 'relative_density':
if value is None:
value = 1
self.relative_density = Parameter.default(
value, name=self.name+" rel. density", limits=(0, inf))
self.density = self.formula.density*self.relative_density
## packing factor code should be correct, but radii are unreliable
#elif type is 'packing_factor':
# max_density = self.formula.mass/self.formula.volume(packing_factor=1)
# if value is None:
# value = self.formula.density/max_density
# self.packing_factor = Parameter.default(
# value, name=self.name+" packing factor")
# self.density = self.packing_factor * max_density
elif type == 'cell_volume':
# Density is in grams/cm^3.
# Mass is in grams.
# Volume is in A^3 = 1e24*cm^3.
if value is None:
value = (1e24*self.formula.molecular_mass)/self.formula.density
self.cell_volume = Parameter.default(
value, name=self.name+" cell volume", limits=(0, inf))
self.density = (1e24*self.formula.molecular_mass)/self.cell_volume
else:
raise ValueError("Unknown density calculation type '%s'"%type)
[docs]
def parameters(self):
return {'density': self.density}
[docs]
def to_dict(self):
return to_dict({
'type': type(self).__name__,
'name': self.name,
'formula': str(self.formula),
'density': self.density,
'use_incoherent': self.use_incoherent,
# TODO: what about fitby, natural_density and cell_volume?
})
[docs]
def sld(self, probe):
rho, irho, incoh = probe.scattering_factors(
self.formula, density=self.density.value)
if self.use_incoherent:
raise NotImplementedError("incoherent scattering not supported")
#irho += incoh
return rho, irho
def __str__(self):
return self.name
def __repr__(self):
return "Material(%s)"%self.name
class Compound(Scatterer):
"""
Chemical formula with variable composition.
:Parameters:
*parts* : [M1, F1, M2, F2, ...]
Unlike a simple material which has a chemical formula and a density,
the formula itself can be varied by fitting the number of atoms of
each component in the unit cell in addition to the overall density.
An individual component can be a chemical formula, not just an element.
"""
def __init__(self, parts=None):
# Split [M1, N1, M2, N2, ...] into [M1, M2, ...], [N1, N2, ...]
formula = [parts[i] for i in range(0, len(parts), 2)]
count = [parts[i] for i in range(1, len(parts), 2)]
# Convert M1, M2, ... to materials if necessary
formula = [periodictable.formula(p) for p in formula]
count = [Parameter.default(w, limits=(0, inf), name=str(f)+" count")
for w, f in zip(count, formula)]
self.parts = formula
self.count = count
def parameters(self):
"""
Adjustable parameters are the fractions associated with each
constituent and the relative scale fraction used to tweak
the overall density.
"""
return {'count': self.count}
def to_dict(self):
return {
'type': type(self).__name__,
'parts': to_dict(self.parts),
'count': to_dict(self.count),
}
def formula(self):
return tuple((c.value, f) for c, f in zip(self.count, self.parts))
def __str__(self):
return "<%s>"%(", ".join(str(M) for M in self.formula()))
def __repr__(self):
return "Compound([%s])"%(", ".join(repr(M) for M in self.formula()))
# ============================ Alloys =============================
class _VolumeFraction(object):
"""
Returns the relative volume for each component in the system given
the volume percentages.
"""
def __init__(self, base, material):
pass
def __call__(self, fraction):
return 0.01*np.asarray(fraction)
class _MassFraction(object):
"""
Returns the relative volume for each component in the system given
the relative masses.
"""
def __init__(self, base, material):
self._material = [base] + material
def __call__(self, fraction):
density = np.array([m.density.value for m in self._material])
volume = fraction/density
return volume/sum(volume)
[docs]
class Mixture(Scatterer):
"""
Mixed block of material.
The components of the mixture can vary relative to each other, either
by mass, by volume or by number::
Mixture.bymass(base, M2, F2..., name='mixture name')
Mixture.byvolume(base, M2, F2..., name='mixture name')
The materials *base*, *M2*, *M3*, ... can be chemical formula
strings including @density or from material objects. Use natural_density
to change from bulk values if the formula has isotope substitutions.
The fractions F2, F3, ... are percentages in [0, 100]. The implicit
fraction for the base material is 100 - (F2+F3+...). The SLD is NaN
then *F1 < 0*.
name defaults to M2.name+M3.name+...
"""
[docs]
@classmethod
def bymass(cls, base, *parts, **kw):
"""
Returns an alloy defined by relative mass of the constituents.
Mixture.bymass(base, M2, F2, ..., name='mixture name')
"""
return cls(base, parts, by='mass', **kw)
[docs]
@classmethod
def byvolume(cls, base, *parts, **kw):
"""
Returns an alloy defined by relative volume of the constituents.
Mixture.byvolume(base, M2, F2, ..., name='mixture name')
"""
return cls(base, parts, by='volume', **kw)
def __init__(self, base, parts, by='volume', name=None, use_incoherent=False):
# Split [M1, M2, F2, ...] into [M1, M2, ...], [F2, ...]
material = [parts[i] for i in range(0, len(parts), 2)]
fraction = [parts[i] for i in range(1, len(parts), 2)]
# Convert M1, M2, ... to materials if necessary
if not isinstance(base, Scatterer):
base = Material(base)
material = [p if isinstance(p, Scatterer) else Material(p)
for p in material]
# Specify the volume calculator based on the type of fraction
if by == 'volume':
_volume = _VolumeFraction(base, material)
elif by == 'mass':
_volume = _MassFraction(base, material)
else:
raise ValueError('fraction must be one of volume, mass or number')
# Name defaults to names of individual components
if name is None:
name = "+".join(p.name for p in material)
# Make the fractions into fittable parameters
fraction = [Parameter.default(w, limits=(0, 100), name=f.name+"% "+by)
for w, f in zip(fraction, material)]
self._volume = _volume
self.base = base
self.material = material
self.fraction = fraction
self.name = name
self.use_incoherent = use_incoherent
[docs]
def parameters(self):
"""
Adjustable parameters are the fractions associated with each
constituent and the relative scale fraction used to tweak
the overall density.
"""
return {
'base':self.base.parameters(),
'material':[m.parameters() for m in self.material],
'fraction':self.fraction,
}
[docs]
def to_dict(self):
return {
'type': type(self).__name__,
'base': to_dict(self.base),
'material': to_dict(self.material),
'fraction': to_dict(self.fraction),
}
def _density(self):
"""
Compute the density of the mixture from the density and proportion
of the individual components.
"""
fraction = np.array([0.]+[m.value for m in self.fraction])
# TODO: handle invalid fractions using penalty functions
# S = sum(fraction)
# scale = S/100 if S > 100 else 1
# fraction[0] = 100 - S/scale
# penalty = scale - 1
fraction[0] = 100 - sum(fraction)
if (fraction < 0).any():
return NaN
volume = self._volume(fraction)
density = np.array([m.density() for m in [self.base]+self.material])
return np.sum(volume*density)
density = property(_density, doc=_density.__doc__)
[docs]
def sld(self, probe):
"""
Return the scattering length density and absorption of the mixture.
"""
# Convert fractions into an array, with the final fraction
fraction = np.array([0.]+[f.value for f in self.fraction])
fraction[0] = 100 - sum(fraction)
# TODO: handle invalid fractions using penalty functions
# S = sum(fraction)
# scale = S/100 if S > 100 else 1
# fraction[0] = 100 - S/scale
# penalty = scale - 1
if (fraction < 0).any():
return NaN, NaN
# Lookup SLD
slds = [c.sld(probe) for c in [self.base] + self.material]
rho, irho = [np.asarray(v) for v in zip(*slds)]
# Use calculator to convert individual SLDs to overall SLD
volume_fraction = self._volume(fraction)
rho = np.sum(rho*extend(volume_fraction, rho))
irho = np.sum(irho*extend(volume_fraction, irho))
if self.use_incoherent:
raise NotImplementedError("incoherent scattering not supported")
#print "Mixture", self.name, coh, absorp
return rho, irho
def __str__(self):
return "<%s>"%(", ".join(str(M) for M in [self.base]+self.material))
def __repr__(self):
return "Mixture(%s)"%(", ".join(repr(M) for M in [self.base]+self.material))
# ============================ SLD cache =============================
[docs]
class ProbeCache(object):
"""
Probe proxy for materials properties.
A caching probe which only looks up scattering factors for materials
which it hasn't seen before. Note that caching is based on object
id, and will fail if the material object is updated with a new atomic
structure.
*probe* is the probe to use when looking up the scattering length density.
The scattering factors need to be retrieved each time the probe
or the composition changes. This can be done either by deleting
an individual material from probe (using del probe[material]) or
by clearing the entire cash.
"""
def __init__(self, probe=None):
self._probe = probe
self._cache = {}
[docs]
def clear(self):
self._cache = {}
def __delitem__(self, material):
if material in self._cache:
del self._cache[material]
[docs]
def scattering_factors(self, material, density):
"""
Return the scattering factors for the material, retrieving them from
the cache if they have already been looked up.
"""
h = id(material)
if h not in self._cache:
# lookup density of 1, and scale to actual density on retrieval
self._cache[h] = self._probe.scattering_factors(material,
density=1.0)
return [v*density for v in self._cache[h]]
def extend(a, b):
"""
Extend *a* to match the number of dimensions of *b*.
This adds dimensions to the end of *a* rather than the beginning. It is
equivalent to *a[..., None, None]* with the right number of None elements
to make the number of dimensions match (or np.newaxis if you prefer).
For example::
from numpy.random import rand
a, b = rand(3, 4), rand(3, 4, 2)
a + b
==> ValueError: operands could not be broadcast together with shapes (3,4) (3,4,2)
c = extend(a, b) + b
c.shape
==> (3, 4, 2)
Numpy broadcasting rules automatically extend arrays to the beginning,
so the corresponding *lextend* function is not needed::
c = rand(3, 4) + rand(2, 3, 4)
c.shape
==> (2, 3, 4)
"""
if np.isscalar(a):
return a
# CRUFT: python 2.7 support
extra_dims = (1,)*(b.ndim-a.ndim)
return a.reshape(a.shape + extra_dims)
# python 3 uses
#extra_dims = (np.newaxis,)*(b.ndim-a.ndim)
#return a[(..., *extra_dims)]