"""Calculates quantities used in other modules.
Uses tp units by default, but can read the tprc.yaml to customise.
"""
#Functions
#---------
#
# cumulate:
# sorts and cumulates.
# gaussian:
# gaussian distribution of height 1.
# lorentzian:
# area conserved lorentzian distribution.
# lifetime:
# particle lifetime.
# mfp:
# particle mean free path.
# fd_occupation:
# fermion occupation.
# be_occupation:
# boson occupation.
# dfdde:
# derivative of the Fermi-Dirac distribution.
# thermal_conductivity:
# total thermal conductivity.
# power_factor:
# power factor.
# zt:
# ZT.
# kl:
# lattice thermal conductivity for target ZT.
#
# power_factor_fromdict:
# adds power factor to dictionary.
# zt_fromdict:
# adds zt to dictionary.
# kl_fromdict:
# adds lattice thermal conductivity for target ZT to dictionary.
#
# to_tp:
# converts quantities to tp defaults from tprc.yaml.
# from_tp:
# converts quantities from tp defaults from tprc.yaml.
# interpolate:
# shrinks to smallest data size and interpolates.
#"""
import numpy as np
import tp
from scipy.constants import physical_constants
[docs]def cumulate(x, y):
"""Sorts by x and cumulates y.
Arguments
---------
x : array-like
x-values.
y : array-like
y-values.
Returns
-------
np.array
sorted x-values.
np.array
cumulated y-values.
"""
x = np.ravel(x)
xsort = x[x.argsort()]
y = np.ravel(y)
ysort = y[x.argsort()]
ycum = np.cumsum(np.ma.masked_invalid(ysort))
return xsort, ycum
[docs]def gaussian(x, x0=0, sigma=1):
"""Gaussian function with height 1 centered on x0.
Arguments
---------
x : array-like
x-values.
x0 : float
origin of function.
sigma :float
standard deviation.
Returns
-------
np.array
Gaussian
"""
x = np.array(x)
return np.exp(-np.power(x - x0, 2) / (2 * sigma**2))
[docs]def lorentzian(x, x0=0, fwhm=1):
"""Area conserved Lorentzian function centered on x0.
Arguments
---------
x : array-like
x-values.
x0 : float
origin of function.
fwhm : float
full-width at half-maximum.
Returns
-------
np.array
Lorentzian
"""
x = np.array(x)
return 0.5 * fwhm / (np.pi * ((x - x0)**2 + (0.5 * fwhm)**2))
[docs]def lifetime(gamma, use_tprc=True):
"""Calculates lifetime from imaginary self-energy (Gamma).
Arguments
---------
gamma : array-like or float
frequencies (by default in THz).
use_tprc : bool, optional
use custom unit conversions. Default: True.
Returns
-------
np.array
lifetimes (by default in s).
"""
if use_tprc:
gamma = to_tp('gamma', gamma)
with np.errstate(divide='ignore', invalid='ignore'):
lifetime = np.reciprocal(np.multiply(2e12 * 2 * np.pi, gamma))
infmax = np.nanmax(np.where(np.isinf(lifetime), 0, lifetime))
lifetime = np.where(np.isinf(lifetime), infmax, lifetime)
lifetime = np.where(np.isnan(lifetime), np.nanmin(lifetime), lifetime)
if use_tprc:
lifetime = from_tp('lifetime', lifetime)
return lifetime
[docs]def mfp(gamma, group_velocity, use_tprc=True):
"""Calculates mean free path from imaginary self-energy (Gamma) and group velocity.
Arguments
---------
gamma : array-like or float
frequencies (by default in THz).
group_velocity : array-like or float
group velocities (by default in m s-1).
use_tprc : bool, optional
use custom unit conversions. Default: True.
Returns
-------
array-like
mean free paths (by default in m).
"""
if use_tprc:
group_velocity = to_tp('group_velocity', group_velocity)
tau = lifetime(gamma, use_tprc=use_tprc)
mfp = np.multiply(np.transpose([tau,] * 3, (1,2,3,0)), group_velocity)
if use_tprc:
mfp = from_tp('mean_free_path', mfp)
return mfp
[docs]def fd_occupation(energy, temperature, fermi_level, use_tprc=True):
"""Calculates Fermi-Dirac occupation.
Assumes singly-occupied bands: double as necessary.
Arguments
---------
energy : array-like or float
energies (by default in eV).
temperature : array-like or float
temperature (by default in K).
fermi_level : array-like or float
fermi level (by default in eV).
use_tprc : bool, optional
use custom unit conversions. Default: True.
Returns
-------
array-like
occupations.
"""
kb = physical_constants['Boltzmann constant in eV/K'][0]
if use_tprc:
energy = to_tp('energy', energy)
temperature = to_tp('temperature', temperature)
fermi_level = to_tp('fermi_level', fermi_level)
occupation = (np.exp(np.divide(np.add.outer(-np.array(fermi_level),
energy),
kb * np.array(temperature)[None, :, None, None])) + 1) ** -1
if use_tprc:
occupation = from_tp('occupation', occupation)
return occupation
[docs]def be_occupation(frequency, temperature, use_tprc=True):
"""Calculates Bose-Einstein occupation.
Arguments
---------
frequency : array-like or float
frequencies (by default in THz).
temperature : array-like or float
temperature (by default in K).
use_tprc : bool, optional
use custom unit conversions. Default: True.
Returns
-------
array-like
occupations.
"""
hbar = physical_constants['Planck constant over 2 pi in eV s'][0]
kb = physical_constants['Boltzmann constant in eV/K'][0]
if use_tprc:
frequency = to_tp('frequency', frequency)
temperature = to_tp('temperature', temperature)
frequency = np.array(frequency)
occupation = np.expm1(np.divide.outer(kb * np.array(temperature),
frequency * 1e12 * hbar) ** -1) ** -1
if use_tprc:
occupation = from_tp('occupation', occupation)
return occupation
[docs]def dfdde(energy, fermi_level, temperature, amset_order=False,
use_tprc=True):
"""Derivative of the Fermi-Dirac distribution wrt energy.
Arguments
---------
energy : array-like
energies per band and k-point (by default in eV).
fermi_level : array-like
fermi level per temperature and dopant (by default in eV).
temperature : array-like
temperatures (by default in K).
amset_order : bool, optional
doping index before temperature index. Default: False.
use_tprc : bool, optional
use custom unit conversions. Default: True.
Returns
-------
np.array
derivative of the Fermi-Dirac distribution.
"""
kb = physical_constants['Boltzmann constant in eV/K'][0]
if use_tprc:
energy = to_tp('energy', energy)
fermi_level = to_tp('energy', fermi_level)
temperature = to_tp('temperature', temperature)
kbt = np.multiply(kb, temperature)
de = -np.subtract.outer(fermi_level, energy)
if amset_order:
weights = -0.25 / np.cosh(0.5 * de / kbt[None, :, None, None]) ** 2
weights = weights / kbt[None, :, None, None]
else:
weights = -0.25 / np.cosh(0.5 * de / kbt[:, None, None, None]) ** 2
weights = weights / kbt[:, None, None, None]
return weights
[docs]def thermal_conductivity(etc, ltc, use_tprc=True):
"""Calculates ZT.
Arguments
---------
etc : array-like
electronic thermal conductivities (by default in W m-1 K-1).
ltc : array-like
lattice thermal conductivities (by default in W m-1 K-1).
use_tprc : bool, optional
use custom unit conversions. Default: True.
Returns
-------
np.array
thermal conductivity.
"""
if use_tprc:
etc = to_tp('electronic_thermal_conductivity', etc)
ltc = to_tp('lattice_thermal_conductivity', ltc)
if np.ndim(etc) in [0, 1]:
tc = np.add(etc, np.array(ltc))
elif np.ndim(etc) == 2:
tc = np.add(etc, np.array(ltc)[:, None])
elif np.ndim(etc) == 3:
tc = np.add(etc, np.array(ltc)[:, :3, :3])
elif np.ndim(etc) == 4:
tc = np.add(etc, np.array(ltc)[:, None, :3, :3])
else:
raise Exception('Unexpectedly dimensionous electronic thermal conductivity!\n'
'Abort!')
if use_tprc:
tc = from_tp('thermal_conductivity', tc)
return tc
[docs]def power_factor(conductivity, seebeck, use_tprc=True):
"""Calculates power factor.
Arguments
---------
conductivity : array-like
conductivities (by default in S m-1).
seebeck : array-like
seebeck coefficients (by default in muV K-1).
use_tprc : bool, optional
use custom unit conversions. Default: True.
Returns
-------
np.array
power factors (by default in W m-1 K-2).
"""
if use_tprc:
conductivity = to_tp('conductivity', conductivity)
seebeck = to_tp('seebeck', seebeck)
pf = np.multiply(conductivity, 1e-12 * np.square(seebeck))
if use_tprc:
pf = from_tp('power_factor', pf)
return pf
[docs]def zt(conductivity, seebeck, electronic_thermal_conductivity,
lattice_thermal_conductivity, temperature, use_tprc=True):
"""Calculates ZT.
Arguments
---------
conductivity : array-like
conductivities (by default in S m-1).
seebeck : array-like
seebeck coefficients (by default in muV K-1).
electronic_thermal_conductivity : array-like
electronic thermal conductivities (by default in W m-1 K-1).
lattice_thermal_conductivity : array-like
lattice thermal conductivities (by default in W m-1 K-1).
temperature : array-like
temperatures (by default in K).
use_tprc : bool, optional
use custom unit conversions. Default: True.
Returns
-------
np.array
ZT.
"""
pf = power_factor(conductivity, seebeck, use_tprc=use_tprc)
tc = thermal_conductivity(electronic_thermal_conductivity,
lattice_thermal_conductivity, use_tprc=use_tprc)
if use_tprc:
pf = to_tp('power_factor', pf)
tc = to_tp('thermal_conductivity', tc)
temperature = to_tp('temperature', temperature)
if isinstance(pf, (float, int)):
zt = np.multiply(pf, temperature) / tc
else:
zt = np.apply_along_axis(np.multiply, 0, pf, temperature) / tc
if use_tprc:
zt = from_tp('zt', zt)
return zt
[docs]def kl(conductivity, seebeck, electronic_thermal_conductivity, zt, temperature,
use_tprc=True):
"""Calculates lattice thermal conductivity.
Arguments
---------
conductivity : array-like
conductivities (by default in S m-1).
seebeck : array-like
seebeck coefficients (by default in muV K-1).
electronic_thermal_conductivity : array-like
electronic thermal conductivities (by default in W m-1 K-1).
zt : array-like
zt.
temperature : array-like
temperatures (by default in K).
use_tprc : bool, optional
use custom unit conversions. Default: True.
Returns
-------
np.array
lattice thermal conductivities (by default in W m-1 K-1).
"""
pf = power_factor(conductivity, seebeck, use_tprc=use_tprc)
if use_tprc:
pf = to_tp('power_factor', pf)
electronic_thermal_conductivity = to_tp(
'electronic_thermal_conductivity', electronic_thermal_conductivity)
zt = to_tp('zt', zt)
temperature = to_tp('temperature', temperature)
mid = np.apply_along_axis(np.multiply, 0, pf, temperature) / zt
kl = np.subtract(mid, electronic_thermal_conductivity)
if use_tprc:
kl = from_tp('lattice_thermal_conductivity', kl)
return kl
[docs]def power_factor_fromdict(data, use_tprc=True):
"""Convenience wrapper to calculate power factor from a dictionary.
Arguments
---------
data : dict
dictionary containing:
conductivity : array-like
conductivities (by default in S m-1).
seebeck : array-like
seebeck coefficients (by default in muV K-1).
use_tprc : bool, optional
use custom unit conversions. Default: True.
Returns
-------
dict
dictionary with power factors (by default in W m-1 K-2).
"""
data['power_factor'] = power_factor(data['conductivity'],
data['seebeck'], use_tprc=use_tprc)
data['meta']['units']['power_factor'] = \
tp.settings.units(use_tprc=use_tprc)['power_factor']
data['meta']['dimensions']['power_factor'] = \
tp.settings.dimensions()['power_factor']
return data
[docs]def zt_fromdict(data, use_tprc=True):
"""Convenience wrapper to calculate ZT from a dictionary.
Arguments
---------
data : dict
dictionary containing:
conductivity : array-like
conductivities (by default in S m-1).
seebeck : array-like
seebeck coefficients (by default in muV K-1).
electronic_thermal_conductivity : array-like
electronic thermal conductivities (by default
in W m-1 K-1).
lattice_thermal_conductivity : array-like
lattice thermal conductivities (by default
in W m-1 K-1).
temperature : array-like
temperatures (by default in K).
use_tprc : bool, optional
use custom unit conversions. Default: True.
Returns
-------
dict
dictionary with ZTs.
"""
tc = 'thermal_conductivity'
etc, ltc = 'electronic_' + tc, 'lattice_' + tc
if 'power_factor' not in data:
data = power_factor_fromdict(data, use_tprc=use_tprc)
if 'thermal_conductivity' not in data:
data[tc] = thermal_conductivity(data[etc], data[ltc])
data['meta']['units'][tc] = tp.settings.units(use_tprc=use_tprc)[tc]
data['meta']['dimensions'][tc] = tp.settings.dimensions()[tc]
data['zt'] = zt(data['conductivity'], data['seebeck'], data[etc],
data[ltc], data['temperature'], use_tprc=use_tprc)
data['meta']['units']['zt'] = tp.settings.units(use_tprc=use_tprc)['zt']
data['meta']['dimensions']['zt'] = tp.settings.dimensions()['zt']
return data
[docs]def kl_fromdict(data, use_tprc=True):
"""Convenience wrapper to calculate k_latt from a dictionary.
Arguments
---------
data : dict
dictionary containing:
conductivity : array-like
conductivities (by default in S m-1).
seebeck : array-like
seebeck coefficients (by default in muV K-1).
electronic_thermal_conductivity : array-like
electronic thermal conductivities (by default
in W m-1 K-1).
zt : array-like
zt.
temperature : array-like
temperatures (by default in K).
use_tprc : bool, optional
use custom unit conversions. Default: True.
Returns
-------
dict
dictionary with lattice thermal conductivities (by default
in W m-1 K-1).
"""
q = 'lattice_thermal_conductivity'
data[q] = kl(data['conductivity'], data['seebeck'],
data['electronic_thermal_conductivity'], data['zt'],
data['temperature'], use_tprc=use_tprc)
data['meta']['units'][q] = tp.settings.units(use_tprc=use_tprc)[q]
data['meta']['dimensions'][q] = data['meta']['dimensions']['seebeck']
return data
[docs]def to_tp(name, value):
"""Converts quantity to tp default units using tprc.yaml.
Arguments
---------
value : array-like
values to be converted.
name : str
name in tprc.yaml.
Returns
-------
np.array
converted values.
"""
conversions = tp.settings.conversions()
if name in conversions and conversions[name] is not None:
value = np.divide(value, float(conversions[name]))
return value
[docs]def from_tp(name, value):
"""Converts quantity from tp default units using tprc.yaml.
Arguments
---------
value : array-like
values to be converted.
name : str
name in tprc.yaml.
Returns
-------
np.array
converted values.
"""
conversions = tp.settings.conversions()
if name in conversions and conversions[name] is not None:
value = np.multiply(value, float(conversions[name]))
return value
[docs]def interpolate(data1, data2, dependent, keys1, keys2, axis1=0, axis2=0,
kind='linear'):
"""Shrinks datasets to smallest common size and interpolates.
Arguments
---------
data(1,2) : dict
input data.
dependent : str
variable to interpolate against.
keys(1,2) : array-like or str
data keys to interpolate
axis(1,2) : int, optional
axis of the dependent variable wrt the keys. Default: 0.
kind : str, optional
interpolation kind
Returns
-------
dict
shrunk data1
dict
shrunk and interpolated data2
"""
# Future: could be rewritten to auto-detect axis like resolve
from copy import deepcopy
from scipy.interpolate import interp1d
data1 = deepcopy(data1)
data2 = deepcopy(data2)
if isinstance(keys1, str):
keys1 = [keys1]
if isinstance(keys2, str):
keys2 = [keys2]
# interpolate onto the densest dataset covered by both
dmin = np.nanmax([np.nanmin(data1[dependent]), np.nanmin(data2[dependent])])
dmax = np.nanmin([np.nanmax(data1[dependent]), np.nanmax(data2[dependent])])
invert = len(np.where((data1[dependent]<=dmax) & (data1[dependent]>=dmin))[0]) < \
len(np.where((data2[dependent]<=dmax) & (data2[dependent]>=dmin))[0])
if invert:
data1, data2 = data2, data1
keys1, keys2 = keys2, keys1
axis1, axis2 = axis2, axis1
index = np.where((np.array(data1[dependent])>=data2[dependent][0]) & \
(np.array(data1[dependent])<=data2[dependent][-1]))[0]
data1[dependent] = np.array(data1[dependent])[index]
for key in keys1:
data1[key] = np.swapaxes(data1[key], 0, axis1)
data1[key] = np.array(data1[key])[index]
data1[key] = np.swapaxes(data1[key], 0, axis1)
for key in keys2:
interp = interp1d(data2[dependent], data2[key], kind=kind, axis=axis2)
data2[key] = interp(data1[dependent])
data2[dependent] = data1[dependent]
if invert:
data1, data2 = data2, data1
keys1, keys2 = keys2, keys1
axis1, axis2 = axis2, axis1
return data1, data2