Compiled with numpy 2.0

Makefile

@@ -31,7 +31,7 @@ test:
 
 types:
 	${python} -m monkeytype run $$(which ${pytest}) ./tests
-	${python} -m monkeytype list-modules | grep ${pkg} | parallel -j${j} "${python} -m monkeytype apply {}"
+	${python} -m monkeytype list-modules | grep ${pkg} | parallel -j${j} "${python} -m monkeytype apply {} > /dev/null && echo {}"
 
 cov:
 	${python} -m pytest -vrs --cov=${pkg} --cov-report html tests

README.rst

@@ -23,11 +23,18 @@ Breaking changes from ``qml``:
 
 * FCHL representations callable interface to be consistent with other representations (e.i. atoms, coordinates)
 
-
 ==============
 How to install
 ==============
 
+You need a fortran compiler and math library. Default is `gfortran` and `openblas`.
+
+
+.. code-block:: bash
+
+    sudo apt install libopenblas-dev gcc
+
+
 A proper pip-package is on the way, for now
 
 .. code-block:: bash

_compile.py

@@ -5,67 +5,85 @@
 import sys
 from pathlib import Path
 
+import numpy as np
+
 f90_modules = {
     "representations/frepresentations": ["frepresentations.f90"],
     "representations/facsf": ["facsf.f90"],
     "representations/fslatm": ["fslatm.f90"],
     "representations/arad/farad_kernels": ["farad_kernels.f90"],
     "representations/fchl/ffchl_module": [
-        "ffchl_module.f90",
-        "ffchl_scalar_kernels.f90",
         "ffchl_kernel_types.f90",
+        "ffchl_module.f90",
+        "ffchl_module_ef.f90",
         "ffchl_kernels.f90",
-        "ffchl_electric_field_kernels.f90",
+        "ffchl_scalar_kernels.f90",
+        "ffchl_kernels_ef.f90",
         "ffchl_force_kernels.f90",
     ],
     "solvers/fsolvers": ["fsolvers.f90"],
     "kernels/fdistance": ["fdistance.f90"],
-    "kernels/fkernels": ["fkernels.f90", "fkpca.f90"],
+    "kernels/fkernels": [
+        "fkernels.f90",
+        "fkpca.f90",
+        "fkwasserstein.f90",
+    ],
     "kernels/fgradient_kernels": ["fgradient_kernels.f90"],
+    "utils/fsettings": ["fsettings.f90"],
 }
 
 
 def find_mkl():
-
-    return
+    raise NotImplementedError()
 
 
-def find_flags(fcc: str):
+def find_env() -> dict[str, str]:
     """Find compiler flags"""
 
+    # TODO Find ifort flags, choose from FCC
     # TODO Find math lib
     # TODO Find os
 
-    # -lgomp", "-lpthread", "-lm", "-ldl
-    # ["-L${MKLROOT}/lib/intel64", "-lmkl_rt"]
-
-    # COMPILER_FLAGS = ["-O3", "-fopenmp", "-m64", "-march=native", "-fPIC",
-    #                 "-Wno-maybe-uninitialized", "-Wno-unused-function", "-Wno-cpp"]
-    # LINKER_FLAGS = ["-lgomp"]
+    COMPILER_FLAGS = [
+        "-O3",
+        "-fopenmp",
+        "-m64",
+        "-march=native",
+        "-fPIC",
+        "-Wno-maybe-uninitialized",
+        "-Wno-unused-function",
+        "-Wno-cpp",
+    ]
 
     extra_flags = ["-lgomp", "-lpthread", "-lm", "-ldl"]
+    math_flags = ["-L/usr/lib/", "-lblas", "-llapack"]
 
-    flags = ["-L/usr/lib/", "-lblas", "-llapack"] + extra_flags
-
-    return flags
-
-
-def find_fcc():
-    """Find the fortran compiler. Either gnu or intel"""
-
-    # fcc = "ifort"
+    fflags = [] + COMPILER_FLAGS
+    ldflags = [] + extra_flags + math_flags
     fcc = "gfortran"
 
