IPKM/nmreval
17
2
Fork 1
Code Issues 42 Pull Requests Actions Packages Projects Releases Wiki Activity

BUGFIX: VFT;

Browse Source 
change to src layout
This commit is contained in:
dominik
2022-10-20 17:23:15 +02:00
parent 89ce4bab9f
commit 8d148b639b
445 changed files with 1387 additions and 1920 deletions
Download Patch File Download Diff File Expand all files Collapse all files

28
src/nmreval/distributions/__init__.py Normal file
View File

@ -0,0 +1,28 @@
"""
=================================================================
Distributions of correlation times (:mod:`nmreval.distributions`)
=================================================================
Each class provides a distribution of correlation times, correlation
function, spectral density, susceptibility, and calculation of different
types of correlation times.
.. autosummary::
:toctree: generated/
:nosignatures:
Debye
ColeCole
ColeDavidson
HavriliakNegami
KWW
LogGaussian
"""
from .havriliaknegami import HavriliakNegami
from .colecole import ColeCole
from .coledavidson import ColeDavidson
from .debye import Debye
from .kww import KWW
from .loggaussian import LogGaussian

89
src/nmreval/distributions/base.py Normal file
View File

@ -0,0 +1,89 @@
from __future__ import annotations
import abc
from typing import Any
import numpy as np
from ..lib.utils import ArrayLike
class Distribution(abc.ABC):
name = 'Abstract distribution'
parameter = None
bounds = None
@classmethod
def __repr__(cls):
return cls.name
@staticmethod
@abc.abstractmethod
def distribution(tau, tau0, *args):
pass
@staticmethod
@abc.abstractmethod
def correlation(t, tau, *args):
pass
@staticmethod
@abc.abstractmethod
def susceptibility(omega, tau, *args):
pass
@classmethod
def specdens(cls, omega: ArrayLike, tau: ArrayLike, *args: Any) -> np.ndarray | float:
return cls.susceptibility(omega, tau, *args).imag / omega
@classmethod
def mean_value(cls, *args, mode='mean'):
if mode == 'mean':
return cls.mean(*args)
if mode == 'logmean':
return np.exp(cls.logmean(*args))
if mode in ['max', 'peak']:
return cls.max(*args)
return args[0]
@staticmethod
def mean(*args):
return args[0]
@staticmethod
def logmean(*args):
"""
Return mean logarithmic tau
Args:
*args:
Returns:
"""
return np.log(args[0])
@staticmethod
def max(*args):
return args[0]
@classmethod
def convert(cls, *args, from_='raw', to_='mean'):
if from_ == to_:
return args[0]
if from_ == 'raw':
return cls.mean_value(*args, mode=to_)
else:
# makes only sense as long as mean/peak is given by <x> = factor * x
factor = cls.mean_value(1, *args[1:], mode=from_)
try:
raw = args[0] / factor
except ZeroDivisionError:
raise ZeroDivisionError('Conversion between mean values impossible')
if to_ == 'raw':
return raw
else:
return cls.mean_value(raw, *args[1:], mode=to_)

73
src/nmreval/distributions/colecole.py Normal file
View File

@ -0,0 +1,73 @@
from __future__ import annotations
import numpy as np
from .base import Distribution
from ..lib.decorator import adjust_dims
from ..math.mittagleffler import mlf
class ColeCole(Distribution):
"""
Functions based on Cole-Cole distribution
"""
name = 'Cole-Cole'
parameter = [r'\alpha']
bounds = [(0, 1)]
@staticmethod
def distribution(tau, tau0, alpha):
"""
Distribution of correlation times
Args:
tau (array_like):
tau0 (array_like):
alpha (float):
"""
z = np.log(tau/tau0)
return np.sin(np.pi*alpha)/(np.cosh(z*alpha) + np.cos(alpha*np.pi))/(2*np.pi)
@staticmethod
@adjust_dims
def susceptibility(omega, tau, alpha: float):
r"""
Complex susceptibility
Args:
omega (array_like): Frequency axis in 1/s (not Hz).
tau (array_like): Correlation times :math:`\tau_\text{CC}` in s.
alpha (float): Shape parameter.
"""
val = 1/(1+(1j*omega*tau)**alpha)
return np.conjugate(val)
@staticmethod
@adjust_dims
def specdens(omega, tau, alpha):
"""
Spectral density
Args:
omega (array_like):
tau (array_like):
alpha (float):
"""
omtau = (omega*tau)**alpha
return np.sin(alpha*np.pi/2) * omtau / (1 + omtau**2 + 2*np.cos(alpha*np.pi/2)*omtau) / omega
@staticmethod
@adjust_dims
def correlation(t, tau0, alpha):
"""
Correlation function :math:`C(t)`
with Mittag-Leffler function :math:`E_a(x)`
Args:
t (array_like):
tau0 (array_like):
alpha (float):
"""
return mlf(-(t/tau0)**alpha, alpha)

