From 2042148d0f80d466880223d96b66e2da48931db0 Mon Sep 17 00:00:00 2001 From: Dominik Demuth Date: Mon, 17 Apr 2023 20:17:54 +0200 Subject: [PATCH] C function for energy distribution spectral density --- src/nmreval/clib/integrate.c | 27 ++++++++++++++++- src/nmreval/distributions/energy.py | 45 ++++++++++++++++++++++------- 2 files changed, 60 insertions(+), 12 deletions(-) diff --git a/src/nmreval/clib/integrate.c b/src/nmreval/clib/integrate.c index 4654530..1d221cf 100644 --- a/src/nmreval/clib/integrate.c +++ b/src/nmreval/clib/integrate.c @@ -1,6 +1,8 @@ /* integrands used in quadrature integration with scipy's LowLevelCallables */ #include +#define KB 8.617333262145179e-05 + /* FFHS functions */ double ffhsSD(double x, void *user_data) { double *c = (double *)user_data; @@ -18,7 +20,7 @@ double ffhsSD(double x, void *user_data) { /* log-gaussian functions */ double logNormalDist(double tau, double tau0, double sigma) { - return exp(- pow((log(tau/tau0) / sigma), 2) / 2) / sqrt(2*M_PI)/sigma; + return exp(- pow((log(tau/tau0) / sigma), 2) / 2.) / sqrt(2*M_PI)/sigma; } double logGaussianSD_high(double u, void *user_data) { @@ -61,3 +63,26 @@ double logGaussianCorrelation(double x, void *user_data) { return dist * exp(-t/uu); } + + +double normalDist(double x, double x0, double sigma) { + return exp(- pow((x-x0) / sigma, 2) / 2.) / sqrt(2 * M_PI) / sigma; +} + +double rate(double tau0, double ea, double t) { + return exp(-ea / t / KB) / tau0; +} + +double energyDist_SD(double x, void *user_data) { + double *c = (double *)user_data; + + double omega = c[0]; + double tau0 = c[1]; + double e_m = c[2]; + double e_b = c[3]; + double temp = c[4]; + + double r = rate(tau0, x, temp); + + return r/(pow(r, 2) + pow(omega, 2)) * normalDist(x, e_m, e_b); +} diff --git a/src/nmreval/distributions/energy.py b/src/nmreval/distributions/energy.py index a97a462..3ef11df 100644 --- a/src/nmreval/distributions/energy.py +++ b/src/nmreval/distributions/energy.py @@ -1,13 +1,18 @@ +import ctypes from itertools import product import numpy as np +from scipy import LowLevelCallable from scipy.integrate import quad, simps as simpson from .base import Distribution from ..lib.utils import ArrayLike from ..utils.constants import kB +from .helper import HAS_C_FUNCS, lib + +# noinspection PyMethodOverriding class EnergyBarriers(Distribution): name = 'Energy barriers' parameter = [r'\tau_{0}', r'E_{m}', r'\Delta E'] @@ -51,7 +56,7 @@ class EnergyBarriers(Distribution): 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) + 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) - @@ -60,23 +65,41 @@ class EnergyBarriers(Distribution): return np.array(ret_val) @staticmethod - def specdens(omega, temperature, *args): + def specdens(omega, temperature, tau0: float, e_m: float, e_b: float): # 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) + if HAS_C_FUNCS: + ret_val = EnergyBarriers.spec_dens_c(omega, temperature, tau0, e_m, e_b) + else: + ret_val = EnergyBarriers.spec_dens_py(omega, temperature, tau0, e_m, e_b) - 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.squeeze() - return ret_val + @staticmethod + def spec_dens_c(omega: np.ndarray, temperature: np.ndarray, tau0: float, e_m: float, e_b: float) -> np.ndarray: + res = [] + for o, t in product(omega, temperature): + c = (ctypes.c_double * 5)(o, tau0, e_m, e_b, t) + user_data = ctypes.cast(ctypes.pointer(c), ctypes.c_void_p) + area = quad(LowLevelCallable(lib.energyDist_SD, user_data), 0, np.infty, epsabs=1e-10)[0] + res.append(area) + + return np.array(res) + + @staticmethod + def spec_dens_py(omega: np.ndarray, temperature: np.ndarray, tau0: float, e_m: float, e_b: float) -> np.ndarray: + 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) + + e_axis = np.linspace(max(0., e_m-50*e_b), e_m+50*e_b, num=5001) + + ret_val = [simpson(integrand(e_axis, o, tau0, e_m, e_b, tt), e_axis) for o in omega for tt in temperature] + + return np.array(ret_val) @staticmethod def mean(*args):