-    return fcc
+    env = {"FFLAGS": " ".join(fflags), "LDFLAGS": " ".join(ldflags), "FCC": fcc}
+
+    return env
 
 
 def main():
     """Compile f90 in src/qmllib"""
 
-    fcc = find_fcc()
-    flags = find_flags(fcc)
+    print(f"Using numpy {np.__version__}")
+
+    # Find and set Fortran compiler, compiler flags and linker flags
+    env = find_env()
+    for key, value in env.items():
+        print(f"export {key}='{value}'")
+        os.environ[key] = value
+
+    f2py = [sys.executable, "-m", "numpy.f2py"]
 
-    os.environ["FCC"] = fcc
+    meson_flags = [
+        "--backend",
+        "meson",
+    ]
 
     for module_name, module_sources in f90_modules.items():
 
@@ -74,9 +92,7 @@ def main():
         stem = path.stem
 
         cwd = Path("src/qmllib") / parent
-        cmd = (
-            [sys.executable, "-m", "numpy.f2py", "-c"] + flags + module_sources + ["-m", str(stem)]
-        )
+        cmd = f2py + ["-c"] + module_sources + ["-m", str(stem)] + meson_flags
         print(cwd, " ".join(cmd))
 
         proc = subprocess.run(cmd, cwd=cwd, capture_output=True, text=True)
@@ -85,9 +101,9 @@ def main():
         exitcode = proc.returncode
 
         if exitcode > 0:
-            print(stdout)
-            print()
             print(stderr)
+            print()
+            print(stdout)
             exit(exitcode)
 
 

pyproject.toml

@@ -1,5 +1,5 @@
 [build-system]
-requires = ["setuptools", "numpy", "meson", "ninja"]
+requires = ["setuptools", "numpy>2.0", "meson", "ninja"]
 build-backend = "setuptools.build_meta"
 
 [project]
@@ -19,7 +19,7 @@ classifiers = [
     "Topic :: Scientific/Engineering :: Chemistry",
 ]
 keywords = ["qml", "quantum chemistry", "machine learning"]
-dependencies=["numpy", "scipy"]
+dependencies=["numpy>2.0", "scipy"]
 
 [project.urls]
 Homepage = "https://qmlcode.org"

src/qmllib/kernels/fkernels.f90

@@ -1,83 +1,4 @@
 
-module searchtools
-
-   implicit none
-
-contains
-
-   function searchsorted(all_values, sorted) result(cdf_idx)
-
-      implicit none
-
-      double precision, dimension(:), intent(in) :: all_values
-      double precision, dimension(:), intent(in) :: sorted
-
-      integer, allocatable, dimension(:) :: cdf_idx
-
-      double precision :: val
-
-      integer :: i, j, n, m
-
-      n = size(all_values) - 1
-      m = size(sorted)
-
-      allocate (cdf_idx(n))
-
-      cdf_idx(:) = 0
-
-      do i = 1, n
-
-         val = all_values(i)
-
-         do j = 1, m
-
-            !write (*,*) i, j, sorted(j), val
-
-            ! if ((sorted(j) <= val) .and. (val < sorted(j+1)))  then
-            if (sorted(j) > val) then
-
-               cdf_idx(i) = j - 1
-               !write(*,*) "found"
-               exit
-
-               !  endif
-            else !iif (val > maxval(sorted)) then
-               cdf_idx(i) = m
-            end if
-
-         end do
-
-      end do
-
-   end function searchsorted
-
-   recursive subroutine quicksort(a, first, last)
-      implicit none
-      double precision ::  a(*), x, t
-      integer first, last
-      integer i, j
-
-      x = a((first + last)/2)
-      i = first
-      j = last
-      do
-         do while (a(i) < x)
-            i = i + 1
-         end do
-         do while (x < a(j))
-            j = j - 1
-         end do
-         if (i >= j) exit
-         t = a(i); a(i) = a(j); a(j) = t
-         i = i + 1
-         j = j - 1
-      end do
-      if (first < i - 1) call quicksort(a, first, i - 1)
-      if (j + 1 < last) call quicksort(a, j + 1, last)
-   end subroutine quicksort
-
-end module searchtools
-
 subroutine fget_local_kernels_gaussian(q1, q2, n1, n2, sigmas, &
         & nm1, nm2, nsigmas, kernels)
 