164
src/nmreval/distributions/coledavidson.py Normal file
View File

@ -0,0 +1,164 @@
from __future__ import annotations
import numpy as np
from scipy.special import psi, gammaincc
from .base import Distribution
from ..lib.decorator import adjust_dims
from ..utils.constants import Eu
class ColeDavidson(Distribution):
"""
Functions based on Cole-Davidson distribution
"""
name = 'Cole-Davidson'
parameter = [r'\gamma']
bounds = [(0, 1)]
@staticmethod
def distribution(tau, tau0, gamma):
"""
Distribution of correlation times
.. math::
G(\ln\\tau) =
\\begin{cases}
\\frac{\sin(\pi\gamma)}{\pi} \\Big(\\frac{\\tau}{\\tau_0 - \\tau}\\Big)^\gamma & \\tau < \\tau_0 \\\\
0 & \\text{else}
\end{cases}
Args:
tau:
tau0:
gamma:
"""
ratio = tau0/tau
ret_val = np.zeros_like(ratio)
ret_val[ratio > 1] = np.sin(np.pi*gamma) / np.pi * (1 / (ratio[ratio>1] - 1))**gamma
return ret_val
@staticmethod
@adjust_dims
def susceptibility(omega: float | np.ndarray, tau: float | np.ndarray, gamma: float) -> float | np.ndarray:
r"""
Complex susceptibility
.. math::
\chi(\omega, \tau, \alpha, \gamma) = \frac{1}{[1-i\omega\tau]^\gamma}
Args:
omega (array_like): Frequency axis in 1/s (not Hz).
tau (array_like): Correlation times in s.
gamma (float): Shape parameter.
"""
return (1/(1+1j*omega*tau)**gamma).conjugate()
@staticmethod
@adjust_dims
def specdens(omega, tau, gamma):
"""
Spectral density
.. math::
J(\omega, \\tau, \\gamma) = \\frac{\sin[\\gamma\\arctan(\\omega\\tau)]}{\omega[1+(\omega\\tau)^2]^{\\gamma/2}}
Args:
omega:
tau:
gamma:
"""
omtau = omega*tau
ret_val = np.sin(gamma*np.arctan2(omtau, 1)) / (1 + omtau**2)**(gamma/2.) / omega
ret_val[np.argwhere(omega == 0)] = gamma*tau
return ret_val
@staticmethod
def mean(tau, gamma):
"""
Mean tau
.. math::
\\langle \\tau \\rangle = \\gamma \\tau
Args:
tau:
gamma:
"""
return tau*gamma
@staticmethod
def logmean(tau, gamma):
"""
Mean logarithm of tau
.. math::
\\langle \ln \\tau \\rangle = \ln\\tau + \psi(\gamma) + \mathrm{Eu}
with the Euler's constant :math:`\mathrm{Eu} \\approx 0.57721567\ldots` and the digamma function
.. math::
\psi(x) = \\frac{d}{dx} \ln(\Gamma(x))
Args:
tau:
gamma:
References:
R. Zorn: Logartihmic moments of relaxation time distribution. J. Chem. Phys. (2002). http://dx.doi.org/10.1063/1.1446035
"""
return np.log(tau) + psi(gamma) + Eu
@staticmethod
def max(tau, gamma):
"""
Calculate :math:`\\tau_\\text{peak}`, i.e. the inverse of the maximum frequency of the imaginary part of the susceptibility
.. math::
\\tau_\\text{peak} = \\frac{\\tau}{\\tan[\\pi/(2\gamma+2)]}
Args:
tau:
gamma:
References:
R. Dı́az-Calleja: Comment on the Maximum in the Loss Permittivity for the Havriliak-Negami Equation.
Macromolecules (2000). https://doi.org/10.1021/ma991082i
"""
return tau/np.tan(np.pi/(2*gamma+2))
@staticmethod
@adjust_dims
def correlation(t, tau, gamma):
r"""
Correlation function
.. math::
C(t, \tau_0, \gamma) = \frac{\Gamma(\gamma, t/\tau_0)}{\Gamma(\gamma)}
with incomplete Gamma function
.. math::
\Gamma(\gamma, t/\tau_0) = \frac{1}{\Gamma(\gamma)}\int_{t/\tau_0}^\infty x^{\gamma-1}\text{e}^{-x}\,\mathrm{d}x
Args:
t:
tau:
gamma:
References:
R. Hilfer: H-function representations for stretched exponential relaxation and non-Debye susceptibilities in glassy systems. Phys. Rev. E (2002) https://doi.org/10.1103/PhysRevE.65.061510
"""
return gammaincc(gamma, t / tau)

46
src/nmreval/distributions/debye.py Normal file
View File

