lower epsabs for integration of energy barrier spec dens

This commit is contained in:
Dominik Demuth 2023-04-19 18:24:29 +02:00
parent 2042148d0f
commit 2d472bd44e
4 changed files with 18 additions and 16 deletions

View File

@ -23,7 +23,7 @@ LDFLAGS = -shared
C_DIR = src/nmreval/clib C_DIR = src/nmreval/clib
all : ui all : ui compile
ui : $(PYQT_UI) ui : $(PYQT_UI)

View File

@ -1,5 +1,5 @@
import ctypes
from itertools import product from itertools import product
from ctypes import c_double, cast, pointer, c_void_p
import numpy as np import numpy as np
from scipy import LowLevelCallable from scipy import LowLevelCallable
@ -82,9 +82,9 @@ class EnergyBarriers(Distribution):
def spec_dens_c(omega: np.ndarray, temperature: np.ndarray, tau0: float, e_m: float, e_b: float) -> np.ndarray: def spec_dens_c(omega: np.ndarray, temperature: np.ndarray, tau0: float, e_m: float, e_b: float) -> np.ndarray:
res = [] res = []
for o, t in product(omega, temperature): for o, t in product(omega, temperature):
c = (ctypes.c_double * 5)(o, tau0, e_m, e_b, t) c = (c_double * 5)(o, tau0, e_m, e_b, t)
user_data = ctypes.cast(ctypes.pointer(c), ctypes.c_void_p) user_data = cast(pointer(c), c_void_p)
area = quad(LowLevelCallable(lib.energyDist_SD, user_data), 0, np.infty, epsabs=1e-10)[0] area = quad(LowLevelCallable(lib.energyDist_SD, user_data), 0, np.infty, epsabs=1e-13)[0]
res.append(area) res.append(area)
return np.array(res) return np.array(res)

View File

@ -1,24 +1,26 @@
from pathlib import Path from pathlib import Path
import ctypes from ctypes import CDLL, c_double, c_void_p
from ..lib.logger import logger from ..lib.logger import logger
lib = None lib = None
try: try:
lib = ctypes.CDLL(str(Path(__file__).parents[1] / 'clib' / 'integrate.so')) lib = CDLL(str(Path(__file__).parents[1] / 'clib' / 'integrate.so'))
# FFHS integrand # FFHS integrand
lib.ffhsSD.restype = ctypes.c_double lib.ffhsSD.restype = c_double
lib.ffhsSD.argtypes = (ctypes.c_double, ctypes.c_void_p) lib.ffhsSD.argtypes = (c_double, c_void_p)
# Log-Gaussian integrands # Log-Gaussian integrands
lib.logGaussianSD_high.restype = ctypes.c_double lib.logGaussianSD_high.restype = c_double
lib.logGaussianSD_high.argtypes = (ctypes.c_double, ctypes.c_void_p) lib.logGaussianSD_high.argtypes = (c_double, c_void_p)
lib.logGaussianSD_low.restype = ctypes.c_double lib.logGaussianSD_low.restype = c_double
lib.logGaussianSD_low.argtypes = (ctypes.c_double, ctypes.c_void_p) lib.logGaussianSD_low.argtypes = (c_double, c_void_p)
lib.energyDist_SD.restype = c_double
lib.energyDist_SD.argtypes = (c_double, c_void_p)
HAS_C_FUNCS = True HAS_C_FUNCS = True
logger.info('Use C functions') logger.info('Use C functions')

View File

@ -86,7 +86,7 @@ class EnergyFC(_AbstractFC):
name = 'Energy distribution' name = 'Energy distribution'
params = ['C', 'T'] + EnergyBarriers.parameter params = ['C', 'T'] + EnergyBarriers.parameter
bounds = [(0, None), (0, None), (0, None), (0, None)] bounds = [(0, None), (0, None), (0, None), (0, None)]
ralax = Relaxation(distribution=EnergyBarriers) relax = Relaxation(distribution=EnergyBarriers)
class _AbstractFCDipolar(_AbstractFC): class _AbstractFCDipolar(_AbstractFC):