add C function for FFHS integration

src/nmreval/distributions/intermolecular.py

@@ -1,8 +1,24 @@
+import ctypes
+import os.path
+from pathlib import Path
+
 import numpy as np
+from scipy import LowLevelCallable
 from scipy.integrate import quad
 
 from .base import Distribution
 
+try:
+    print(str(Path(__file__).parents[1] / 'clib' / 'integrate.so'))
+    lib = ctypes.CDLL(str(Path(__file__).parents[1] / 'clib' / 'integrate.so'))
+    lib.ffhs_sd.restype = ctypes.c_double
+    lib.ffhs_sd.argtypes = (ctypes.c_double, ctypes.c_void_p)
+    HAS_C_FUNCS = True
+    print('Use C functions')
+except OSError:
+    HAS_C_FUNCS = False
+    print('Use python functions')
+
 
 class FFHS(Distribution):
     name = 'Intermolecular (FFHS)'
@@ -24,7 +40,7 @@ class FFHS(Distribution):
         return ret_val
 
     @staticmethod
-    def specdens(omega, tau0, *args):
+    def specdens_py(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)
@@ -33,6 +49,17 @@ class FFHS(Distribution):
 
         return ret_val * 54 / np.pi
 
+    @staticmethod
+    def specdens_c(omega, tau0, *args):
+        res = []
+        for o in omega:
+            c = (ctypes.c_double * 2)(o, tau0)
+            user_data = ctypes.cast(ctypes.pointer(c), ctypes.c_void_p)
+            func = LowLevelCallable(lib.ffhs_sd, user_data)
+            res.append(quad(func, 0, np.infty)[0])
+
+        return np.array(res) * 54 / np.pi
+
     @staticmethod
     def susceptibility(omega, tau0, *args):
         def integrand_real(u, o, tau0):
@@ -48,6 +75,9 @@ class FFHS(Distribution):
         return ret_val
 
 
+FFHS.specdens = FFHS.specdens_c if HAS_C_FUNCS else FFHS.specdens_py
+
+
 # class Bessel(Distribution):
 #     name = 'Intermolecular (Bessel)'
 #     parameter = None