@ -0,0 +1,46 @@
import numbers
import numpy as np
from .base import Distribution
from ..lib.decorator import adjust_dims
from ..lib.utils import ArrayLike
class Debye(Distribution):
name = 'Debye'
@staticmethod
def distribution(tau: ArrayLike, tau0: float) -> ArrayLike:
if isinstance(tau, numbers.Number):
return 1 if tau == tau0 else 0
ret_val = np.zeros_like(tau)
ret_val[tau == tau0] = 1.
return ret_val
@staticmethod
@adjust_dims
def correlation(t, tau, *args):
return np.exp(-t/tau)
@staticmethod
@adjust_dims
def susceptibility(omega: ArrayLike, tau: ArrayLike) -> ArrayLike:
r"""
Complex susceptibility
.. math::
\chi(\omega, \tau) = \frac{1}{1-i\omega\tau}
Args:
omega (array_like): Frequency axis in 1/s (not Hz).
tau (array_like): Correlation times in s.
"""
return 1/(1 - 1j*omega*tau)
@staticmethod
@adjust_dims
def specdens(omega: ArrayLike, tau: ArrayLike) -> ArrayLike:
return tau / (1 + (omega*tau)**2)

120
src/nmreval/distributions/energy.py Normal file
View File

@ -0,0 +1,120 @@
from itertools import product
import numpy as np
from scipy.integrate import quad, simps as simpson
from .base import Distribution
from ..lib.utils import ArrayLike
from ..utils.constants import kB
class EnergyBarriers(Distribution):
name = 'Energy barriers'
parameter = [r'\tau_{0}', r'E_{m}', r'\Delta E']
@staticmethod
def distribution(tau, temperature, *args):
t0, em, sigma = args
m = np.exp(em / (kB * temperature)) * t0
return temperature * np.exp(-0.5 * (np.log(tau/m)*kB*temperature/sigma)**2) / np.sqrt(2*np.pi)/sigma
@staticmethod
def rate(t0, e_a, te):
return np.exp(-e_a / (kB * te)) / t0
@staticmethod
def energydistribution(e_a, mu, sigma):
return np.exp(-0.5 * ((mu-e_a) / sigma) ** 2) / (np.sqrt(2 * np.pi) * sigma)
@staticmethod
def correlation(t, temperature, *args):
tau0, e_m, e_b = args
def integrand(e_a, ti, t0, mu, sigma, te):
# correlation time would go to inf for higher energies, so we use rate
return np.exp(-ti*EnergyBarriers.rate(t0, e_a, te)) * EnergyBarriers.energydistribution(e_a, mu, sigma)
t = np.atleast_1d(t)
temperature = np.atleast_1d(temperature)
e_axis = np.linspace(max(0, e_m-50*e_b), e_m+50*e_b, num=5001)
ret_val = np.array([simpson(integrand(e_axis, o, tau0, e_m, e_b, tt), e_axis)
for o in t for tt in temperature])
return ret_val
@staticmethod
def susceptibility(omega: ArrayLike, temperature: ArrayLike, tau0: float, e_m: float, e_b: float) -> ArrayLike:
# in contrast to other spectral densities, it's omega and temperature
omega = np.atleast_1d(omega)
temperature = np.atleast_1d(temperature)
e_axis = np.linspace(max(0, e_m-50*e_b), e_m+50*e_b, num=5001)
ret_val = []
for o, tt in product(omega, temperature):
ret_val.append(simpson(_integrand_freq_real(e_axis, o, tau0, e_m, e_b, tt), e_axis) -
1j * simpson(_integrand_freq_imag(e_axis, o, tau0, e_m, e_b, tt), e_axis))
return np.array(ret_val)
@staticmethod
def specdens(omega, temperature, *args):
# in contrast to other spectral densities, it's omega and temperature
tau0, e_m, e_b = args
def integrand(e_a, w, t0, mu, sigma, t):
r = EnergyBarriers.rate(t0, e_a, t)
return r/(r**2 + w**2) * EnergyBarriers.energydistribution(e_a, mu, sigma)
omega = np.atleast_1d(omega)
temperature = np.atleast_1d(temperature)
e_axis = np.linspace(max(0, e_m-50*e_b), e_m+50*e_b, num=5001)
ret_val = np.array([simpson(integrand(e_axis, o, tau0, e_m, e_b, tt), e_axis)
for o in omega for tt in temperature])
return ret_val
@staticmethod
def mean(*args):
return args[1]*np.exp(args[2]/(kB*args[0]))
@staticmethod
def logmean(*args):
return args[1] + args[2] / (kB * args[0])
@staticmethod
def max(*args):
return args[1] * np.exp(args[2] / (kB * args[0]))
def _integrate_process_imag(args):
pass
def _integrate_process_real(args):
omega_i, t, tau0, mu, sigma, temp_j = args
return quad(_integrand_freq_real(), 0, 10, args=(omega_i, t, tau0, mu, sigma, temp_j))[0]
def _integrate_process_time(args):
omega_i, t, tau0, mu, sigma, temp_j = args
return quad(_integrand_time, 0, 10, args=(omega_i, t, tau0, mu, sigma, temp_j))[0]
def _integrand_freq_real(u, omega, tau0, mu, sigma, temp):
r = EnergyBarriers.rate(tau0, u, temp)
return 1 / (r**2 + omega**2) * EnergyBarriers.energydistribution(u, mu, sigma)
def _integrand_freq_imag(u, omega, tau0, mu, sigma, temp):
rate = EnergyBarriers.rate(tau0, u, temp)
return omega * rate / (rate**2 + omega**2) * EnergyBarriers.energydistribution(u, mu, sigma)
def _integrand_time(u, t, tau0, mu, sigma, temp):
rate = EnergyBarriers.rate(tau0, u, temp)
return EnergyBarriers.energydistribution(u, mu, sigma) * np.exp(-t*rate)