@@ -796,95 +717,3 @@ subroutine fsargan_kernel(a, na, b, nb, k, sigma, gammas, ng)
    deallocate (prefactor)
 
 end subroutine fsargan_kernel
-
-subroutine fwasserstein_kernel(a, na, b, nb, k, sigma, p, q)
-
-   use searchtools
-   implicit none
-
-   double precision, dimension(:, :), intent(in) :: a
-   double precision, dimension(:, :), intent(in) :: b
-
-   double precision, allocatable, dimension(:, :) :: asorted
-   double precision, allocatable, dimension(:, :) :: bsorted
-
-   double precision, allocatable, dimension(:) :: rep
-
-   integer, intent(in) :: na, nb
-
-   double precision, dimension(:, :), intent(inout) :: k
-   double precision, intent(in) :: sigma
-
-   integer, intent(in) :: p
-   integer, intent(in) :: q
-
-   double precision :: inv_sigma
-
-   integer :: i, j, l
-   integer :: rep_size
-
-   double precision, allocatable, dimension(:) :: deltas
-   double precision, allocatable, dimension(:) :: all_values
-
-   double precision, allocatable, dimension(:) :: a_cdf
-   double precision, allocatable, dimension(:) :: b_cdf
-   integer, allocatable, dimension(:) :: a_cdf_idx
-   integer, allocatable, dimension(:) :: b_cdf_idx
-
-   rep_size = size(a, dim=1)
-   allocate (asorted(rep_size, na))
-   allocate (bsorted(rep_size, nb))
-   allocate (rep(rep_size))
-
-   allocate (all_values(rep_size*2))
-   allocate (deltas(rep_size*2 - 1))
-
-   allocate (a_cdf(rep_size*2 - 1))
-   allocate (b_cdf(rep_size*2 - 1))
-
-   allocate (a_cdf_idx(rep_size*2 - 1))
-   allocate (b_cdf_idx(rep_size*2 - 1))
-
-   asorted(:, :) = a(:, :)
-   bsorted(:, :) = b(:, :)
-
-   do i = 1, na
-      rep(:) = asorted(:, i)
-      call quicksort(rep, 1, rep_size)
-      asorted(:, i) = rep(:)
-   end do
-
-   do i = 1, nb
-      rep(:) = bsorted(:, i)
-      call quicksort(rep, 1, rep_size)
-      bsorted(:, i) = rep(:)
-   end do
-
-   !$OMP PARALLEL DO PRIVATE(all_values,a_cdf_idx,b_cdf_idx,a_cdf,b_cdf,deltas)
-   do j = 1, nb
-      do i = 1, na
-
-         all_values(:rep_size) = asorted(:, i)
-         all_values(rep_size + 1:) = bsorted(:, j)
-
-         call quicksort(all_values, 1, 2*rep_size)
-
-         do l = 1, 2*rep_size - 1
-            deltas(l) = all_values(l + 1) - all_values(l)
-         end do
-
-         a_cdf_idx = searchsorted(all_values, asorted(:, i))
-         b_cdf_idx = searchsorted(all_values, bsorted(:, j))
-
-         a_cdf(:) = a_cdf_idx(:)
-         b_cdf(:) = b_cdf_idx(:)
-         a_cdf(:) = a_cdf(:)/rep_size
-         b_cdf(:) = b_cdf(:)/rep_size
-
-         ! k(i,j) = exp(-sum(abs(a_cdf-b_cdf)*deltas)/sigma)
-         k(i, j) = exp(-(sum((abs(a_cdf - b_cdf)**p)*deltas)**(1.0d0/p))**q/sigma)
-
-      end do
-   end do
-   !$OMP END PARALLEL DO
-end subroutine fwasserstein_kernel