178
src/nmreval/distributions/gengamma.py Normal file
View File

@ -0,0 +1,178 @@
from __future__ import annotations
from abc import ABC
import numpy as np
try:
from scipy.integrate import simpson
except ImportError:
from scipy.integrate import simps as simpson
from scipy.special import gammaln
from nmreval.distributions.base import Distribution
from nmreval.math.logfourier import logft
class AbstractGG(Distribution, ABC):
"""
Base class for general gamma functions
"""
@classmethod
def correlation(cls, t, tau0, *args):
tt = np.asanyarray(t)
taus, ln_tau = AbstractGG._prepare_integration(tau0)
g_tau = cls.distribution(taus, tau0, *args)
ret_val = np.array([simpson(np.exp(-t_i/taus) * g_tau, ln_tau) for t_i in tt])
return ret_val
@classmethod
def susceptibility(cls, omega, tau0, *args):
r"""
Calculate spectral density \int G(ln(tau) tau/(1+(w*tau)^2) dln(tau)
"""
w = np.asanyarray(omega)
taus, ln_tau = AbstractGG._prepare_integration(tau0)
g_tau = cls.distribution(taus, tau0, *args)
ret_val = np.array([simpson(g_tau / (1 - 1j*w_i*taus), ln_tau) for w_i in w])
return ret_val
@classmethod
def specdens(cls, omega, tau0, *args):
r"""
Calculate spectral density \int G(ln(tau) tau/(1+(w*tau)^2) dln(tau)
"""
w = np.asanyarray(omega)
taus, ln_tau = AbstractGG._prepare_integration(tau0)
g_tau = cls.distribution(taus, tau0, *args)
ret_val = np.array([simpson(g_tau * taus / (1 + (w_i*taus)**2), ln_tau) for w_i in w])
return ret_val
@staticmethod
def _prepare_integration(
tau0: float, limits: tuple[int, int] = (20, 20), num_steps: int = 4001
) -> tuple[np.ndarray, np.ndarray]:
"""
Create array of correlation times for integration over ln(tau)
Args:
tau0 (float): center of correlation times
limits (tuple of int): Define limits ln(tau0)-limits[0] and ln(tau0)+limits[1] of array. Default is (20, 20)
num_steps (int): length of array. Default is 4001
Returns:
array of taus and array of ln(tau)
"""
ln_tau0 = np.log(tau0)
ln_tau = np.linspace(ln_tau0-limits[0], ln_tau0+limits[1], num=num_steps)
return np.exp(ln_tau), ln_tau
# noinspection PyMethodOverriding
class GGAlpha(AbstractGG):
name = r'General \Gamma (\alpha)'
parameter = [r'\tau', r'\alpha', r'\beta']
@staticmethod
def distribution(taus: float | np.ndarray, tau: float, alpha: float, beta: float) -> float | np.ndarray:
b_to_a = beta / alpha
norm = np.exp(gammaln(b_to_a) - b_to_a * np.log(b_to_a)) / alpha
t_to_t0 = taus / tau
ret_val = np.exp(-b_to_a * t_to_t0**alpha) * t_to_t0**beta
return ret_val / norm
# noinspection PyMethodOverriding
class GGAlphaEW(AbstractGG):
name = r'General \Gamma (\alpha + EW)'
parameter = [r'\tau', r'\alpha', r'\beta', r'\sigma', r'\gamma']
@staticmethod
def distribution(tau: float | np.ndarray, tau0: float,
alpha: float, beta: float, sigma: float, gamma: float) -> float | np.ndarray:
if gamma == beta or sigma == 0:
return GGAlpha.distribution(tau, tau0, alpha, beta)
b_to_a = beta / alpha
ln_b_to_a = np.log(b_to_a)
g_to_a = gamma / alpha
t_to_t0 = tau / tau0
norm = (np.exp(gammaln(b_to_a)) + sigma**(gamma-beta) * np.exp(gammaln(g_to_a) + (b_to_a-g_to_a)*ln_b_to_a)) / \
np.exp(b_to_a*ln_b_to_a) / alpha
ret_val = np.exp(-b_to_a * t_to_t0**alpha) * t_to_t0**beta * (1 + (t_to_t0*sigma)**(gamma-beta))
return ret_val / norm
# noinspection PyMethodOverriding
class GGBeta(AbstractGG):
name = r'General \Gamma (\beta)'
parameter = [r'\tau', 'a', 'b']
@staticmethod
def distribution(tau: float | np.ndarray, tau0: float, a: float, b: float) -> float | np.ndarray:
norm = a * (1+b) * np.sin(np.pi*b/(1+b)) * b**(b/(1+b)) / np.pi
ret_val = b * (tau/tau0)**a + (tau/tau0)**(-a*b)
return norm / ret_val
class GGAlphaEWBeta(Distribution):
name = r'General Gamma (WW (\alpha + EW), \beta)'
parameter = [
r'\tau_{\alpha}', r'\alpha', r'\beta', r'\sigma', r'\gamma',
r'\tau_{\beta}', 'a', 'b', 'R'
]
@staticmethod
def distribution(tau, tau0, *args):
pass
@staticmethod
def correlation(t, *args):
tau_alpha, alpha, beta, sigma, gamma, tau_beta, a, b, r = args
if r == 0:
return GGAlphaEW.correlation(t, tau_alpha, alpha, beta, sigma, gamma)
elif r == 1:
return GGBeta.correlation(t, tau_beta, a, b)
gg_alpha = GGAlphaEW.correlation(t, tau_alpha, alpha, beta, sigma, gamma)
gg_beta = GGBeta.correlation(t, tau_beta, a, b)
return ((1-r)*gg_beta + r) * gg_alpha
@staticmethod
def susceptibility(omega, *args):
t_axis, corr = GGAlphaEWBeta._correlation_for_ft(omega, *args)
ft = logft(t_axis, corr, new_x=omega, mode='complex')[1]
return 1 + omega * ft.imag + 1j * omega * ft.real
@staticmethod
def specdens(omega, *args):
t_axis, corr = GGAlphaEWBeta._correlation_for_ft(omega, *args)
return logft(t_axis, corr, new_x=omega, mode='real')[1]
@staticmethod
def _correlation_for_ft(omega, *args):
# time axis for correlation function
additional_decades = 5
lower_limit = np.log10(1 / np.max(omega)) - additional_decades
upper_limit = np.log10(1 / np.min(omega)) + additional_decades
num_steps = int(len(omega) * additional_decades)
t_axis = np.logspace(lower_limit, upper_limit, num=num_steps)
return t_axis, GGAlphaEWBeta.correlation(t_axis, *args)

168
src/nmreval/distributions/havriliaknegami.py Normal file
View File

@ -0,0 +1,168 @@
import numpy as np
from scipy.special import psi
from .base import Distribution
from ..lib.utils import ArrayLike
from ..lib.decorator import adjust_dims
from ..math.mittagleffler import mlf
from ..utils import Eu
class HavriliakNegami(Distribution):
"""
Functions based on Havriliak-Negami distribution
"""
name = 'Havriliak-Negami'
parameter = [r'\alpha', r'\gamma']
bounds = [(0, 1), (0, 1)]
@staticmethod
def distribution(tau: ArrayLike, tau0: float, alpha: float, gamma: float) -> ArrayLike:
r"""
Distribution of correlation times :math:`G(\ln\tau)`.
.. math::
G(\ln\tau) = \frac{(\tau/\tau_0)^{\alpha\gamma}
\sin\left\lbrace\gamma\arctan\left[\frac{\sin\alpha\pi}{(\tau/\tau_0)^\alpha + \cos\alpha\pi}\right] \right\rbrace}{\pi\left[1+2(\tau/\tau_0)^\alpha\cos\alpha\pi + (\tau/\tau_0)^{2\gamma}\right]}
Args:
tau (array_like) :
tau0 (array_like) :
alpha:
gamma:
Returns:
"""
if alpha == 1:
from .coledavidson import ColeDavidson
return ColeDavidson.distribution(tau, tau0, gamma)
elif gamma == 1:
from .colecole import ColeCole
return ColeCole.distribution(tau, tau0, alpha)
else:
_y = tau/tau0
om_y = (1 + 2*np.cos(np.pi*alpha)*(_y**alpha) + _y**(2*alpha))**(-gamma/2)
theta_y = np.arctan2(np.sin(np.pi * alpha), _y**alpha + np.cos(np.pi*alpha))
return np.sin(gamma*theta_y) * om_y * (_y**(alpha*gamma)) / np.pi
@staticmethod
@adjust_dims
def correlation(t: ArrayLike, tau: ArrayLike, alpha: float, gamma: float) -> ArrayLike:
r"""
Correlation function
.. math::
C(t, \tau_\text{HN}, \alpha, \gamma) = 1 - \left(\frac{t}{\tau_\text{HN}}\right)^{\alpha\gamma}
\text{E}_{\alpha,\alpha\gamma+1}^\gamma\left[-\left(\frac{t}{\tau_\text{HN}}\right)^\alpha\right]
with three-parameter Mittag-Leffler function
.. math::
\text{E}_{a,b}^g (z) = \frac{1}{\Gamma(g)} \sum_{k=0}^\infty \frac{\Gamma(g+k) z^k}{k!\Gamma(ak+b)}
Args:
t (array_like): Time axis in s.
tau (array_like): Correlation times :math:`\tau_\text{HN}` in s.
alpha (float): Cole-Cole shape parameter.
gamma (float): Cole-Davidson shape parameter.
References:
R. Garrappa: Models of dielectric relaxation based on completely monotone functions.
Frac. Calc Appl. Anal. 19 (2016) https://arxiv.org/abs/1611.04028
"""
return 1 - (t/tau)**(alpha*gamma) * mlf(-(t/tau)**alpha, alpha, alpha*gamma+1, gamma)
@staticmethod
@adjust_dims
def susceptibility(omega: ArrayLike, tau: ArrayLike, alpha: float, gamma: float) -> ArrayLike:
r"""
Complex susceptibility
.. math::
\chi(\omega, \tau, \alpha, \gamma) = \frac{1}{[1-(i\omega\tau)^\alpha]^\gamma}
Args:
omega (array_like): Frequency axis in 1/s (not Hz).
tau (array_like): Correlation times in s.
alpha (float): Shape parameter.
gamma (float): Even more shape parameter.
"""
return np.conjugate(1/(1 + (1j*omega*tau)**alpha)**gamma)
@staticmethod
@adjust_dims
def specdens(omega: ArrayLike, tau: ArrayLike, alpha: float, gamma: float) -> ArrayLike:
r"""
Spectral density :math:`J(\omega)`
.. math::
J(\omega) = \frac{\sin\left\lbrace \gamma\arctan\left[ \frac{(\omega\tau)^\alpha\sin(\alpha\pi/2)}{1+(\omega\tau)^\alpha\cos(\alpha\pi/2)} \right] \right\rbrace}
{\omega \left[ 1+(\omega\tau)^\alpha \cos(\alpha\pi/2) + (\omega\tau)^{2\alpha} \right]^{\gamma/2}}
Args:
omega (array_like): Frequency axis in 1/s (not Hz).
tau (array_like): Correlation times in s.
alpha (float): Cole-Cole shape parameter.
gamma (float): Cole-Davidson shape parameter.
Returns:
"""
omtau = (omega*tau)**alpha
zaehler = np.sin(gamma * np.arctan2(omtau*np.sin(0.5*alpha*np.pi), 1 + omtau*np.cos(0.5*alpha*np.pi)))
nenner = (1 + 2*omtau * np.cos(0.5*alpha*np.pi) + omtau**2)**(0.5*gamma)
result = (1 / omega) * (zaehler/nenner)
return result
@staticmethod
def mean(tau: ArrayLike, alpha: float, gamma: float) -> ArrayLike:
r"""
.. math::
\langle \tau \rangle \approx \tau_\text{HN} \alpha \gamma
This is an approximation given in `[1]`_
Args:
tau (array_like): Function parameter in s.
alpha (float): Shape parameter.
gamma (float): Even more shape parameter.
Reference:
_`[1]` Bauer, Th., Köhler, M., Lunkenheimer, P., Angell, C.A.:
Relaxation dynamics and ionic conductivity in a fragile plastic crystal.
J. Chem. Phys. 133, 144509 (2010). https://doi.org/10.1063/1.3487521
"""
return tau * alpha*gamma
@staticmethod
def logmean(tau: ArrayLike, alpha: float, gamma: float) -> ArrayLike:
r"""
Calculate mean logarithm of tau
.. math::
\langle \ln \tau \rangle = \ln \tau + \frac{\Psi(\gamma) + \text{Eu}}{\alpha}
Args:
tau:
alpha (float):
gamma (float):
"""
return np.log(tau) + (psi(gamma)+Eu)/alpha
@staticmethod
def max(tau: ArrayLike, alpha: float, gamma: float) -> ArrayLike:
return tau*(np.sin(0.5*np.pi*alpha*gamma/(gamma+1)) / np.sin(0.5*np.pi*alpha/(gamma+1)))**(1/alpha)

84
src/nmreval/distributions/intermolecular.py Normal file
View File

@ -0,0 +1,84 @@
import numpy as np
from scipy.integrate import quad
from .base import Distribution
class FFHS(Distribution):
name = 'Intermolecular (FFHS)'
parameter = None
@staticmethod
def distribution(taus, tau0, *args):
# Distribution over tau, not log(tau), like in other cases
z = taus / tau0
return 54 * np.sqrt(z) / (162*z**3 + 18*z**2 - 4*z + 2) / np.pi / tau0
@staticmethod
def correlation(t, tau0, *args):
def integrand(u, tt, tau0):
return FFHS.distribution(u, tau0) * np.exp(-tt/u) / u
ret_val = np.array([quad(integrand, 0, np.infty, args=(tt, tau0), epsabs=1e-12, epsrel=1e-12)[0] for tt in t])
return ret_val
@staticmethod
def specdens(omega, tau0, *args):
def integrand(u, o, tau0):
return u**4 * tau0 / (81 + 9*u**2 - 2*u**4 + u**6) / (u**4 + (o*tau0)**2)
# return FFHS.distribution(u, tau0) * u / (1+o**2 * u**2)
ret_val = np.array([quad(integrand, 0, np.infty, args=(o, tau0), epsabs=1e-12, epsrel=1e-12)[0] for o in omega])
return ret_val
@staticmethod
def susceptibility(omega, tau0, *args):
def integrand_real(u, o, tau0):
return FFHS.distribution(u, tau0) / (1+o**2 * u**2)
def integrand_imag(u, o, tau0):
return FFHS.distribution(u, tau0) * o*u / (1+o**2 * u**2)
ret_val = np.array([quad(integrand_real, 0, np.infty, args=(o, tau0),
epsabs=1e-12, epsrel=1e-12)[0] for o in omega], dtype=complex)
ret_val.imag += np.array([quad(integrand_imag, 0, np.infty, args=(o, tau0),
epsabs=1e-12, epsrel=1e-12)[0] for o in omega])
return ret_val
# class Bessel(Distribution):
# name = 'Intermolecular (Bessel)'
# parameter = None
#
# @staticmethod
# def distribution(taus, tau0, *args):
# x = np.sqrt(tau0 / taus)
# return special.jv(1.5, x)**2 / tau0 / 2
#
# @staticmethod
# def correlation(t, tau0, *args):
# def integrand(lx, tt):
# x = np.exp(lx)
# return Bessel.distribution(x, tau0)*np.exp(-tt/lx)
#
# logaxis = np.linspace(-20+np.log(tau0), 20+np.log(tau0), num=5001)
#
# ret_val = np.array([simps(integrand(logaxis, tt), logaxis) for tt in t])
#
# return ret_val
#
# @staticmethod
# def susceptibility(omega, tau0, *args):
# def integrand(lx, o):
# x = np.exp(lx)
# return Bessel.distribution(x, tau0) * 1/(1+1j*omega*x)
#
# logaxis = np.linspace(-20 + np.log(tau0), 20 + np.log(tau0), num=5001)
#
# ret_val = np.array([simps(integrand(logaxis, o), logaxis) for o in omega])
#
# return ret_val

63
src/nmreval/distributions/kww.py Normal file
View File

@ -0,0 +1,63 @@
import numpy as np
from scipy.interpolate import interp1d
from scipy.special import gamma
from .base import Distribution
from ..lib.decorator import adjust_dims
from ..math.kww import kww_cos, kww_sin
from ..utils.constants import Eu
class KWW(Distribution):
name = 'KWW'
parameter = [r'\beta']
boounds = [(0, 1)]
@staticmethod
def distribution(taus, tau, beta):
if not (0.1 <= beta <= 0.9):
raise ValueError('KWW distribution is only defined between 0.1 and 0.9')
b_supp = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]
B_supp = [0.145, 0.197, 0.243, 0.285, 0.382, 0.306, 0.360, 0.435, 0.700]
C_supp = [0.89, 0.50, 0.35, 0.25, 0, 0.13, 0.22, 0.4015, 0.32]
B = interp1d(b_supp, B_supp)(beta)
C = interp1d(b_supp, C_supp)(beta)
tt = tau/taus
delta = beta * abs(beta-0.5) / (1-beta)
if beta > 0.5:
f = 1 + C * tt**delta
else:
f = 1 / (1 + C * tt**delta)
ret_val = tau * B * np.exp(-(1-beta) * beta**(beta/(1-beta)) / tt**(beta/(1-beta))) / tt**((1-0.5*beta)/(1-beta))
return ret_val * f / taus
@staticmethod
@adjust_dims
def correlation(t, tau, beta):
return np.exp(-(t/tau)**beta)
@staticmethod
@adjust_dims
def susceptibility(omega, tau, beta):
return 1-omega*kww_sin(omega, tau, beta) + 1j*omega*kww_cos(omega, tau, beta)
@staticmethod
@adjust_dims
def specdens(omega, tau, beta):
return kww_cos(omega, tau, beta)
@staticmethod
def mean(tau, beta):
return tau/beta * gamma(1 / beta)
@staticmethod
def logmean(tau, beta):
return (1-1/beta) * Eu + np.log(tau)
@staticmethod
def max(tau, beta):
return (1.7851 - 0.87052*beta - 0.028836*beta**2 + 0.11391*beta**3) * tau

121
src/nmreval/distributions/loggaussian.py Normal file
View File

@ -0,0 +1,121 @@
from multiprocessing import Pool, cpu_count
from itertools import product
import numpy as np
try:
from scipy.integrate import simpson
except ImportError:
from scipy.integrate import simps as simpson
from scipy.integrate import quad
from .base import Distribution
__all__ = ['LogGaussian']
class LogGaussian(Distribution):
name = 'Log-Gaussian'
parameter = [r'\sigma']
bounds = [(0, 10)]
@staticmethod
def distribution(tau, tau0, sigma: float):
return np.exp(-0.5*(np.log(tau/tau0)/sigma)**2)/np.sqrt(2*np.pi)/sigma
@staticmethod
def correlation(t, tau0, sigma: float):
_t = np.atleast_1d(t)
_tau = np.atleast_1d(tau0)
pool = Pool(processes=min(cpu_count(), 4))
integration_ranges = [(omega_i, tau_j, sigma) for (omega_i, tau_j) in product(_t, _tau)]
with np.errstate(divide='ignore'):
res = np.array(pool.map(_integrate_process_time, integration_ranges))
ret_val = res.reshape((_t.shape[0], _tau.shape[0]))
return ret_val.squeeze()
@staticmethod
def susceptibility(omega, tau0, sigma: float):
_omega = np.atleast_1d(omega)
_tau = np.atleast_1d(tau0)
pool = Pool(processes=min(cpu_count(), 4))
integration_ranges = [(omega_i, tau_j, sigma) for (omega_i, tau_j) in product(_omega, _tau)]
with np.errstate(divide='ignore'):
res_real = np.array(pool.map(_integrate_process_imag, integration_ranges))
res_imag = np.array(pool.map(_integrate_process_real, integration_ranges))
ret_val = (res_real+1j*res_imag).reshape((_omega.shape[0], _tau.shape[0]))
return ret_val.squeeze()
@staticmethod
def specdens(omega, tau0, sigma):
_omega = np.atleast_1d(omega)
_tau = np.atleast_1d(tau0)
pool = Pool(processes=min(cpu_count(), 4))
integration_ranges = [(omega_i, tau_j, sigma) for (omega_i, tau_j) in product(_omega, _tau)]
with np.errstate(divide='ignore'):
res = np.array(pool.map(_integrate_process_imag, integration_ranges))
ret_val = res.reshape((_omega.shape[0], _tau.shape[0]))
ret_val /= _omega[:, None]
return ret_val.squeeze()
def mean(*args):
return args[0]*np.exp(args[1]**2 / 2)
def _integrate_process_imag(args):
omega_i, tau_j, sigma = args
area = quad(_integrand_freq_imag_high, 0, 50, args=(omega_i, tau_j, sigma))[0]
area += quad(_integrand_freq_imag_low, -50, 0, args=(omega_i, tau_j, sigma))[0]
return area
def _integrate_process_real(args):
omega_i, tau_j, sigma = args
area = quad(_integrand_freq_real_high, 0, 50, args=(omega_i, tau_j, sigma))[0]
area += quad(_integrand_freq_real_low, -50, 0, args=(omega_i, tau_j, sigma))[0]
return area
def _integrate_process_time(args):
omega_i, tau_j, sigma = args
return quad(_integrand_time, -50, 50, args=(omega_i, tau_j, sigma))[0]
def _integrand_time(u, t, tau, sigma):
uu = np.exp(u)
return LogGaussian.distribution(uu, tau, sigma) * np.exp(-t/uu)
# integrands
def _integrand_freq_imag_low(u, omega, tau, sigma):
# integrand
uu = np.exp(u)
return LogGaussian.distribution(uu, tau, sigma) * omega * uu / (1 + (omega*uu)**2)
def _integrand_freq_imag_high(u, omega, tau, sigma):
uu = np.exp(-u)
return LogGaussian.distribution(1/uu, tau, sigma) * omega * uu / (uu**2 + omega**2)
def _integrand_freq_real_low(u, omega, tau, sigma):
uu = np.exp(u)
return LogGaussian.distribution(uu, tau, sigma) / (1 + (omega*uu)**2)
def _integrand_freq_real_high(u, omega, tau, sigma):
uu = np.exp(-2*u)
return LogGaussian.distribution(np.exp(u), tau, sigma) * uu / (uu + omega**2)