Fix wrong gradient at nodes of AOs

CHange version
Removed scorecards
2024-06-01 02:45:43 +02:00 · 2024-05-08 13:56:32 +02:00 · 2024-04-10 16:38:10 +02:00 · 2024-03-29 09:27:40 +01:00 · 2024-02-26 09:23:23 +01:00 · 2024-02-24 01:18:06 +02:00
37 changed files with 19225 additions and 13848 deletions
--- a/.github/workflows/devskim.yml
+++ b/.github/workflows/devskim.yml
@ -1,37 +0,0 @@
-# This workflow uses actions that are not certified by GitHub.
-# They are provided by a third-party and are governed by
-# separate terms of service, privacy policy, and support
-# documentation.
-
-name: DevSkim
-
-on:
-  push:
-    branches: [ "master" ]
-  pull_request:
-    branches: [ "master" ]
-  schedule:
-    - cron: '19 5 * * 2'
-
-permissions:
-  contents: read
-
-jobs:
-  lint:
-    name: DevSkim
-    runs-on: ubuntu-20.04
-    permissions:
-      actions: read
-      contents: read
-      security-events: write
-    steps:
-      - name: Checkout code
-        uses: actions/checkout@93ea575cb5d8a053eaa0ac8fa3b40d7e05a33cc8
-
-      - name: Run DevSkim scanner
-        uses: microsoft/DevSkim-Action@a8a9e06bab570db990fe7351ae9d4d444b9489ca
-
-      - name: Upload DevSkim scan results to GitHub Security tab
-        uses: github/codeql-action/upload-sarif@678fc3afe258fb2e0cdc165ccf77b85719de7b3c
-        with:
-          sarif_file: devskim-results.sarif
--- a/.github/workflows/scorecards.yml
+++ b/.github/workflows/scorecards.yml
@ -1,72 +0,0 @@
-# This workflow uses actions that are not certified by GitHub. They are provided
-# by a third-party and are governed by separate terms of service, privacy
-# policy, and support documentation.
-
-name: Scorecards supply-chain security
-on:
-  # For Branch-Protection check. Only the default branch is supported. See
-  # https://github.com/ossf/scorecard/blob/main/docs/checks.md#branch-protection
-  branch_protection_rule:
-  # To guarantee Maintained check is occasionally updated. See
-  # https://github.com/ossf/scorecard/blob/main/docs/checks.md#maintained
-  schedule:
-    - cron: '33 0 * * 5'
-  push:
-    branches: [ "master" ]
-
-# Declare default permissions as read only.
-permissions: read-all
-
-jobs:
-  analysis:
-    name: Scorecards analysis
-    runs-on: ubuntu-latest
-    permissions:
-      # Needed to upload the results to code-scanning dashboard.
-      security-events: write
-      # Needed to publish results and get a badge (see publish_results below).
-      id-token: write
-      # Uncomment the permissions below if installing in a private repository.
-      # contents: read
-      # actions: read
-
-    steps:
-      - name: "Checkout code"
-        uses: actions/checkout@93ea575cb5d8a053eaa0ac8fa3b40d7e05a33cc8 # v3.1.0
-        with:
-          persist-credentials: false
-
-      - name: "Run analysis"
-        uses: ossf/scorecard-action@99c53751e09b9529366343771cc321ec74e9bd3d # v2.0.6
-        with:
-          results_file: results.sarif
-          results_format: sarif
-          # (Optional) Read-only PAT token. Uncomment the `repo_token` line below if:
-          # - you want to enable the Branch-Protection check on a *public* repository, or
-          # - you are installing Scorecards on a *private* repository
-          # To create the PAT, follow the steps in https://github.com/ossf/scorecard-action#authentication-with-pat.
-          # repo_token: ${{ secrets.SCORECARD_READ_TOKEN }}
-
-          # Public repositories:
-          #   - Publish results to OpenSSF REST API for easy access by consumers
-          #   - Allows the repository to include the Scorecard badge.
-          #   - See https://github.com/ossf/scorecard-action#publishing-results.
-          # For private repositories:
-          #   - `publish_results` will always be set to `false`, regardless
-          #     of the value entered here.
-          publish_results: true
-
-      # Upload the results as artifacts (optional). Commenting out will disable uploads of run results in SARIF
-      # format to the repository Actions tab.
-      - name: "Upload artifact"
-        uses: actions/upload-artifact@3cea5372237819ed00197afe530f5a7ea3e805c8 # v3.1.0
-        with:
-          name: SARIF file
-          path: results.sarif
-          retention-days: 5
-
-      # Upload the results to GitHub's code scanning dashboard.
-      - name: "Upload to code-scanning"
-        uses: github/codeql-action/upload-sarif@807578363a7869ca324a79039e6db9c843e0e100 # v2.1.27
-        with:
-          sarif_file: results.sarif
--- a/.github/workflows/test-build.yml
+++ b/.github/workflows/test-build.yml
@ -147,60 +147,48 @@ jobs:
        name: test-report-ubuntu-debug
        path: _build_hpc/test-suite.log

-#  x86_macos:
-#
-#    runs-on: macos-latest
-#    name: x86 MacOS latest
-#
-#    steps:
-#    - uses: actions/checkout@e2f20e631ae6d7dd3b768f56a5d2af784dd54791
-#    - name: install dependencies
-#      run: brew install emacs hdf5 automake pkg-config
-#
-#    - name: Symlink gfortran (macOS)
-#      if: runner.os == 'macOS'
-#      run: |
-#        # make sure gfortran is available
-#        # https://github.com/actions/virtual-environments/issues/2524
-#        # https://github.com/cbg-ethz/dce/blob/master/.github/workflows/pkgdown.yaml
-#        sudo ln -s /usr/local/bin/gfortran-10 /usr/local/bin/gfortran
-#        sudo mkdir /usr/local/gfortran
-#        sudo ln -s /usr/local/Cellar/gcc@10/*/lib/gcc/10 /usr/local/gfortran/lib
-#        gfortran --version
-#
-#    - name: Install the latest TREXIO from the GitHub clone
-#      run: |
-#        git clone https://github.com/TREX-CoE/trexio.git
-#        cd trexio
-#        ./autogen.sh
-#        ./configure --prefix=${PWD}/_install --enable-silent-rules
-#        make -j 4
-#        make install
-#
-#    - name: Test TREXIO
-#      run: make -j 4 check
-#      working-directory: trexio
-#
-#    - name: Archive TREXIO test log file
-#      if: failure()
-#      uses: actions/upload-artifact@82c141cc518b40d92cc801eee768e7aafc9c2fa2
-#      with:
-#        name: test-report-trexio-macos
-#        path: trexio/test-suite.log
-#
-#    - name: Build QMCkl
-#      run: |
-#        export PKG_CONFIG_PATH=${PWD}/trexio/_install/lib/pkgconfig:$PKG_CONFIG_PATH
-#        ./autogen.sh
-#        ./configure CC=gcc-10 FC=gfortran-10 --enable-silent-rules
-#        make -j 4
-#
-#    - name: Run test
-#      run: make -j 4 check
-#
-#    - name: Archive test log file
-#      if: failure()
-#      uses: actions/upload-artifact@82c141cc518b40d92cc801eee768e7aafc9c2fa2
-#      with:
-#        name: test-report-macos
-#        path: test-suite.log
+  macos:
+
+    runs-on: macos-12
+    name: x86 MacOS 12
+
+    steps:
+    - uses: actions/checkout@e2f20e631ae6d7dd3b768f56a5d2af784dd54791
+
+    - name: Install dependencies
+      run: |
+         brew install emacs
+         brew install automake
+         brew install hdf5
+         brew install gcc
+         brew install gfortran
+         brew --prefix hdf5
+
+    - name: Install the latest TREXIO from the GitHub clone
+      run: |
+        git clone https://github.com/TREX-CoE/trexio.git
+        cd trexio
+        ./autogen.sh
+        ./configure FC=gfortran-12 --enable-silent-rules
+        make -j 4
+        sudo make install
+
+    - name: Compile QMCkl in HPC mode
+      run: |
+        ./autogen.sh
+        mkdir _build_hpc
+        cd _build_hpc
+        ../configure --enable-hpc FC=gfortran-12 CC=gcc-12
+        make -j2
+
+    - name: Run test
+      run: make -j2 check
+      working-directory: _build_hpc
+
+    - name: Archive test log file
+      if: failure()
+      uses: actions/upload-artifact@82c141cc518b40d92cc801eee768e7aafc9c2fa2
+      with:
+        name: test-report-macos-x86
+        path: _build_hpc/test-suite.log
+
--- a/Makefile.am
+++ b/Makefile.am
@ -38,7 +38,7 @@ VERSION_MINOR   = @VERSION_MINOR@
 VERSION_PATCH   = @VERSION_PATCH@

 SUBDIRS =
-CLEANFILES = qmckl.mod qmckl_verificarlo_f.mod
+CLEANFILES = qmckl.mod qmckl_verificarlo_f.mod qmckl_constants.mod
 EXTRA_DIST =

 pkgconfigdir = $(libdir)/pkgconfig
@ -46,6 +46,7 @@ pkgconfig_DATA = pkgconfig/qmckl.pc

 qmckl_h = include/qmckl.h
 qmckl_f = include/qmckl_f.F90
+qmckl_fo = include/qmckl_f.o
 include_HEADERS = $(qmckl_h) $(qmckl_f)

 header_tests  = tests/chbrclf.h tests/n2.h
@ -61,7 +62,7 @@ lib_LTLIBRARIES = src/libqmckl.la
 src_libqmckl_la_SOURCES = $(qmckl_h) $(qmckl_f) $(C_FILES) $(F_FILES) $(H_PRIVATE_FUNC_FILES) $(H_PRIVATE_TYPE_FILES)
 src_libqmckl_la_LDFLAGS = $(LDFLAGS)

-CLEANFILES+=$(test_qmckl_fo) $(src_qmckl_fo) $(test_qmckl_o) $(src_qmckl_o) $(FH_TYPE_FILES) $(FH_FUNC_FILES)
+CLEANFILES+=$(qmckl_fo) $(test_qmckl_o) $(FH_TYPE_FILES) $(FH_FUNC_FILES)


 include generated.mk
--- a/README.md
+++ b/README.md
@ -2,8 +2,7 @@

 <img src="https://trex-coe.eu/sites/default/files/styles/responsive_no_crop/public/2022-01/QMCkl%20code.png?itok=UvOUClA5" width=200>

-![Build Status](https://github.com/TREX-CoE/qmckl/workflows/test-build/badge.svg?branch=master)
-
+![Build Status](https://github.com/TREX-CoE/qmckl/actions/workflows/test-build.yml/badge.svg?branch=master)
 The domain of quantum chemistry needs a library in which the main
 kernels of Quantum Monte Carlo (QMC) methods are implemented. In the
 library proposed in this project, we expose the main algorithms in a
--- a/autogen.sh
+++ b/autogen.sh
@ -1,5 +1,5 @@
-#!/bin/bash
+#!/bin/sh

 export srcdir="."
 python3 ${srcdir}/tools/build_makefile.py
-autoreconf -i -Wall --no-recursive
+autoreconf -vi -Wall --no-recursive
--- a/configure.ac
+++ b/configure.ac
@ -35,7 +35,7 @@

 AC_PREREQ([2.69])

-AC_INIT([qmckl],[0.3.1],[https://github.com/TREX-CoE/qmckl/issues],[],[https://trex-coe.github.io/qmckl/index.html])
+AC_INIT([qmckl],[1.0.0],[https://github.com/TREX-CoE/qmckl/issues],[],[https://trex-coe.github.io/qmckl/index.html])
 AC_CONFIG_AUX_DIR([tools])
 AM_INIT_AUTOMAKE([subdir-objects color-tests parallel-tests silent-rules 1.11])

@ -58,6 +58,22 @@ AS_IF([test "x$with_ifort" = "xyes"], [
  FCFLAGS="-march=native -ip -Ofast -ftz -finline -g -mkl=sequential" ])

 # Intel C compiler
+AC_ARG_WITH([icx],
+            [AS_HELP_STRING([--with-icx],
+                            [Use Intel C compiler])],
+            [with_icx=$withval],
+            [with_icx=no])
+
+AS_IF([test "x$with_icx" = "xyes"], [
+  CC=icx
+  CFLAGS="-march=native -Ofast -ftz -finline -g -qmkl=sequential" ])
+
+AS_IF([test "x$with_icx.$with_ifort" = "xyes.yes"], [
+  ax_blas_ok="yes"
+  ax_lapack_ok="yes"
+  BLAS_LIBS=""
+  LAPACK_LIBS=""])
+
 AC_ARG_WITH([icc],
            [AS_HELP_STRING([--with-icc],
                            [Use Intel C compiler])],
@ -215,31 +231,6 @@ PKG_CFLAGS2="$PKG_CFLAGS2 $TREXIO_CFLAGS"
 PKG_LIBS2="$PKG_LIBS2 $TREXIO_LIBS"


-# Checking SIMD
-AC_MSG_CHECKING([for SIMD length])
-AC_RUN_IFELSE(
-    [AC_LANG_PROGRAM([],[
-int simd=1;
-#if defined(__AVX512F__)
-   simd=8;
-#elif defined(__AVX__)
-   simd=4;
-#elif defined(__SSE2__)
-   simd=2;
-#elif defined(__ARM_NEON__)
-   simd=2;
-#endif
-  return simd;
-])],
-[ SIMD_LENGTH=1 ],
-[ SIMD_LENGTH=$? ],
-[ SIMD_LENGTH=1 ]
-)
-AC_MSG_RESULT([$SIMD_LENGTH])
-AC_DEFINE_UNQUOTED([SIMD_LENGTH], [$SIMD_LENGTH], [Length of SIMD vectors])
-
-
-
 # QMCKLDGEMM
 AC_ARG_WITH([qmckldgemm],
            [AS_HELP_STRING([--with-qmckldgemm],
@ -307,6 +298,18 @@ AC_ARG_ENABLE([hpc],
 AS_IF([test "x$enable_hpc" = "xyes"],
      [AC_DEFINE([HAVE_HPC], [1], [Activate HPC routines])])

+
+AC_ARG_ENABLE([fpe],
+              [AS_HELP_STRING([--enable-fpe],
+                              [Enable floating-point exceptions])],
+              [enable_fpe=$enableval],
+              [enable_fpe=no])
+
+AS_IF([test "x$enable_fpe" = "xyes"],
+      [AC_DEFINE([HAVE_FPE], [1], [Activate floating-point exceptions])])
+
+
+
 AC_ARG_ENABLE([doc],
            [AS_HELP_STRING([--disable-doc],
                            [Disable documentation])],
@ -336,62 +339,6 @@ AS_IF([test "$FC" = "verificarlo-f"], [
  FCFLAGS="-Mpreprocess $FCFLAGS"
 ])

-## Enable GPU offloading
-
-# GPU offloading
-AC_ARG_ENABLE(gpu, [AS_HELP_STRING([--enable-gpu],[openmp|openacc : Use GPU-offloaded functions])], enable_gpu=$enableval, enable_gpu=no)
-AS_IF([test "x$enable_gpu" = "xyes"], [enable_gpu="openmp"])
-
-# OpenMP offloading
-HAVE_OPENMP_OFFLOAD="no"
-AS_IF([test "x$enable_gpu" = "xopenmp"], [
-  AC_DEFINE([HAVE_OPENMP_OFFLOAD], [1], [If defined, activate OpenMP-offloaded routines])
-  HAVE_OPENMP_OFFLOAD="yes"
-  AS_CASE([$CC],
-    [*gcc*], [CFLAGS="$CFLAGS -fopenmp"],
-    [*nvc*], [CFLAGS="$CFLAGS -mp=gpu"],
-    [])
-  AS_CASE([$FC],
-    [*gfortran*],  [FCFLAGS="$FCFLAGS -fopenmp"],
-    [*nvfortran*], [FCFLAGS="$FCFLAGS -mp=gpu"],
-    [])
-  ]
-)
-
-# OpenMP offloading
-HAVE_OPENACC_OFFLOAD="no"
-AS_IF([test "x$enable_gpu" = "xopenacc"], [
-  AC_DEFINE([HAVE_OPENACC_OFFLOAD], [1], [If defined, activate OpenACC-offloaded routines])
-  HAVE_OPENACC_OFFLOAD="yes"
-  AS_CASE([$CC],
-    [*gcc*], [CFLAGS="$CFLAGS -fopenacc"],
-    [*nvc*], [CFLAGS="$CFLAGS -acc=gpu"],
-    [])
-  AS_CASE([$FC],
-    [*gfortran*],  [FCFLAGS="$FCFLAGS -fopenacc"],
-    [*nvfortran*], [FCFLAGS="$FCFLAGS -acc=gpu"],
-    [])
-  ]
-
-])
-
-# cuBLAS offloading
-AC_ARG_WITH(cublas, [AS_HELP_STRING([--with-cublas],[Use cuBLAS-offloaded functions])], HAVE_CUBLAS_OFFLOAD=$withval, HAVE_CUBLAS_OFFLOAD=no)
-AS_IF([test "x$HAVE_CUBLAS_OFFLOAD" = "xyes"], [
-  AC_DEFINE([HAVE_CUBLAS_OFFLOAD], [1], [If defined, activate cuBLAS-offloaded routines])
-  HAVE_OPENACC_OFFLOAD="yes"
-  AS_CASE([$CC],
-    [*gcc*], [CFLAGS="$CFLAGS -fopenmp"
-              LDFLAGS="-lcublas"],
-    [*nvc*], [CFLAGS="$CFLAGS -mp=gpu -cudalib=cublas"],
-    [])
-  AS_CASE([$FC],
-    [*gfortran*],  [FCFLAGS="$FCFLAGS -fopenmp"],
-    [*nvfortran*], [FCFLAGS="$FCFLAGS -mp=gpu -cudalib=cublas"],
-    [])
-  ]
-])
-
 AC_ARG_ENABLE(malloc-trace, [AS_HELP_STRING([--enable-malloc-trace],[use debug malloc/free])], ok=$enableval, ok=no)
 AS_IF([test "x$ok" = "xyes"], [
 	AC_DEFINE(MALLOC_TRACE,"malloc_trace.dat",[Define to use debugging malloc/free])
@ -419,12 +366,12 @@ AC_ARG_ENABLE(debug, [AS_HELP_STRING([--enable-debug],[compile for debugging])],
 AS_IF([test "x$ok" = "xyes"], [
        AS_IF([test "x$GCC" = "xyes"], [
                CPPFLAGS="-Wdate-time -D_FORTIFY_SOURCE=2"
-	        CFLAGS="$CFLAGS \
-g -Wall -W -Wbad-function-cast -Wcast-qual -Warray-bounds -Wdisabled-optimization \
-fsanitize=address -fno-omit-frame-pointer -fstack-protector-strong -Wformat -Werror=format-security \
+	        CFLAGS="$CFLAGS -g \
+-Wall -W -Wbad-function-cast -Wcast-qual -Warray-bounds -Wdisabled-optimization \
+-fno-omit-frame-pointer -fstack-protector-strong -Wformat -Werror=format-security \
 -Wpointer-arith -Wcast-align -Wpedantic -Wextra -Walloc-zero -Werror \
 "
-		LDFLAGS="$LDFLAGS -fsanitize=address"
+		LDFLAGS="$LDFLAGS"
        ])
        AS_IF([test "x$GFC" = "xyes"], [
 	        FCFLAGS="$FCFLAGS \
@ -498,6 +445,144 @@ AS_IF([test "x${QMCKL_DEVEL}" != "x"], [
  ] )
 ])

+# Checking SIMD
+AC_MSG_CHECKING([for SIMD length])
+SIMD_LENGTH=1
+AC_RUN_IFELSE(
+    [AC_LANG_PROGRAM([],[
+int simd=1;
+#if defined(__AVX512F__)
+   simd=8;
+#elif defined(__AVX2__)
+   simd=4;
+#elif defined(__AVX__)
+   simd=4;
+#elif defined(__SSE2__)
+   simd=2;
+#elif defined(__ARM_NEON__)
+   simd=2;
+#endif
+  return simd;
+])], [SIMD_LENGTH=1],
+[ AS_CASE([$?],
+  [1], [SIMD_LENGTH=1],
+  [2], [SIMD_LENGTH=2],
+  [4], [SIMD_LENGTH=4],
+  [8], [SIMD_LENGTH=8],
+  [16], [SIMD_LENGTH=16],
+  [SIMD_LENGTH=1])],
+  [SIMD_LENGTH=1]
+)
+AC_MSG_RESULT([$SIMD_LENGTH])
+AC_DEFINE_UNQUOTED([SIMD_LENGTH], [$SIMD_LENGTH], [Length of SIMD vectors])
+
+# Checking IVDEP
+ivdep=""
+AC_MSG_CHECKING([for ivdep pragma])
+AC_COMPILE_IFELSE([AC_LANG_PROGRAM([[
+    #include <stdio.h>
+]], [[
+    int main() {
+        #pragma ivdep
+        for (int i = 0; i < 10; ++i) {
+            printf("Testing: %d\n", i);
+        }
+        return 0;
+    }
+]])],
+    [ivdep='_Pragma("ivdep")'], [
+])
+
+AS_IF([test "x$ivdep" = "x"], [
+  AC_COMPILE_IFELSE([AC_LANG_PROGRAM([[
+      #include <stdio.h>
+  ]], [[
+      int main() {
+          #pragma clang loop vectorize(enable)
+          for (int i = 0; i < 10; ++i) {
+              printf("Testing: %d\n", i);
+          }
+          return 0;
+      }
+  ]])],
+      [ivdep='_Pragma("clang loop vectorize(enable)")'], [
+  ])
+])
+
+AS_IF([test "x$ivdep" = "x"], [
+  AC_COMPILE_IFELSE([AC_LANG_PROGRAM([[
+      #include <stdio.h>
+  ]], [[
+      int main() {
+          #pragma GCC ivdep
+          for (int i = 0; i < 10; ++i) {
+              printf("Testing: %d\n", i);
+          }
+          return 0;
+      }
+  ]])],
+      [ivdep='_Pragma("GCC ivdep")'], [
+  ])
+])
+
+AC_DEFINE_UNQUOTED([IVDEP], [$ivdep], [IVDEP pragma])
+AS_IF([test "x$ivdep" = "x"], [
+  ivdep="no"
+])
+AC_MSG_RESULT([$ivdep])
+
+
+# Checking ALIGNED
+
+AC_MSG_CHECKING([for posix_memalign])
+AC_COMPILE_IFELSE([AC_LANG_PROGRAM([[
+#include <stdlib.h>
+#include <stdio.h>
+]], [[
+    int main() {
+      void *ptr;
+      int ret = posix_memalign(&ptr, 64, 1024);
+      if (ret != 0) {
+        return EXIT_FAILURE;
+      }
+      free(ptr);
+      return 0;
+    }
+]])],
+    [have_posix_memalign=yes], [have_posix_memalign=no
+])
+AS_IF([test "x$have_posix_memalign" = "xyes"], [
+  AC_DEFINE([HAVE_POSIX_MEMALIGN], [1], [Define to 1 if you have the posix_memalign function.])
+])
+AC_MSG_RESULT([$have_posix_memalign])
+
+aligned=""
+AC_MSG_CHECKING([for vector aligned pragma])
+AC_COMPILE_IFELSE([AC_LANG_PROGRAM([[
+]], [[
+    int main() {
+        double __attribute__((aligned(8))) a[10] ;
+        #pragma vector aligned
+        for (int i = 0; i < 10; ++i) {
+            a[i] = (double) i;
+        }
+        return 0;
+    }
+]])],
+    [aligned='_Pragma("vector aligned")'], [
+])
+
+AS_IF([test "x$have_posix_memalign" = "xno"], [
+  aligned=""
+])
+
+AC_DEFINE_UNQUOTED([ALIGNED], [$aligned], [VECTOR ALIGNED pragma])
+AS_IF([test "x$aligned" = "x"], [
+  aligned="no"
+])
+AC_MSG_RESULT([$aligned])
+
+

 #PKG-CONFIG
 #mkl-dynamic-lp64-seq
--- a/m4/ax_blas.m4
+++ b/m4/ax_blas.m4
@ -65,6 +65,10 @@

 #serial 12

+# Updates for QMCKl:
+# - sgemm -> dgemm
+# - Include check for ARMPL
+
 AU_ALIAS([ACX_BLAS], [AX_BLAS])
 AC_DEFUN([AX_BLAS], [
 AC_PREREQ(2.50)
@ -91,8 +95,8 @@ LIBS="$LIBS $FLIBS"
 if test $ax_blas_ok = no; then
 if test "x$BLAS_LIBS" != x; then
 	save_LIBS="$LIBS"; LIBS="$BLAS_LIBS $LIBS"
-	AC_MSG_CHECKING([for $sgemm in $BLAS_LIBS])
-	AC_TRY_LINK_FUNC($sgemm, [ax_blas_ok=yes], [BLAS_LIBS=""])
+	AC_MSG_CHECKING([for $dgemm in $BLAS_LIBS])
+	AC_TRY_LINK_FUNC($dgemm, [ax_blas_ok=yes], [BLAS_LIBS=""])
 	AC_MSG_RESULT($ax_blas_ok)
 	LIBS="$save_LIBS"
 fi
@ -101,22 +105,22 @@ fi
 # BLAS linked to by default?  (happens on some supercomputers)
 if test $ax_blas_ok = no; then
 	save_LIBS="$LIBS"; LIBS="$LIBS"
-	AC_MSG_CHECKING([if $sgemm is being linked in already])
-	AC_TRY_LINK_FUNC($sgemm, [ax_blas_ok=yes])
+	AC_MSG_CHECKING([if $dgemm is being linked in already])
+	AC_TRY_LINK_FUNC($dgemm, [ax_blas_ok=yes])
 	AC_MSG_RESULT($ax_blas_ok)
 	LIBS="$save_LIBS"
 fi

 # BLAS in OpenBLAS library? (http://xianyi.github.com/OpenBLAS/)
 if test $ax_blas_ok = no; then
-	AC_CHECK_LIB(openblas, $sgemm, [ax_blas_ok=yes
+	AC_CHECK_LIB(openblas, $dgemm, [ax_blas_ok=yes
 			                BLAS_LIBS="-lopenblas"])
 fi

 # BLAS in ATLAS library? (http://math-atlas.sourceforge.net/)
 if test $ax_blas_ok = no; then
 	AC_CHECK_LIB(atlas, ATL_xerbla,
-		[AC_CHECK_LIB(f77blas, $sgemm,
+		[AC_CHECK_LIB(f77blas, $dgemm,
 		[AC_CHECK_LIB(cblas, cblas_dgemm,
 			[ax_blas_ok=yes
 			 BLAS_LIBS="-lcblas -lf77blas -latlas"],
@ -136,33 +140,41 @@ fi

 # BLAS in Intel MKL library?
 if test $ax_blas_ok = no; then
-	AC_CHECK_LIB(mkl, $sgemm, [ax_blas_ok=yes;BLAS_LIBS="-lmkl -lguide -lpthread"],,[-lguide -lpthread])
+	AC_CHECK_LIB(mkl, $dgemm, [ax_blas_ok=yes;BLAS_LIBS="-lmkl -lguide -lpthread"],,[-lguide -lpthread])
+fi
+
+if test $ax_blas_ok = no; then
+	AC_CHECK_LIB(mkl_gnu_thread, $dgemm, [ax_blas_ok=yes;BLAS_LIBS="-lmkl_gnu_thread -lmkl_core -ldl"],,[-lmkl_core -ldl])
+fi
+
+if test $ax_blas_ok = no; then
+	AC_CHECK_LIB(mkl_sequential, $dgemm, [ax_blas_ok=yes;BLAS_LIBS="-lmkl_sequential -lmkl_core -ldl"],,[-lmkl_core -ldl])
 fi

 # BLAS in Apple vecLib library?
 if test $ax_blas_ok = no; then
 	save_LIBS="$LIBS"; LIBS="-framework vecLib $LIBS"
-	AC_MSG_CHECKING([for $sgemm in -framework vecLib])
-	AC_TRY_LINK_FUNC($sgemm, [ax_blas_ok=yes;BLAS_LIBS="-framework vecLib"])
+	AC_MSG_CHECKING([for $dgemm in -framework vecLib])
+	AC_TRY_LINK_FUNC($dgemm, [ax_blas_ok=yes;BLAS_LIBS="-framework vecLib"])
 	AC_MSG_RESULT($ax_blas_ok)
 	LIBS="$save_LIBS"
 fi

 # BLAS in Alpha CXML library?
 if test $ax_blas_ok = no; then
-	AC_CHECK_LIB(cxml, $sgemm, [ax_blas_ok=yes;BLAS_LIBS="-lcxml"])
+	AC_CHECK_LIB(cxml, $dgemm, [ax_blas_ok=yes;BLAS_LIBS="-lcxml"])
 fi

 # BLAS in Alpha DXML library? (now called CXML, see above)
 if test $ax_blas_ok = no; then
-	AC_CHECK_LIB(dxml, $sgemm, [ax_blas_ok=yes;BLAS_LIBS="-ldxml"])
+	AC_CHECK_LIB(dxml, $dgemm, [ax_blas_ok=yes;BLAS_LIBS="-ldxml"])
 fi

 # BLAS in Sun Performance library?
 if test $ax_blas_ok = no; then
 	if test "x$GCC" != xyes; then # only works with Sun CC
 		AC_CHECK_LIB(sunmath, acosp,
-			[AC_CHECK_LIB(sunperf, $sgemm,
+			[AC_CHECK_LIB(sunperf, $dgemm,
 				[BLAS_LIBS="-xlic_lib=sunperf -lsunmath"
                                 ax_blas_ok=yes],[],[-lsunmath])])
 	fi
@ -170,26 +182,46 @@ fi

 # BLAS in SCSL library?  (SGI/Cray Scientific Library)
 if test $ax_blas_ok = no; then
-	AC_CHECK_LIB(scs, $sgemm, [ax_blas_ok=yes; BLAS_LIBS="-lscs"])
+	AC_CHECK_LIB(scs, $dgemm, [ax_blas_ok=yes; BLAS_LIBS="-lscs"])
 fi

 # BLAS in SGIMATH library?
 if test $ax_blas_ok = no; then
-	AC_CHECK_LIB(complib.sgimath, $sgemm,
+	AC_CHECK_LIB(complib.sgimath, $dgemm,
 		     [ax_blas_ok=yes; BLAS_LIBS="-lcomplib.sgimath"])
 fi

 # BLAS in IBM ESSL library? (requires generic BLAS lib, too)
 if test $ax_blas_ok = no; then
-	AC_CHECK_LIB(blas, $sgemm,
-		[AC_CHECK_LIB(essl, $sgemm,
+	AC_CHECK_LIB(blas, $dgemm,
+		[AC_CHECK_LIB(essl, $dgemm,
 			[ax_blas_ok=yes; BLAS_LIBS="-lessl -lblas"],
 			[], [-lblas $FLIBS])])
 fi

+# BLAS in ARMPL library?
+if test $ax_blas_ok = no; then
+	AC_CHECK_LIB(armpl_lp64_mp, $dgemm, [ax_blas_ok=yes;BLAS_LIBS="-larmpl_lp64_mp"])
+fi
+
+if test $ax_blas_ok = no; then
+	AC_CHECK_LIB(armpl_lp64, $dgemm, [ax_blas_ok=yes;BLAS_LIBS="-larmpl_lp64"])
+fi
+
+
+# BLAS in ACML?
+if test $ax_blas_ok = no; then
+	AC_CHECK_LIB(acml_mp, $dgemm, [ax_blas_ok=yes; BLAS_LIBS="-lacml_mp"])
+fi
+
+if test $ax_blas_ok = no; then
+	AC_CHECK_LIB(acml, $dgemm, [ax_blas_ok=yes; BLAS_LIBS="-lacml"])
+fi
+
+
 # Generic BLAS library?
 if test $ax_blas_ok = no; then
-	AC_CHECK_LIB(blas, $sgemm, [ax_blas_ok=yes; BLAS_LIBS="-lblas"])
+	AC_CHECK_LIB(blas, $dgemm, [ax_blas_ok=yes; BLAS_LIBS="-lblas"])
 fi

 AC_SUBST(BLAS_LIBS)
--- a/m4/ax_openmp.m4
+++ b/m4/ax_openmp.m4
@ -1,5 +1,5 @@
 # ===========================================================================
-#         http://www.gnu.org/software/autoconf-archive/ax_openmp.html
+#        https://www.gnu.org/software/autoconf-archive/ax_openmp.html
 # ===========================================================================
 #
 # SYNOPSIS
@ -38,6 +38,8 @@
 # LICENSE
 #
 #   Copyright (c) 2008 Steven G. Johnson <stevenj@alum.mit.edu>
+#   Copyright (c) 2015 John W. Peterson <jwpeterson@gmail.com>
+#   Copyright (c) 2016 Nick R. Papior <nickpapior@gmail.com>
 #
 #   This program is free software: you can redistribute it and/or modify it
 #   under the terms of the GNU General Public License as published by the
@ -50,7 +52,7 @@
 #   Public License for more details.
 #
 #   You should have received a copy of the GNU General Public License along
-#   with this program. If not, see <http://www.gnu.org/licenses/>.
+#   with this program. If not, see <https://www.gnu.org/licenses/>.
 #
 #   As a special exception, the respective Autoconf Macro's copyright owner
 #   gives unlimited permission to copy, distribute and modify the configure
@ -65,16 +67,19 @@
 #   modified version of the Autoconf Macro, you may extend this special
 #   exception to the GPL to apply to your modified version as well.

-#serial 8
+#serial 14

 AC_DEFUN([AX_OPENMP], [
-AC_PREREQ(2.59) dnl for _AC_LANG_PREFIX
+AC_PREREQ([2.69]) dnl for _AC_LANG_PREFIX

 AC_CACHE_CHECK([for OpenMP flag of _AC_LANG compiler], ax_cv_[]_AC_LANG_ABBREV[]_openmp, [save[]_AC_LANG_PREFIX[]FLAGS=$[]_AC_LANG_PREFIX[]FLAGS
 ax_cv_[]_AC_LANG_ABBREV[]_openmp=unknown
-# Flags to try:  -fopenmp (gcc), -openmp (icc), -mp (SGI & PGI),
-#                -xopenmp (Sun), -omp (Tru64), -qsmp=omp (AIX), none
-ax_openmp_flags="-fopenmp -openmp -mp -xopenmp -omp -qsmp=omp none"
+# Flags to try:  -fopenmp (gcc), -mp (SGI & PGI),
+#                -qopenmp (icc>=15), -openmp (icc),
+#                -xopenmp (Sun), -omp (Tru64),
+#                -qsmp=omp (AIX),
+#                none
+ax_openmp_flags="-fopenmp -openmp -qopenmp -mp -xopenmp -omp -qsmp=omp none"
 if test "x$OPENMP_[]_AC_LANG_PREFIX[]FLAGS" != x; then
  ax_openmp_flags="$OPENMP_[]_AC_LANG_PREFIX[]FLAGS $ax_openmp_flags"
 fi
@ -83,8 +88,27 @@ for ax_openmp_flag in $ax_openmp_flags; do
    none) []_AC_LANG_PREFIX[]FLAGS=$save[]_AC_LANG_PREFIX[] ;;
    *) []_AC_LANG_PREFIX[]FLAGS="$save[]_AC_LANG_PREFIX[]FLAGS $ax_openmp_flag" ;;
  esac
-  AC_TRY_LINK_FUNC(omp_set_num_threads,
-	[ax_cv_[]_AC_LANG_ABBREV[]_openmp=$ax_openmp_flag; break])
+  AC_LINK_IFELSE([AC_LANG_SOURCE([[
+@%:@include <omp.h>
+
+static void
+parallel_fill(int * data, int n)
+{
+  int i;
+@%:@pragma omp parallel for
+  for (i = 0; i < n; ++i)
+    data[i] = i;
+}
+
+int
+main(void)
+{
+  int arr[100000];
+  omp_set_num_threads(2);
+  parallel_fill(arr, 100000);
+  return 0;
+}
+]])],[ax_cv_[]_AC_LANG_ABBREV[]_openmp=$ax_openmp_flag; break],[])
 done
 []_AC_LANG_PREFIX[]FLAGS=$save[]_AC_LANG_PREFIX[]FLAGS
 ])
@ -97,3 +121,4 @@ else
  m4_default([$1], [AC_DEFINE(HAVE_OPENMP,1,[Define if OpenMP is enabled])])
 fi
 ])dnl AX_OPENMP
+
--- a/org/README.org
+++ b/org/README.org
@ -17,24 +17,25 @@ grep TITLE $(cat table_of_contents) | tr ':' ' '

  #+RESULTS: toc
  | qmckl.org                           | #+TITLE | Introduction              |             |         |
-  | qmckl_ao.org                        | #+TITLE | Atomic                    | Orbitals    |         |
-  | qmckl_blas.org                      | #+TITLE | BLAS                      | functions   |         |
  | qmckl_context.org                   | #+TITLE | Context                   |             |         |
-  | qmckl_determinant.org               | #+TITLE | Slater                    | Determinant |         |
-  | qmckl_distance.org                  | #+TITLE | Inter-particle            | distances   |         |
-  | qmckl_electron.org                  | #+TITLE | Electrons                 |             |         |
  | qmckl_error.org                     | #+TITLE | Error                     | handling    |         |
-  | qmckl_jastrow.org                   | #+TITLE | Jastrow                   | Factor      |         |
-  | qmckl_local_energy.org              | #+TITLE | Local                     | Energy      |         |
+  | qmckl_blas.org                      | #+TITLE | BLAS                      | functions   |         |
  | qmckl_memory.org                    | #+TITLE | Memory                    | management  |         |
-  | qmckl_mo.org                        | #+TITLE | Molecular                 | Orbitals    |         |
  | qmckl_numprec.org                   | #+TITLE | Numerical                 | precision   |         |
+  | qmckl_point.org                     | #+TITLE | Point                     |             |         |
  | qmckl_nucleus.org                   | #+TITLE | Nucleus                   |             |         |
+  | qmckl_electron.org                  | #+TITLE | Electrons                 |             |         |
+  | qmckl_distance.org                  | #+TITLE | Inter-particle            | distances   |         |
+  | qmckl_ao.org                        | #+TITLE | Atomic                    | Orbitals    |         |
+  | qmckl_mo.org                        | #+TITLE | Molecular                 | Orbitals    |         |
+  | qmckl_determinant.org               | #+TITLE | Slater                    | Determinant |         |
  | qmckl_sherman_morrison_woodbury.org | #+TITLE | Sherman-Morrison-Woodbury |             |         |
-  | qmckl_utils.org                     | #+TITLE | Utility                   | functions   |         |
+  | qmckl_jastrow_champ.org             | #+TITLE | CHAMP                     | Jastrow     | Factor  |
+  | qmckl_local_energy.org              | #+TITLE | Local                     | Energy      |         |
  | qmckl_trexio.org                    | #+TITLE | TREXIO                    | I/O         | library |
-  | qmckl_verificarlo.org               | #+TITLE | Verificarlo               | CI          |         |
  | qmckl_tests.org                     | #+TITLE | Data                      | for         | Tests   |
+  | qmckl_verificarlo.org               | #+TITLE | Verificarlo               | CI          |         |
+  | qmckl_examples.org                  | #+TITLE | Code                      | examples    |         |

  #+begin_src python :var data=toc :exports results :results raw
 result = []
@ -47,24 +48,25 @@ return '\n'.join(result)

  #+RESULTS:
   - [[./qmckl.html][Introduction]]
-   - [[./qmckl_ao.html][Atomic Orbitals]]
-   - [[./qmckl_blas.html][BLAS functions]]
   - [[./qmckl_context.html][Context]]
-   - [[./qmckl_determinant.html][Slater Determinant]]
-   - [[./qmckl_distance.html][Inter-particle distances]]
-   - [[./qmckl_electron.html][Electrons]]
   - [[./qmckl_error.html][Error handling]]
-   - [[./qmckl_jastrow.html][Jastrow Factor]]
-   - [[./qmckl_local_energy.html][Local Energy]]
+   - [[./qmckl_blas.html][BLAS functions]]
   - [[./qmckl_memory.html][Memory management]]
-   - [[./qmckl_mo.html][Molecular Orbitals]]
   - [[./qmckl_numprec.html][Numerical precision]]
+   - [[./qmckl_point.html][Point]]
   - [[./qmckl_nucleus.html][Nucleus]]
+   - [[./qmckl_electron.html][Electrons]]
+   - [[./qmckl_distance.html][Inter-particle distances]]
+   - [[./qmckl_ao.html][Atomic Orbitals]]
+   - [[./qmckl_mo.html][Molecular Orbitals]]
+   - [[./qmckl_determinant.html][Slater Determinant]]
   - [[./qmckl_sherman_morrison_woodbury.html][Sherman-Morrison-Woodbury]]
-   - [[./qmckl_utils.html][Utility functions]]
+   - [[./qmckl_jastrow_champ.html][CHAMP Jastrow Factor]]
+   - [[./qmckl_local_energy.html][Local Energy]]
   - [[./qmckl_trexio.html][TREXIO I/O library]]
-   - [[./qmckl_verificarlo.html][Verificarlo CI]]
   - [[./qmckl_tests.html][Data for Tests]]
+   - [[./qmckl_verificarlo.html][Verificarlo CI]]
+   - [[./qmckl_examples.html][Code examples]]


 --------------------------------
--- a/org/qmckl_ao.org
+++ b/org/qmckl_ao.org
--- a/org/qmckl_blas.org
+++ b/org/qmckl_blas.org
--- a/org/qmckl_context.org
+++ b/org/qmckl_context.org
@ -33,7 +33,7 @@ int main() {
 #include "qmckl_electron_private_type.h"
 #include "qmckl_ao_private_type.h"
 #include "qmckl_mo_private_type.h"
-#include "qmckl_jastrow_private_type.h"
+#include "qmckl_jastrow_champ_private_type.h"
 #include "qmckl_determinant_private_type.h"
 #include "qmckl_local_energy_private_type.h"
 #include "qmckl_point_private_func.h"
@ -41,6 +41,7 @@ int main() {
 #include "qmckl_electron_private_func.h"
 #include "qmckl_ao_private_func.h"
 #include "qmckl_mo_private_func.h"
+#include "qmckl_jastrow_champ_private_func.h"
 #include "qmckl_determinant_private_func.h"
 #include "qmckl_local_energy_private_func.h"
  #+end_src
@ -130,13 +131,13 @@ typedef struct qmckl_context_struct {
  qmckl_point_struct         point;

 /* -- Molecular system -- */
-  qmckl_nucleus_struct       nucleus;
-  qmckl_electron_struct      electron;
-  qmckl_ao_basis_struct      ao_basis;
-  qmckl_mo_basis_struct      mo_basis;
-  qmckl_jastrow_struct       jastrow;
-  qmckl_determinant_struct   det;
-  qmckl_local_energy_struct  local_energy;
+  qmckl_nucleus_struct        nucleus;
+  qmckl_electron_struct       electron;
+  qmckl_ao_basis_struct       ao_basis;
+  qmckl_mo_basis_struct       mo_basis;
+  qmckl_jastrow_champ_struct  jastrow_champ;
+  qmckl_determinant_struct    det;
+  qmckl_local_energy_struct   local_energy;

  /* To be implemented:
  ,*/
@ -303,7 +304,7 @@ qmckl_context qmckl_context_create() {
    rc = qmckl_init_determinant(context);
    assert (rc == QMCKL_SUCCESS);

-    rc = qmckl_init_jastrow(context);
+    rc = qmckl_init_jastrow_champ(context);
    assert (rc == QMCKL_SUCCESS);
  }

--- a/org/qmckl_determinant.org
+++ b/org/qmckl_determinant.org
@ -1225,7 +1225,7 @@ qmckl_check(context, rc);

 const double  * mo_coefficient          =  &(chbrclf_mo_coef[0]);

-rc = qmckl_set_mo_basis_coefficient(context, mo_coefficient);
+rc = qmckl_set_mo_basis_coefficient(context, mo_coefficient, chbrclf_mo_num*chbrclf_ao_num);
 qmckl_check(context, rc);

 assert(qmckl_mo_basis_provided(context));
--- a/org/qmckl_distance.org
+++ b/org/qmckl_distance.org
@ -90,20 +90,26 @@ int main() {
    #+end_src

    #+begin_src f90 :tangle (eval f)
-integer function qmckl_distance_sq_f(context, transa, transb, m, n, &
+function qmckl_distance_sq(context, transa, transb, m, n, &
     A, LDA, B, LDB, C, LDC) &
-     result(info)
-  use qmckl
+     bind(C) result(info)
+
+  use qmckl_constants
  implicit none
-  integer(qmckl_context)  , intent(in)  :: context
-  character  , intent(in)  :: transa, transb
-  integer*8  , intent(in)  :: m, n
-  integer*8  , intent(in)  :: lda
-  real*8     , intent(in)  :: A(lda,*)
-  integer*8  , intent(in)  :: ldb
-  real*8     , intent(in)  :: B(ldb,*)
-  integer*8  , intent(in)  :: ldc
-  real*8     , intent(out) :: C(ldc,*)
+
+  integer (qmckl_context) , intent(in)  , value :: context
+  character(c_char)   , intent(in)  , value :: transa
+  character(c_char)   , intent(in)  , value :: transb
+  integer (c_int64_t) , intent(in)  , value :: m
+  integer (c_int64_t) , intent(in)  , value :: n
+  integer (c_int64_t) , intent(in)  , value :: lda
+  integer (c_int64_t) , intent(in)  , value :: ldb
+  integer (c_int64_t) , intent(in)  , value :: ldc
+  real    (c_double ) , intent(in)          :: A(lda,*)
+  real    (c_double ) , intent(in)          :: B(ldb,*)
+  real    (c_double ) , intent(out)         :: C(ldc,n)
+
+  integer(qmckl_exit_code) :: info

  integer*8 :: i,j
  real*8    :: x, y, z
@ -157,12 +163,12 @@ integer function qmckl_distance_sq_f(context, transa, transb, m, n, &
     return
  endif

-  if (iand(transab,2) == 0 .and. LDA < 3) then
+  if (iand(transab,2) == 0 .and. LDB < 3) then
     info = QMCKL_INVALID_ARG_7
     return
  endif

-  if (iand(transab,2) == 2 .and. LDA < m) then
+  if (iand(transab,2) == 2 .and. LDB < n) then
     info = QMCKL_INVALID_ARG_7
     return
  endif
@ -216,7 +222,7 @@ integer function qmckl_distance_sq_f(context, transa, transb, m, n, &

  end select

-end function qmckl_distance_sq_f
+end function qmckl_distance_sq
    #+end_src

 *** Performance
@ -224,59 +230,29 @@ end function qmckl_distance_sq_f
    This function is more efficient when ~A~ and ~B~ are
    transposed.

-   #+CALL: generate_c_interface(table=qmckl_distance_sq_args,rettyp=get_value("FRetType"),fname=get_value("Name"))
-
-   #+RESULTS:
-   #+begin_src f90 :tangle (eval f) :comments org :exports none
-   integer(c_int32_t) function qmckl_distance_sq &
-       (context, transa, transb, m, n, A, lda, B, ldb, C, ldc) &
-       bind(C) result(info)
-
-     use, intrinsic :: iso_c_binding
-     implicit none
-
-     integer (c_int64_t) , intent(in)  , value :: context
-     character           , intent(in)  , value :: transa
-     character           , intent(in)  , value :: transb
-     integer (c_int64_t) , intent(in)  , value :: m
-     integer (c_int64_t) , intent(in)  , value :: n
-     real    (c_double ) , intent(in)          :: A(lda,*)
-     integer (c_int64_t) , intent(in)  , value :: lda
-     real    (c_double ) , intent(in)          :: B(ldb,*)
-     integer (c_int64_t) , intent(in)  , value :: ldb
-     real    (c_double ) , intent(out)         :: C(ldc,n)
-     integer (c_int64_t) , intent(in)  , value :: ldc
-
-     integer(c_int32_t), external :: qmckl_distance_sq_f
-     info = qmckl_distance_sq_f &
-            (context, transa, transb, m, n, A, lda, B, ldb, C, ldc)
-
-   end function qmckl_distance_sq
-   #+end_src
-
-   #+CALL: generate_f_interface(table=qmckl_distance_sq_args,rettyp=get_value("FRetType"),fname=get_value("Name"))
+   #+CALL: generate_f_interface(table=qmckl_distance_sq_args,fname=get_value("Name"))

   #+RESULTS:
   #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
   interface
-     integer(c_int32_t) function qmckl_distance_sq &
+     integer(qmckl_exit_code) function qmckl_distance_sq &
         (context, transa, transb, m, n, A, lda, B, ldb, C, ldc) &
         bind(C)
       use, intrinsic :: iso_c_binding
       import
       implicit none

-       integer (c_int64_t) , intent(in)  , value :: context
-       character           , intent(in)  , value :: transa
-       character           , intent(in)  , value :: transb
+       integer (qmckl_context), intent(in)  , value :: context
+       character(c_char  ) , intent(in)  , value :: transa
+       character(c_char  ) , intent(in)  , value :: transb
       integer (c_int64_t) , intent(in)  , value :: m
       integer (c_int64_t) , intent(in)  , value :: n
-       real    (c_double ) , intent(in)          :: A(lda,*)
       integer (c_int64_t) , intent(in)  , value :: lda
-       real    (c_double ) , intent(in)          :: B(ldb,*)
       integer (c_int64_t) , intent(in)  , value :: ldb
-       real    (c_double ) , intent(out)         :: C(ldc,n)
       integer (c_int64_t) , intent(in)  , value :: ldc
+       real    (c_double ) , intent(in)          :: A(lda,*)
+       real    (c_double ) , intent(in)          :: B(ldb,*)
+       real    (c_double ) , intent(out)         :: C(ldc,n)

     end function qmckl_distance_sq
   end interface
@ -288,7 +264,6 @@ integer(qmckl_exit_code) function test_qmckl_distance_sq(context) bind(C)

  use qmckl
  use qmckl_verificarlo_f
-  use iso_c_binding

  implicit none

@ -490,20 +465,23 @@ end function test_qmckl_distance_sq

 *** Source
    #+begin_src f90 :tangle (eval f)
-integer function qmckl_distance_f(context, transa, transb, m, n, &
+function qmckl_distance(context, transa, transb, m, n, &
     A, LDA, B, LDB, C, LDC) &
-     result(info)
-  use qmckl
+     bind(C) result(info)
+  use qmckl_constants
  implicit none
-  integer(qmckl_context)  , intent(in)  :: context
-  character  , intent(in)  :: transa, transb
-  integer*8  , intent(in)  :: m, n
-  integer*8  , intent(in)  :: lda
-  real*8     , intent(in)  :: A(lda,*)
-  integer*8  , intent(in)  :: ldb
-  real*8     , intent(in)  :: B(ldb,*)
-  integer*8  , intent(in)  :: ldc
-  real*8     , intent(out) :: C(ldc,*)
+  integer(qmckl_context), intent(in), value  :: context
+  character(c_char)   , intent(in)  , value :: transa
+  character(c_char)   , intent(in)  , value :: transb
+  integer (c_int64_t) , intent(in)  , value :: m
+  integer (c_int64_t) , intent(in)  , value :: n
+  integer (c_int64_t) , intent(in)  , value :: lda
+  integer (c_int64_t) , intent(in)  , value :: ldb
+  integer (c_int64_t) , intent(in)  , value :: ldc
+  real    (c_double ) , intent(in)          :: A(lda,*)
+  real    (c_double ) , intent(in)          :: B(ldb,*)
+  real    (c_double ) , intent(out)         :: C(ldc,n)
+  integer (qmckl_exit_code) :: info

  integer*8 :: i,j
  real*8    :: x, y, z
@ -558,16 +536,6 @@ integer function qmckl_distance_f(context, transa, transb, m, n, &
     return
  endif

-  if (iand(transab,2) == 0 .and. LDA < 3) then
-     info = QMCKL_INVALID_ARG_7
-     return
-  endif
-
-  if (iand(transab,2) == 2 .and. LDA < m) then
-     info = QMCKL_INVALID_ARG_7
-     return
-  endif
-
  ! check for LDB
  if (iand(transab,1) == 0 .and. LDB < 3) then
     info = QMCKL_INVALID_ARG_9
@ -579,16 +547,6 @@ integer function qmckl_distance_f(context, transa, transb, m, n, &
     return
  endif

-  if (iand(transab,2) == 0 .and. LDB < 3) then
-     info = QMCKL_INVALID_ARG_9
-     return
-  endif
-
-  if (iand(transab,2) == 2 .and. LDB < n) then
-     info = QMCKL_INVALID_ARG_9
-     return
-  endif
-
  ! check for LDC
  if (LDC < m) then
     info = QMCKL_INVALID_ARG_11
@ -648,73 +606,41 @@ integer function qmckl_distance_f(context, transa, transb, m, n, &

  end select

-end function qmckl_distance_f
+end function qmckl_distance
    #+end_src

-*** Performance
-
-    This function is more efficient when ~A~ and ~B~ are transposed.
-
-** C interface                                                     :noexport:
-
-   #+CALL: generate_c_interface(table=qmckl_distance_args,rettyp=get_value("FRetType"),fname=get_value("Name"))
-
-   #+RESULTS:
-   #+begin_src f90 :tangle (eval f) :comments org :exports none
-   integer(c_int32_t) function qmckl_distance &
-       (context, transa, transb, m, n, A, lda, B, ldb, C, ldc) &
-       bind(C) result(info)
-
-     use, intrinsic :: iso_c_binding
-     implicit none
-
-     integer (c_int64_t) , intent(in)  , value :: context
-     character           , intent(in)  , value :: transa
-     character           , intent(in)  , value :: transb
-     integer (c_int64_t) , intent(in)  , value :: m
-     integer (c_int64_t) , intent(in)  , value :: n
-     real    (c_double ) , intent(in)          :: A(lda,*)
-     integer (c_int64_t) , intent(in)  , value :: lda
-     real    (c_double ) , intent(in)          :: B(ldb,*)
-     integer (c_int64_t) , intent(in)  , value :: ldb
-     real    (c_double ) , intent(out)         :: C(ldc,n)
-     integer (c_int64_t) , intent(in)  , value :: ldc
-
-     integer(c_int32_t), external :: qmckl_distance_f
-     info = qmckl_distance_f &
-            (context, transa, transb, m, n, A, lda, B, ldb, C, ldc)
-
-   end function qmckl_distance
-   #+end_src
-
-   #+CALL: generate_f_interface(table=qmckl_distance_args,rettyp=get_value("FRetType"),fname=get_value("Name"))
+   #+CALL: generate_f_interface(table=qmckl_distance_args,fname="qmckl_distance")

   #+RESULTS:
   #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
   interface
-     integer(c_int32_t) function qmckl_distance &
+     integer(qmckl_exit_code) function qmckl_distance &
         (context, transa, transb, m, n, A, lda, B, ldb, C, ldc) &
         bind(C)
       use, intrinsic :: iso_c_binding
       import
       implicit none

-       integer (c_int64_t) , intent(in)  , value :: context
-       character           , intent(in)  , value :: transa
-       character           , intent(in)  , value :: transb
+       integer (qmckl_context), intent(in)  , value :: context
+       character(c_char  ) , intent(in)  , value :: transa
+       character(c_char  ) , intent(in)  , value :: transb
       integer (c_int64_t) , intent(in)  , value :: m
       integer (c_int64_t) , intent(in)  , value :: n
-       real    (c_double ) , intent(in)          :: A(lda,*)
       integer (c_int64_t) , intent(in)  , value :: lda
-       real    (c_double ) , intent(in)          :: B(ldb,*)
       integer (c_int64_t) , intent(in)  , value :: ldb
-       real    (c_double ) , intent(out)         :: C(ldc,n)
       integer (c_int64_t) , intent(in)  , value :: ldc
+       real    (c_double ) , intent(in)          :: A(lda,*)
+       real    (c_double ) , intent(in)          :: B(ldb,*)
+       real    (c_double ) , intent(out)         :: C(ldc,n)

     end function qmckl_distance
   end interface
   #+end_src

+*** Performance
+
+    This function is more efficient when ~A~ and ~B~ are transposed.
+
 *** Test                                                           :noexport:
    #+begin_src f90 :tangle (eval f_test)

@ -722,7 +648,6 @@ integer(qmckl_exit_code) function test_qmckl_dist(context) bind(C)

  use qmckl
  use qmckl_verificarlo_f
-  use iso_c_binding

  implicit none

@ -875,7 +800,7 @@ end function test_qmckl_dist
   pairs of points in two sets, one point within each set:

   \[
-   C_{ij} = \left( 1 - \exp \left(-\kappa C_{ij} \right) \right)/\kappa
+   C_{ij} = \frac{ 1 - e^{-\kappa r_{ij}}}{\kappa}
   \]

   If the input array is normal (~'N'~), the xyz coordinates are in
@ -934,21 +859,24 @@ end function test_qmckl_dist

 *** Source
    #+begin_src f90 :tangle (eval f)
-integer function qmckl_distance_rescaled_f(context, transa, transb, m, n, &
+function qmckl_distance_rescaled(context, transa, transb, m, n, &
     A, LDA, B, LDB, C, LDC, rescale_factor_kappa) &
-     result(info)
-  use qmckl
+     bind(C) result(info)
+  use qmckl_constants
  implicit none
-  integer(qmckl_context)  , intent(in)  :: context
-  character  , intent(in)  :: transa, transb
-  integer*8  , intent(in)  :: m, n
-  integer*8  , intent(in)  :: lda
-  real*8     , intent(in)  :: A(lda,*)
-  integer*8  , intent(in)  :: ldb
-  real*8     , intent(in)  :: B(ldb,*)
-  integer*8  , intent(in)  :: ldc
-  real*8     , intent(out) :: C(ldc,*)
-  real*8     , intent(in)  :: rescale_factor_kappa
+  integer (qmckl_context), intent(in)  , value :: context
+  character(c_char  ) , intent(in)  , value :: transa
+  character(c_char  ) , intent(in)  , value :: transb
+  integer (c_int64_t) , intent(in)  , value :: m
+  integer (c_int64_t) , intent(in)  , value :: n
+  integer (c_int64_t) , intent(in)  , value :: lda
+  integer (c_int64_t) , intent(in)  , value :: ldb
+  integer (c_int64_t) , intent(in)  , value :: ldc
+  real    (c_double ) , intent(in)  , value :: rescale_factor_kappa
+  real    (c_double ) , intent(in)          :: A(lda,*)
+  real    (c_double ) , intent(in)          :: B(ldb,*)
+  real    (c_double ) , intent(out)         :: C(ldc,n)
+  integer(qmckl_exit_code) :: info

  integer*8 :: i,j
  real*8    :: x, y, z, dist, rescale_factor_kappa_inv
@ -1005,27 +933,7 @@ integer function qmckl_distance_rescaled_f(context, transa, transb, m, n, &
     return
  endif

-  if (iand(transab,2) == 0 .and. LDA < 3) then
-     info = QMCKL_INVALID_ARG_7
-     return
-  endif
-
-  if (iand(transab,2) == 2 .and. LDA < m) then
-     info = QMCKL_INVALID_ARG_7
-     return
-  endif
-
  ! check for LDB
-  if (iand(transab,1) == 0 .and. LDB < 3) then
-     info = QMCKL_INVALID_ARG_9
-     return
-  endif
-
-  if (iand(transab,1) == 1 .and. LDB < n) then
-     info = QMCKL_INVALID_ARG_9
-     return
-  endif
-
  if (iand(transab,2) == 0 .and. LDB < 3) then
     info = QMCKL_INVALID_ARG_9
     return
@ -1095,7 +1003,7 @@ integer function qmckl_distance_rescaled_f(context, transa, transb, m, n, &

  end select

-end function qmckl_distance_rescaled_f
+end function qmckl_distance_rescaled
    #+end_src

 *** Performance
@ -1104,120 +1012,267 @@ end function qmckl_distance_rescaled_f

 ** C interface                                                     :noexport:

-   #+CALL: generate_c_interface(table=qmckl_distance_rescaled_args,rettyp=get_value("FRetType"),fname=get_value("Name"))
-
-   #+RESULTS:
-   #+begin_src f90 :tangle (eval f) :comments org :exports none
-   integer(c_int32_t) function qmckl_distance_rescaled &
-       (context, transa, transb, m, n, A, lda, B, ldb, C, ldc, rescale_factor_kappa) &
-       bind(C) result(info)
-
-     use, intrinsic :: iso_c_binding
-     implicit none
-
-     integer (c_int64_t) , intent(in)  , value :: context
-     character           , intent(in)  , value :: transa
-     character           , intent(in)  , value :: transb
-     integer (c_int64_t) , intent(in)  , value :: m
-     integer (c_int64_t) , intent(in)  , value :: n
-     real    (c_double ) , intent(in)          :: A(lda,*)
-     integer (c_int64_t) , intent(in)  , value :: lda
-     real    (c_double ) , intent(in)          :: B(ldb,*)
-     integer (c_int64_t) , intent(in)  , value :: ldb
-     real    (c_double ) , intent(out)         :: C(ldc,n)
-     integer (c_int64_t) , intent(in)  , value :: ldc
-     real    (c_double ) , intent(in)  , value :: rescale_factor_kappa
-
-     integer(c_int32_t), external :: qmckl_distance_rescaled_f
-     info = qmckl_distance_rescaled_f &
-            (context, transa, transb, m, n, A, lda, B, ldb, C, ldc, rescale_factor_kappa)
-
-   end function qmckl_distance_rescaled
-   #+end_src
-
-   #+CALL: generate_f_interface(table=qmckl_distance_rescaled_args,rettyp=get_value("FRetType"),fname=get_value("Name"))
+   #+CALL: generate_f_interface(table=qmckl_distance_rescaled_args,fname="qmckl_distance_rescaled")

   #+RESULTS:
   #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
   interface
-     integer(c_int32_t) function qmckl_distance_rescaled &
+     integer(qmckl_exit_code) function qmckl_distance_rescaled &
         (context, transa, transb, m, n, A, lda, B, ldb, C, ldc, rescale_factor_kappa) &
         bind(C)
       use, intrinsic :: iso_c_binding
       import
       implicit none

-       integer (c_int64_t) , intent(in)  , value :: context
-       character           , intent(in)  , value :: transa
-       character           , intent(in)  , value :: transb
+       integer (qmckl_context), intent(in)  , value :: context
+       character(c_char  ) , intent(in)  , value :: transa
+       character(c_char  ) , intent(in)  , value :: transb
       integer (c_int64_t) , intent(in)  , value :: m
       integer (c_int64_t) , intent(in)  , value :: n
-       real    (c_double ) , intent(in)          :: A(lda,*)
       integer (c_int64_t) , intent(in)  , value :: lda
-       real    (c_double ) , intent(in)          :: B(ldb,*)
       integer (c_int64_t) , intent(in)  , value :: ldb
-       real    (c_double ) , intent(out)         :: C(ldc,n)
       integer (c_int64_t) , intent(in)  , value :: ldc
       real    (c_double ) , intent(in)  , value :: rescale_factor_kappa
+       real    (c_double ) , intent(in)          :: A(lda,*)
+       real    (c_double ) , intent(in)          :: B(ldb,*)
+       real    (c_double ) , intent(out)         :: C(ldc,n)

     end function qmckl_distance_rescaled
   end interface
   #+end_src

 *** Test                                                           :noexport:
+
+   #+BEGIN_SRC python :results output :exports none
+import numpy as np
+
+kappa     = 0.6
+kappa_inv = 1./kappa
+
+# For H2O we have the following data:
+elec_num     = 10
+nucl_num     = 2
+up_num       = 5
+down_num     = 5
+nucl_coord = np.array([ [0.000000,  0.000000 ],
+   [0.000000,  0.000000 ],
+   [0.000000,  2.059801 ] ])
+
+elec_coord = np.array(   [[[-0.250655104764153      ,  0.503070975550133      ,  -0.166554344502303],
+     [-0.587812193472177      , -0.128751981129274      ,   0.187773606533075],
+     [ 1.61335569047166       , -0.615556732874863      ,  -1.43165470979934 ],
+     [-4.901239896295210E-003 , -1.120440036458986E-002 ,   1.99761909330422 ],
+     [ 0.766647499681200      , -0.293515395797937      ,   3.66454589201239 ],
+     [-0.127732483187947      , -0.138975497694196      ,  -8.669850480215846E-002],
+     [-0.232271834949124      , -1.059321673434182E-002 ,  -0.504862241464867],
+     [ 1.09360863531826       , -2.036103063808752E-003 ,  -2.702796910818986E-002],
+     [-0.108090166832043      ,  0.189161729653261      ,   2.15398313919894],
+     [ 0.397978144318712      , -0.254277292595981      ,   2.54553335476344]]])
+
+ee_distance_rescaled = \
+  np.array([ [(1.-np.exp(-kappa*np.linalg.norm(elec_coord[0,j,:]-elec_coord[0,i,:])))/kappa \
+             for i in range(elec_num) ]
+             for j in range(elec_num) ])
+
+en_distance_rescaled = \
+  np.array([ [(1.-np.exp(-kappa*np.linalg.norm(elec_coord[0,j,:]-nucl_coord[:,i])))/kappa \
+             for j in range(elec_num) ]
+             for i in range(nucl_num) ])
+
+print(ee_distance_rescaled)
+print()
+print(en_distance_rescaled)
+   #+END_SRC
+
+   #+RESULTS:
+   #+begin_example
+   [[0.         0.63475074 1.29816415 1.23147027 1.51933127 0.54402406
+     0.51452479 0.96538731 1.25878564 1.3722995 ]
+    [0.63475074 0.         1.35148664 1.13524156 1.48940503 0.4582292
+     0.62758076 1.06560856 1.179133   1.30763703]
+    [1.29816415 1.35148664 0.         1.50021375 1.59200788 1.23488312
+     1.20844259 1.0355537  1.52064535 1.53049239]
+    [1.23147027 1.13524156 1.50021375 0.         1.12016142 1.19158954
+     1.29762585 1.24824277 0.25292267 0.58609336]
+    [1.51933127 1.48940503 1.59200788 1.12016142 0.         1.50217017
+     1.54012828 1.48753895 1.10441805 0.84504381]
+    [0.54402406 0.4582292  1.23488312 1.19158954 1.50217017 0.
+     0.39417354 0.87009603 1.23838502 1.33419121]
+    [0.51452479 0.62758076 1.20844259 1.29762585 1.54012828 0.39417354
+     0.         0.95118809 1.33068934 1.41097406]
+    [0.96538731 1.06560856 1.0355537  1.24824277 1.48753895 0.87009603
+     0.95118809 0.         1.29422213 1.33222493]
+    [1.25878564 1.179133   1.52064535 0.25292267 1.10441805 1.23838502
+     1.33068934 1.29422213 0.         0.62196802]
+    [1.3722995  1.30763703 1.53049239 0.58609336 0.84504381 1.33419121
+     1.41097406 1.33222493 0.62196802 0.        ]]
+
+   [[0.49421587 0.52486545 1.23280503 1.16396998 1.49156627 0.1952946
+     0.4726453  0.80211227 1.21198526 1.31411513]
+    [1.24641375 1.15444238 1.50565145 0.06218339 1.10153184 1.20919677
+     1.3111856  1.26122875 0.22122563 0.55669168]]
+   #+end_example
+
+
+    #+begin_src f90 :tangle (eval f_test)
+integer(qmckl_exit_code) function test_qmckl_dist_rescaled(context) bind(C)
+
+  use qmckl
+
+  implicit none
+
+  integer(qmckl_context), intent(in), value :: context
+  integer*8                     :: m, n, LDA, LDB, LDC
+  double precision              :: x
+  integer*8                     :: i,j
+
+  double precision, parameter :: kappa = 0.6d0
+  double precision, parameter :: kappa_inv = 1.d0/kappa
+
+  integer*8, parameter :: elec_num = 10_8
+  integer*8, parameter :: nucl_num = 2_8
+
+  double precision :: nucl_coord(nucl_num,3) = reshape( (/ &
+        0.0d0,  0.0d0 ,  &
+        0.0d0,  0.0d0 ,  &
+        0.0d0,  2.059801d0 /), shape(nucl_coord) )
+
+  double precision :: elec_coord(3,elec_num) = reshape( (/ &
+        -0.250655104764153d0    ,  0.503070975550133d0    ,  -0.166554344502303d0 , &
+        -0.587812193472177d0    , -0.128751981129274d0    ,   0.187773606533075d0 , &
+         1.61335569047166d0     , -0.615556732874863d0    ,  -1.43165470979934d0  , &
+        -4.901239896295210d-003 , -1.120440036458986d-002 ,   1.99761909330422d0  , &
+         0.766647499681200d0    , -0.293515395797937d0    ,   3.66454589201239d0  , &
+        -0.127732483187947d0    , -0.138975497694196d0    ,  -8.669850480215846d-002 , &
+        -0.232271834949124d0    , -1.059321673434182d-002 ,  -0.504862241464867d0 , &
+         1.09360863531826d0     , -2.036103063808752d-003 ,  -2.702796910818986d-002 , &
+        -0.108090166832043d0    ,  0.189161729653261d0    ,   2.15398313919894d0 , &
+         0.397978144318712d0    , -0.254277292595981d0    ,   2.54553335476344d0  /), &
+        shape(elec_coord))
+
+  double precision :: ref_ee(elec_num,elec_num) = reshape( (/ &
+       0.d0, 0.63475074d0, 1.29816415d0, 1.23147027d0, 1.51933127d0, &
+       0.54402406d0, 0.51452479d0, 0.96538731d0, 1.25878564d0, 1.3722995d0 , &
+       0.63475074d0, 0.d0, 1.35148664d0, 1.13524156d0, 1.48940503d0, &
+       0.4582292d0, 0.62758076d0, 1.06560856d0, 1.179133d0, 1.30763703d0 , &
+       1.29816415d0, 1.35148664d0, 0.d0, 1.50021375d0, 1.59200788d0, &
+       1.23488312d0, 1.20844259d0, 1.0355537d0, 1.52064535d0, 1.53049239d0 , &
+       1.23147027d0, 1.13524156d0, 1.50021375d0, 0.d0, 1.12016142d0, &
+       1.19158954d0, 1.29762585d0, 1.24824277d0, 0.25292267d0, 0.58609336d0 , &
+       1.51933127d0, 1.48940503d0, 1.59200788d0, 1.12016142d0, 0.d0, &
+       1.50217017d0, 1.54012828d0, 1.48753895d0, 1.10441805d0, 0.84504381d0 , &
+       0.54402406d0, 0.4582292d0, 1.23488312d0, 1.19158954d0, 1.50217017d0, &
+       0.d0, 0.39417354d0, 0.87009603d0, 1.23838502d0, 1.33419121d0 , &
+       0.51452479d0, 0.62758076d0, 1.20844259d0, 1.29762585d0, 1.54012828d0, &
+       0.39417354d0, 0.d0, 0.95118809d0, 1.33068934d0, 1.41097406d0 , &
+       0.96538731d0, 1.06560856d0, 1.0355537d0, 1.24824277d0, 1.48753895d0, &
+       0.87009603d0, 0.95118809d0, 0.d0, 1.29422213d0, 1.33222493d0 , &
+       1.25878564d0, 1.179133d0, 1.52064535d0, 0.25292267d0, 1.10441805d0, &
+       1.23838502d0, 1.33068934d0, 1.29422213d0, 0.d0, 0.62196802d0 , &
+       1.3722995d0, 1.30763703d0, 1.53049239d0, 0.58609336d0, 0.84504381d0, &
+       1.33419121d0, 1.41097406d0, 1.33222493d0, 0.62196802d0, 0.d0 /), shape(ref_ee) )
+
+  double precision :: ref_en(elec_num, nucl_num) = reshape( (/ &
+       0.49421587d0, 0.52486545d0, 1.23280503d0, 1.16396998d0, 1.49156627d0, &
+       0.1952946d0, 0.4726453d0, 0.80211227d0, 1.21198526d0, 1.31411513d0, &
+       1.24641375d0, 1.15444238d0, 1.50565145d0, 0.06218339d0, 1.10153184d0, &
+       1.20919677d0, 1.3111856d0, 1.26122875d0, 0.22122563d0, 0.55669168d0 /), shape(ref_en) )
+
+  double precision, allocatable :: distance_ee(:,:), distance_en(:,:)
+
+  allocate( distance_ee(elec_num,elec_num), distance_en(elec_num,nucl_num) )
+
+  print *, 'ee'
+  test_qmckl_dist_rescaled = &
+       qmckl_distance_rescaled(context, 'N', 'N', elec_num, elec_num, elec_coord, &
+       size(elec_coord,1)*1_8, elec_coord, size(elec_coord,1)*1_8, &
+       distance_ee, size(distance_ee,1)*1_8, kappa)
+
+  if (test_qmckl_dist_rescaled /= QMCKL_SUCCESS) return
+
+  test_qmckl_dist_rescaled = QMCKL_FAILURE
+
+  do j=1,elec_num
+     do i=1,elec_num
+        print *, i,j,real(distance_ee(i,j)), real(ref_ee(i,j))
+        if (dabs(distance_ee(i,j) - ref_ee(i,j)) > 1.d-7) then
+           return
+        endif
+     end do
+  end do
+
+  print *, 'en'
+  test_qmckl_dist_rescaled = &
+       qmckl_distance_rescaled(context, 'N', 'T', elec_num, nucl_num, elec_coord, &
+       size(elec_coord,1)*1_8, nucl_coord, size(nucl_coord,1)*1_8, &
+       distance_en, size(distance_en,1)*1_8, kappa)
+
+  if (test_qmckl_dist_rescaled /= QMCKL_SUCCESS) return
+
+  test_qmckl_dist_rescaled = QMCKL_FAILURE
+
+  do j=1,nucl_num
+     do i=1,elec_num
+        print *, i,j,real(distance_en(i,j)), real(ref_en(i,j))
+        if (dabs(distance_en(i,j) - ref_en(i,j)) > 1.d-7) then
+           return
+        endif
+     end do
+  end do
+
+  test_qmckl_dist_rescaled = QMCKL_SUCCESS
+end function test_qmckl_dist_rescaled
+    #+end_src
+
+    #+begin_src c :comments link :tangle (eval c_test)
+    qmckl_exit_code test_qmckl_dist_rescaled(qmckl_context context);
+    assert(test_qmckl_dist_rescaled(context) == QMCKL_SUCCESS);
+    #+end_src
 * Rescaled Distance Derivatives

-** ~qmckl_distance_rescaled_deriv_e~
+** ~qmckl_distance_rescaled_gl~
   :PROPERTIES:
-   :Name:     qmckl_distance_rescaled_deriv_e
+   :Name:     qmckl_distance_rescaled_gl
   :CRetType: qmckl_exit_code
   :FRetType: qmckl_exit_code
   :END:

- ~qmckl_distance_rescaled_deriv_e~ computes the matrix of the gradient and laplacian of the
+ ~qmckl_distance_rescaled_gl~ computes the matrix of the gradient and Laplacian of the
   rescaled distance with respect to the electron coordinates. The derivative is a rank 3 tensor.
   The first dimension has a dimension of 4 of which the first three coordinates
-   contains the gradient vector and the last index is the laplacian.
+   contains the gradient vector and the last index is the Laplacian.


   \[
-   C_{ij} = \left( 1 - \exp{-\kappa C_{ij}}\right)/\kappa
+   C(r_{ij}) = \left( 1 - \exp(-\kappa\, r_{ij})\right)/\kappa
   \]

   Here the gradient is defined as follows:

   \[
-   \nabla (C_{ij}(\mathbf{r}_{ee})) = \left(\frac{\delta C_{ij}(\mathbf{r}_{ee})}{\delta x},\frac{\delta C_{ij}(\mathbf{r}_{ee})}{\delta y},\frac{\delta C_{ij}(\mathbf{r}_{ee})}{\delta z} \right)
+   \nabla_i C(r_{ij}) = \left(\frac{\partial C(r_{ij})}{\partial x_i},\frac{\partial C(r_{ij})}{\partial y_i},\frac{\partial C(r_{ij})}{\partial z_i} \right)
   \]
-   and the laplacian is defined as follows:
+   and the Laplacian is defined as follows:

   \[
-   \triangle (C_{ij}(r_{ee})) = \frac{\delta^2}{\delta x^2} + \frac{\delta^2}{\delta y^2} + \frac{\delta^2}{\delta z^2}
+   \Delta_i C(r_{ij}) = \frac{\partial^2}{\partial x_i^2} + \frac{\partial^2}{\partial y_i^2} + \frac{\partial^2}{\partial z_i^2}
   \]

-   Using the above three formulae, the expression for the gradient and laplacian is
-   as follows:
+   Using the above three formulas, the expression for the gradient and Laplacian is:

   \[
-   \frac{\delta  C_{ij}(\mathbf{r}_{ee})}{\delta x} = \frac{|(x_i - x_j)|}{r_{ij}} (1 - \kappa R_{ij})
+   \frac{\partial  C(r_{ij})}{\partial x_i} = \frac{|(x_i -
+   x_j)|}{r_{ij}} \exp (- \kappa \, r_{ij})
   \]

   \[
-   \frac{\delta  C_{ij}(\mathbf{r}_{ee})}{\delta y} = \frac{|(y_i - y_j)|}{r_{ij}} (1 - \kappa R_{ij})
-   \]
-
-   \[
-   \frac{\delta  C_{ij}(\mathbf{r}_{ee})}{\delta z} = \frac{|(z_i - z_j)|}{r_{ij}} (1 - \kappa R_{ij})
-   \]
-
-   \[
-   \Delta(C_{ij}(r_{ee}) = \left[ \frac{2}{r_{ij}} - \kappa  \right] (1-\kappa R_{ij})
+   \Delta C_{ij}(r_{ij}) = \left[ \frac{2}{r_{ij}} - \kappa  \right] \exp (- \kappa \, r_{ij})
   \]

   If the input array is normal (~'N'~), the xyz coordinates are in
   the leading dimension: ~[n][3]~ in C and ~(3,n)~ in Fortran.

-   #+NAME: qmckl_distance_rescaled_deriv_e_args
+   #+NAME: qmckl_distance_rescaled_gl_args
   | Variable               | Type                | In/Out | Description                                           |
   |------------------------+---------------------+--------+-------------------------------------------------------|
   | ~context~              | ~qmckl_context~     | in     | Global state                                          |
@ -1229,7 +1284,7 @@ end function qmckl_distance_rescaled_f
   | ~lda~                  | ~int64_t~           | in     | Leading dimension of array ~A~                        |
   | ~B~                    | ~double[][ldb]~     | in     | Array containing the $n \times 3$ matrix $B$          |
   | ~ldb~                  | ~int64_t~           | in     | Leading dimension of array ~B~                        |
-   | ~C~                    | ~double[4][n][ldc]~ | out    | Array containing the $4 \times m \times n$ matrix $C$ |
+   | ~C~                    | ~double[n][ldc][4]~ | out    | Array containing the $4 \times m \times n$ matrix $C$ |
   | ~ldc~                  | ~int64_t~           | in     | Leading dimension of array ~C~                        |
   | ~rescale_factor_kappa~ | ~double~            | in     | Factor for calculating rescaled distances derivatives |

@ -1246,12 +1301,12 @@ end function qmckl_distance_rescaled_f
    - ~A~ is allocated with at least $3 \times m \times 8$ bytes
    - ~B~ is allocated with at least $3 \times n \times 8$ bytes
    - ~C~ is allocated with at least $4 \times m \times n \times 8$ bytes
-      
-    #+CALL: generate_c_header(table=qmckl_distance_rescaled_deriv_e_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
+
+    #+CALL: generate_c_header(table=qmckl_distance_rescaled_gl_args,rettyp=get_value("CRetType"),fname=get_value("Name"))

    #+RESULTS:
    #+begin_src c :tangle (eval h_func) :comments org
-    qmckl_exit_code qmckl_distance_rescaled_deriv_e (
+    qmckl_exit_code qmckl_distance_rescaled_gl (
          const qmckl_context context,
          const char transa,
          const char transb,
@ -1263,33 +1318,35 @@ end function qmckl_distance_rescaled_f
          const int64_t ldb,
          double* const C,
          const int64_t ldc,
-          const double rescale_factor_kappa ); 
+          const double rescale_factor_kappa );
    #+end_src

    #+begin_src f90 :tangle (eval f)
-integer function qmckl_distance_rescaled_deriv_e_f(context, transa, transb, m, n, &
+function qmckl_distance_rescaled_gl(context, transa, transb, m, n, &
     A, LDA, B, LDB, C, LDC, rescale_factor_kappa) &
-     result(info)
-  use qmckl
+     bind(C) result(info)
+  use qmckl_constants
  implicit none
-  integer(qmckl_context)  , intent(in)  :: context
-  character  , intent(in)  :: transa, transb
-  integer*8  , intent(in)  :: m, n
-  integer*8  , intent(in)  :: lda
-  real*8     , intent(in)  :: A(lda,*)
-  integer*8  , intent(in)  :: ldb
-  real*8     , intent(in)  :: B(ldb,*)
-  integer*8  , intent(in)  :: ldc
-  real*8     , intent(out) :: C(4,ldc,*)
-  real*8     , intent(in)  :: rescale_factor_kappa
+
+  integer(qmckl_exit_code) :: info
+  integer (qmckl_context), intent(in)  , value :: context
+  character(c_char  ) , intent(in)  , value :: transa
+  character(c_char  ) , intent(in)  , value :: transb
+  integer (c_int64_t) , intent(in)  , value :: m
+  integer (c_int64_t) , intent(in)  , value :: n
+  integer (c_int64_t) , intent(in)  , value :: lda
+  integer (c_int64_t) , intent(in)  , value :: ldb
+  integer (c_int64_t) , intent(in)  , value :: ldc
+  real    (c_double ) , intent(in)  , value :: rescale_factor_kappa
+  real    (c_double ) , intent(in)          :: A(lda,*)
+  real    (c_double ) , intent(in)          :: B(ldb,*)
+  real    (c_double ) , intent(out)         :: C(4,ldc,n)

  integer*8 :: i,j
  real*8    :: x, y, z, dist, dist_inv
-  real*8    :: rescale_factor_kappa_inv, rij
+  real*8    :: rij
  integer   :: transab

-  rescale_factor_kappa_inv = 1.0d0/rescale_factor_kappa;
-
  info = QMCKL_SUCCESS

  if (context == QMCKL_NULL_CONTEXT) then
@ -1339,27 +1396,7 @@ integer function qmckl_distance_rescaled_deriv_e_f(context, transa, transb, m, n
     return
  endif

-  if (iand(transab,2) == 0 .and. LDA < 3) then
-     info = QMCKL_INVALID_ARG_7
-     return
-  endif
-
-  if (iand(transab,2) == 2 .and. LDA < m) then
-     info = QMCKL_INVALID_ARG_7
-     return
-  endif
-
  ! check for LDB
-  if (iand(transab,1) == 0 .and. LDB < 3) then
-     info = QMCKL_INVALID_ARG_9
-     return
-  endif
-
-  if (iand(transab,1) == 1 .and. LDB < n) then
-     info = QMCKL_INVALID_ARG_9
-     return
-  endif
-
  if (iand(transab,2) == 0 .and. LDB < 3) then
     info = QMCKL_INVALID_ARG_9
     return
@ -1386,13 +1423,13 @@ integer function qmckl_distance_rescaled_deriv_e_f(context, transa, transb, m, n
           x = A(1,i) - B(1,j)
           y = A(2,i) - B(2,j)
           z = A(3,i) - B(3,j)
-           dist = dsqrt(x*x + y*y + z*z)
+           dist = max(1.d-20, dsqrt(x*x + y*y + z*z))
           dist_inv = 1.0d0/dist
-           rij = (1.0d0 - dexp(-rescale_factor_kappa * dist)) * rescale_factor_kappa_inv
-           C(1,i,j) = x * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
-           C(2,i,j) = y * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
-           C(3,i,j) = z * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
-           C(4,i,j) = (2.0d0 * dist_inv - rescale_factor_kappa_inv) * ( 1.0d0 - rescale_factor_kappa_inv * rij)
+           rij = dexp(-rescale_factor_kappa * dist)
+           C(1,i,j) = x * dist_inv * rij
+           C(2,i,j) = y * dist_inv * rij
+           C(3,i,j) = z * dist_inv * rij
+           C(4,i,j) = (2.0d0 * dist_inv - rescale_factor_kappa) * rij
        end do
     end do

@ -1403,13 +1440,13 @@ integer function qmckl_distance_rescaled_deriv_e_f(context, transa, transb, m, n
           x = A(i,1) - B(1,j)
           y = A(i,2) - B(2,j)
           z = A(i,3) - B(3,j)
-           dist = dsqrt(x*x + y*y + z*z)
+           dist = max(1.d-20, dsqrt(x*x + y*y + z*z))
           dist_inv = 1.0d0/dist
-           rij = (1.0d0 - dexp(-rescale_factor_kappa * dist)) * rescale_factor_kappa_inv
-           C(1,i,j) = x * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
-           C(2,i,j) = y * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
-           C(3,i,j) = z * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
-           C(4,i,j) = (2.0d0 * dist_inv - rescale_factor_kappa_inv) * ( 1.0d0 - rescale_factor_kappa_inv * rij)
+           rij = dexp(-rescale_factor_kappa * dist)
+           C(1,i,j) = x * dist_inv * rij
+           C(2,i,j) = y * dist_inv * rij
+           C(3,i,j) = z * dist_inv * rij
+           C(4,i,j) = (2.0d0 * dist_inv - rescale_factor_kappa) * rij
        end do
     end do

@ -1420,13 +1457,13 @@ integer function qmckl_distance_rescaled_deriv_e_f(context, transa, transb, m, n
           x = A(1,i) - B(j,1)
           y = A(2,i) - B(j,2)
           z = A(3,i) - B(j,3)
-           dist = dsqrt(x*x + y*y + z*z)
+           dist = max(1.d-20, dsqrt(x*x + y*y + z*z))
           dist_inv = 1.0d0/dist
-           rij = (1.0d0 - dexp(-rescale_factor_kappa * dist)) * rescale_factor_kappa_inv
-           C(1,i,j) = x * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
-           C(2,i,j) = y * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
-           C(3,i,j) = z * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
-           C(4,i,j) = (2.0d0 * dist_inv - rescale_factor_kappa_inv) * ( 1.0d0 - rescale_factor_kappa_inv * rij)
+           rij = dexp(-rescale_factor_kappa * dist)
+           C(1,i,j) = x * dist_inv * rij
+           C(2,i,j) = y * dist_inv * rij
+           C(3,i,j) = z * dist_inv * rij
+           C(4,i,j) = (2.0d0 * dist_inv - rescale_factor_kappa) * rij
        end do
     end do

@ -1437,80 +1474,50 @@ integer function qmckl_distance_rescaled_deriv_e_f(context, transa, transb, m, n
           x = A(i,1) - B(j,1)
           y = A(i,2) - B(j,2)
           z = A(i,3) - B(j,3)
-           dist = dsqrt(x*x + y*y + z*z)
+           dist = max(1.d-20, dsqrt(x*x + y*y + z*z))
           dist_inv = 1.0d0/dist
-           rij = (1.0d0 - dexp(-rescale_factor_kappa * dist)) * rescale_factor_kappa_inv
-           C(1,i,j) = x * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
-           C(2,i,j) = y * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
-           C(3,i,j) = z * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
-           C(4,i,j) = (2.0d0 * dist_inv - rescale_factor_kappa_inv) * ( 1.0d0 - rescale_factor_kappa_inv * rij)
+           rij = dexp(-rescale_factor_kappa * dist)
+           C(1,i,j) = x * dist_inv * rij
+           C(2,i,j) = y * dist_inv * rij
+           C(3,i,j) = z * dist_inv * rij
+           C(4,i,j) = (2.0d0 * dist_inv - rescale_factor_kappa) * rij
        end do
     end do

  end select

-end function qmckl_distance_rescaled_deriv_e_f
+end function qmckl_distance_rescaled_gl
    #+end_src

    This function is more efficient when ~A~ and ~B~ are transposed.
-    
-    #+CALL: generate_c_interface(table=qmckl_distance_rescaled_deriv_e_args,fname=get_value("Name"))

-   #+RESULTS:
-   #+begin_src f90 :tangle (eval f) :comments org :exports none
-   integer(c_int32_t) function qmckl_distance_rescaled_deriv_e &
-       (context, transa, transb, m, n, A, lda, B, ldb, C, ldc, rescale_factor_kappa) &
-       bind(C) result(info)

-     use, intrinsic :: iso_c_binding
-     implicit none
-
-     integer (c_int64_t) , intent(in)  , value :: context
-     character           , intent(in)  , value :: transa
-     character           , intent(in)  , value :: transb
-     integer (c_int64_t) , intent(in)  , value :: m
-     integer (c_int64_t) , intent(in)  , value :: n
-     real    (c_double ) , intent(in)          :: A(lda,*)
-     integer (c_int64_t) , intent(in)  , value :: lda
-     real    (c_double ) , intent(in)          :: B(ldb,*)
-     integer (c_int64_t) , intent(in)  , value :: ldb
-     real    (c_double ) , intent(out)         :: C(ldc,n,4)
-     integer (c_int64_t) , intent(in)  , value :: ldc
-     real    (c_double ) , intent(in)  , value :: rescale_factor_kappa
-
-     integer(c_int32_t), external :: qmckl_distance_rescaled_deriv_e_f
-     info = qmckl_distance_rescaled_deriv_e_f &
-            (context, transa, transb, m, n, A, lda, B, ldb, C, ldc, rescale_factor_kappa)
-
-   end function qmckl_distance_rescaled_deriv_e
-   #+end_src
-
-   #+CALL: generate_f_interface(table=qmckl_distance_rescaled_deriv_e_args,rettyp=get_value("FRetType"),fname=get_value("Name"))
+   #+CALL: generate_f_interface(table=qmckl_distance_rescaled_gl_args,rettyp=get_value("FRetType"),fname=get_value("Name"))

   #+RESULTS:
   #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
   interface
-     integer(c_int32_t) function qmckl_distance_rescaled_deriv_e &
+     integer(qmckl_exit_code) function qmckl_distance_rescaled_gl &
         (context, transa, transb, m, n, A, lda, B, ldb, C, ldc, rescale_factor_kappa) &
         bind(C)
       use, intrinsic :: iso_c_binding
       import
       implicit none

-       integer (c_int64_t) , intent(in)  , value :: context
-       character           , intent(in)  , value :: transa
-       character           , intent(in)  , value :: transb
+       integer (qmckl_context), intent(in)  , value :: context
+       character(c_char  ) , intent(in)  , value :: transa
+       character(c_char  ) , intent(in)  , value :: transb
       integer (c_int64_t) , intent(in)  , value :: m
       integer (c_int64_t) , intent(in)  , value :: n
-       real    (c_double ) , intent(in)          :: A(lda,*)
       integer (c_int64_t) , intent(in)  , value :: lda
-       real    (c_double ) , intent(in)          :: B(ldb,*)
       integer (c_int64_t) , intent(in)  , value :: ldb
-       real    (c_double ) , intent(out)         :: C(ldc,n,4)
       integer (c_int64_t) , intent(in)  , value :: ldc
       real    (c_double ) , intent(in)  , value :: rescale_factor_kappa
+       real    (c_double ) , intent(in)          :: A(lda,*)
+       real    (c_double ) , intent(in)          :: B(ldb,*)
+       real    (c_double ) , intent(out)         :: C(4,ldc,n)

-     end function qmckl_distance_rescaled_deriv_e
+     end function qmckl_distance_rescaled_gl
   end interface
   #+end_src

--- a/org/qmckl_electron.org
+++ b/org/qmckl_electron.org
@ -14,6 +14,8 @@ the components is =[ (x,y,z), (x,y,z), ... ]=
 Using the ='T'= flage, it is =[ (x,x,x,...), (y,y,y,...), (z,z,z,...) ]=

 # TODO: replace the qmckl_matrix by qmckl_point data structures.
+# TODO: in provide funcions, replace the way to check that dimensions
+# have changed

 * Headers                                                          :noexport:
  #+begin_src elisp :noexport :results none
@ -91,16 +93,16 @@ int main() {

  Computed data:

-  | Variable                            | Type                                   | Description                                                          |
-  |-------------------------------------+----------------------------------------+----------------------------------------------------------------------|
-  | ~ee_distance~                       | ~double[walker.num][num][num]~         | Electron-electron distances                                          |
-  | ~ee_distance_date~                  | ~uint64_t~                             | Last modification date of the electron-electron distances            |
-  | ~en_distance~                       | ~double[walker.num][nucl_num][num]~    | Electron-nucleus distances                                           |
-  | ~en_distance_date~                  | ~uint64_t~                             | Last modification date of the electron-electron distances            |
-  | ~ee_potential~                      | ~double[walker.num]~                   | Electron-electron potential energy                                   |
-  | ~ee_potential_date~                 | ~uint64_t~                             | Last modification date of the electron-electron potential            |
-  | ~en_potential~                      | ~double[walker.num]~                   | Electron-nucleus potential energy                                    |
-  | ~en_potential_date~                 | ~int64_t~                              | Date when the electron-nucleus potential energy was computed         |
+  | Variable            | Type                           | Description                                                  |
+  |---------------------+--------------------------------+--------------------------------------------------------------|
+  | ~ee_distance~       | ~double[walker.num][num][num]~ | Electron-electron distances                                  |
+  | ~ee_distance_date~  | ~uint64_t~                     | Last modification date of the electron-electron distances    |
+  | ~en_distance~       | ~double[num][nucl_num]~        | Electron-nucleus distances                                   |
+  | ~en_distance_date~  | ~uint64_t~                     | Last modification date of the electron-electron distances    |
+  | ~ee_potential~      | ~double[walker.num]~           | Electron-electron potential energy                           |
+  | ~ee_potential_date~ | ~uint64_t~                     | Last modification date of the electron-electron potential    |
+  | ~en_potential~      | ~double[walker.num]~           | Electron-nucleus potential energy                            |
+  | ~en_potential_date~ | ~int64_t~                      | Date when the electron-nucleus potential energy was computed |

 ** Data structure

@ -110,7 +112,7 @@ typedef struct qmckl_walker_struct {
  int64_t num;
  qmckl_point_struct point;
 } qmckl_walker;
-  
+
 typedef struct qmckl_electron_struct {
  int64_t        num;
  int64_t        up_num;
@ -285,7 +287,7 @@ qmckl_set_electron_coord(qmckl_context context,
 {

  int32_t mask = 0;
-  
+
  <<pre2>>

  if (transp != 'N' && transp != 'T') {
@ -332,6 +334,10 @@ qmckl_set_electron_coord(qmckl_context context,
  ctx->electron.walker.num = walk_num;
  memcpy(&(ctx->electron.walker.point), &(ctx->point), sizeof(qmckl_point_struct));

+  // If it is the first time we set the electrons, we set also walkers_old.
+  if (ctx->electron.walker_old.num == 0) {
+    qmckl_set_electron_coord(context, transp, walk_num, coord, size_max);
+  }
  return QMCKL_SUCCESS;

 }
@ -345,9 +351,9 @@ interface
    implicit none

    integer (c_int64_t) , intent(in)  , value :: context
-    character           , intent(in)  , value :: transp
+    character(c_char)   , intent(in)  , value :: transp
    integer (c_int64_t) , intent(in)  , value :: walk_num
-    double precision    , intent(in)          :: coord(*)
+    real(c_double)      , intent(in)          :: coord(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function
 end interface
@ -562,6 +568,21 @@ qmckl_get_electron_coord (const qmckl_context context,
 }

      #+end_src
+
+   #+begin_src f90 :comments org :tangle (eval fh_func) :noweb yes
+interface
+  integer(c_int32_t) function qmckl_get_electron_coord(context, transp, coord, size_max) bind(C)
+    use, intrinsic :: iso_c_binding
+    import
+    implicit none
+
+    integer (c_int64_t) , intent(in)  , value :: context
+    character(c_char)   , intent(in)  , value :: transp
+    real(c_double)      , intent(in)          :: coord(*)
+    integer (c_int64_t) , intent(in)  , value :: size_max
+  end function
+end interface
+   #+end_src
 ** Test

   #+begin_src python :results output :exports none
@ -715,13 +736,13 @@ qmckl_exit_code qmckl_provide_ee_distance(qmckl_context context)


 /* Compute if necessary */
-  if (ctx->electron.walker.point.date > ctx->electron.ee_distance_date) {
+  if (ctx->point.date > ctx->electron.ee_distance_date) {

    if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
      free(ctx->electron.ee_distance);
      ctx->electron.ee_distance = NULL;
    }
-    
+
 /* Allocate array */
    if (ctx->electron.ee_distance == NULL) {

@ -965,7 +986,7 @@ qmckl_exit_code qmckl_provide_ee_potential(qmckl_context context)
  if (rc != QMCKL_SUCCESS) return rc;

 /* Compute if necessary */
-  if (ctx->electron.walker.point.date > ctx->electron.ee_potential_date) {
+  if (ctx->point.date > ctx->electron.ee_potential_date) {

    if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
      free(ctx->electron.ee_potential);
@ -1144,7 +1165,7 @@ qmckl_exit_code qmckl_get_electron_en_distance(qmckl_context context, double* di
  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
  assert (ctx != NULL);

-  size_t sze = ctx->electron.num * ctx->nucleus.num * ctx->electron.walker.num;
+  size_t sze = ctx->point.num * ctx->nucleus.num;
  memcpy(distance, ctx->electron.en_distance, sze * sizeof(double));

  return QMCKL_SUCCESS;
@ -1176,19 +1197,28 @@ qmckl_exit_code qmckl_provide_en_distance(qmckl_context context)
  }

 /* Compute if necessary */
-  if (ctx->electron.walker.point.date > ctx->electron.en_distance_date) {
+  if (ctx->point.date > ctx->electron.en_distance_date) {

-    if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
-      free(ctx->electron.en_distance);
-      ctx->electron.en_distance = NULL;
+    qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
+    mem_info.size = ctx->point.num * ctx->nucleus.num * sizeof(double);
+
+    if (ctx->electron.en_distance != NULL) {
+      qmckl_memory_info_struct mem_info_test = qmckl_memory_info_struct_zero;
+      qmckl_exit_code rc = qmckl_get_malloc_info(context, ctx->electron.en_distance, &mem_info_test);
+
+      /* if rc != QMCKL_SUCCESS, we are maybe in an _inplace function because the
+         memory was not allocated with qmckl_malloc */
+
+      if ((rc == QMCKL_SUCCESS) && (mem_info_test.size != mem_info.size)) {
+        rc = qmckl_free(context, ctx->electron.en_distance);
+        assert (rc == QMCKL_SUCCESS);
+        ctx->electron.en_distance = NULL;
+      }
    }

 /* Allocate array */
    if (ctx->electron.en_distance == NULL) {

-      qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
-      mem_info.size = ctx->electron.num * ctx->nucleus.num *
-        ctx->electron.walker.num * sizeof(double);
      double* en_distance = (double*) qmckl_malloc(context, mem_info);

      if (en_distance == NULL) {
@ -1202,10 +1232,9 @@ qmckl_exit_code qmckl_provide_en_distance(qmckl_context context)

    qmckl_exit_code rc =
      qmckl_compute_en_distance(context,
-                                ctx->electron.num,
+                                ctx->point.num,
                                ctx->nucleus.num,
-                                ctx->electron.walker.num,
-                                ctx->electron.walker.point.coord.data,
+                                ctx->point.coord.data,
                                ctx->nucleus.coord.data,
                                ctx->electron.en_distance);
    if (rc != QMCKL_SUCCESS) {
@ -1227,28 +1256,26 @@ qmckl_exit_code qmckl_provide_en_distance(qmckl_context context)
   :END:

    #+NAME: qmckl_en_distance_args
-    | Variable      | Type                                   | In/Out | Description                |
-    |---------------+----------------------------------------+--------+----------------------------|
-    | ~context~     | ~qmckl_context~                        | in     | Global state               |
-    | ~elec_num~    | ~int64_t~                              | in     | Number of electrons        |
-    | ~nucl_num~    | ~int64_t~                              | in     | Number of nuclei           |
-    | ~walk_num~    | ~int64_t~                              | in     | Number of walkers          |
-    | ~elec_coord~  | ~double[3][walk_num][elec_num]~        | in     | Electron coordinates       |
-    | ~nucl_coord~  | ~double[3][elec_num]~                  | in     | Nuclear coordinates        |
-    | ~en_distance~ | ~double[walk_num][nucl_num][elec_num]~ | out    | Electron-nucleus distances |
+    | Variable      | Type                          | In/Out | Description                |
+    |---------------+-------------------------------+--------+----------------------------|
+    | ~context~     | ~qmckl_context~               | in     | Global state               |
+    | ~point_num~   | ~int64_t~                     | in     | Number of points           |
+    | ~nucl_num~    | ~int64_t~                     | in     | Number of nuclei           |
+    | ~elec_coord~  | ~double[3][point_num]~        | in     | Electron coordinates       |
+    | ~nucl_coord~  | ~double[3][nucl_num]~         | in     | Nuclear coordinates        |
+    | ~en_distance~ | ~double[point_num][nucl_num]~ | out    | Electron-nucleus distances |

    #+begin_src f90 :comments org :tangle (eval f) :noweb yes
-integer function qmckl_compute_en_distance_f(context, elec_num, nucl_num, walk_num, elec_coord, nucl_coord, en_distance) &
+integer function qmckl_compute_en_distance_f(context, point_num, nucl_num, elec_coord, nucl_coord, en_distance) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
-  integer*8             , intent(in)  :: elec_num
+  integer*8             , intent(in)  :: point_num
  integer*8             , intent(in)  :: nucl_num
-  integer*8             , intent(in)  :: walk_num
-  double precision      , intent(in)  :: elec_coord(elec_num,walk_num,3)
+  double precision      , intent(in)  :: elec_coord(point_num,3)
  double precision      , intent(in)  :: nucl_coord(nucl_num,3)
-  double precision      , intent(out) :: en_distance(elec_num,nucl_num,walk_num)
+  double precision      , intent(out) :: en_distance(nucl_num,point_num)

  integer*8 :: k

@ -1259,7 +1286,7 @@ integer function qmckl_compute_en_distance_f(context, elec_num, nucl_num, walk_n
     return
  endif

-  if (elec_num <= 0) then
+  if (point_num <= 0) then
     info = QMCKL_INVALID_ARG_2
     return
  endif
@ -1269,20 +1296,10 @@ integer function qmckl_compute_en_distance_f(context, elec_num, nucl_num, walk_n
     return
  endif

-  if (walk_num <= 0) then
-     info = QMCKL_INVALID_ARG_4
-     return
-  endif
-
-  do k=1,walk_num
-     info = qmckl_distance(context, 'T', 'T', elec_num, nucl_num, &
-          elec_coord(1,k,1), elec_num * walk_num, &
+  info = qmckl_distance(context, 'T', 'T', nucl_num, point_num, &
          nucl_coord, nucl_num, &
-          en_distance(1,1,k), elec_num)
-     if (info /= QMCKL_SUCCESS) then
-        exit
-     endif
-  end do
+          elec_coord, point_num, &
+          en_distance, nucl_num)

 end function qmckl_compute_en_distance_f
    #+end_src
@ -1290,9 +1307,8 @@ end function qmckl_compute_en_distance_f
    #+begin_src c :tangle (eval h_private_func) :comments org :exports none
 qmckl_exit_code qmckl_compute_en_distance (
          const qmckl_context context,
-          const int64_t elec_num,
+          const int64_t point_num,
          const int64_t nucl_num,
-          const int64_t walk_num,
          const double* elec_coord,
          const double* nucl_coord,
          double* const en_distance );
@ -1303,23 +1319,22 @@ qmckl_exit_code qmckl_compute_en_distance (
    #+RESULTS:
    #+begin_src f90 :tangle (eval f) :comments org :exports none
    integer(c_int32_t) function qmckl_compute_en_distance &
-        (context, elec_num, nucl_num, walk_num, elec_coord, nucl_coord, en_distance) &
+        (context, point_num, nucl_num, elec_coord, nucl_coord, en_distance) &
        bind(C) result(info)

      use, intrinsic :: iso_c_binding
      implicit none

      integer (c_int64_t) , intent(in)  , value :: context
-      integer (c_int64_t) , intent(in)  , value :: elec_num
+      integer (c_int64_t) , intent(in)  , value :: point_num
      integer (c_int64_t) , intent(in)  , value :: nucl_num
-      integer (c_int64_t) , intent(in)  , value :: walk_num
-      real    (c_double ) , intent(in)          :: elec_coord(elec_num,walk_num,3)
-      real    (c_double ) , intent(in)          :: nucl_coord(elec_num,3)
-      real    (c_double ) , intent(out)         :: en_distance(elec_num,nucl_num,walk_num)
+      real    (c_double ) , intent(in)          :: elec_coord(point_num,3)
+      real    (c_double ) , intent(in)          :: nucl_coord(nucl_num,3)
+      real    (c_double ) , intent(out)         :: en_distance(nucl_num,point_num)

      integer(c_int32_t), external :: qmckl_compute_en_distance_f
      info = qmckl_compute_en_distance_f &
-             (context, elec_num, nucl_num, walk_num, elec_coord, nucl_coord, en_distance)
+             (context, point_num, nucl_num, elec_coord, nucl_coord, en_distance)

    end function qmckl_compute_en_distance
    #+end_src
@ -1368,7 +1383,7 @@ qmckl_check(context, rc);

 assert(qmckl_nucleus_provided(context));

-double en_distance[walk_num][nucl_num][elec_num];
+double en_distance[walk_num][elec_num][nucl_num];

 rc = qmckl_get_electron_en_distance(context, &(en_distance[0][0][0]));
 qmckl_check(context, rc);
@ -1378,19 +1393,19 @@ qmckl_check(context, rc);
 assert(fabs(en_distance[0][0][0] - 7.546738741619978) < 1.e-12);

 // (1,2,1)
-assert(fabs(en_distance[0][1][0] - 8.77102435246984) < 1.e-12);
+assert(fabs(en_distance[0][0][1] - 8.77102435246984) < 1.e-12);

 // (2,1,1)
-assert(fabs(en_distance[0][0][1] - 3.698922010513608) < 1.e-12);
+assert(fabs(en_distance[0][1][0] - 3.698922010513608) < 1.e-12);

 // (1,1,2)
 assert(fabs(en_distance[1][0][0] - 5.824059436060509) < 1.e-12);

 // (1,2,2)
-assert(fabs(en_distance[1][1][0] - 7.080482110317645) < 1.e-12);
+assert(fabs(en_distance[1][0][1] - 7.080482110317645) < 1.e-12);

 // (2,1,2)
-assert(fabs(en_distance[1][0][1] - 3.1804527583077356) < 1.e-12);
+assert(fabs(en_distance[1][1][0] - 3.1804527583077356) < 1.e-12);

     #+end_src

@ -1456,7 +1471,7 @@ qmckl_exit_code qmckl_provide_en_potential(qmckl_context context)
  if (rc != QMCKL_SUCCESS) return rc;

 /* Compute if necessary */
-  if (ctx->electron.walker.point.date > ctx->electron.en_potential_date) {
+  if (ctx->point.date > ctx->electron.en_potential_date) {

    if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
      free(ctx->electron.en_potential);
@ -1512,7 +1527,7 @@ qmckl_exit_code qmckl_provide_en_potential(qmckl_context context)
    | ~nucl_num~     | ~int64_t~                              | in     | Number of nuclei                     |
    | ~walk_num~     | ~int64_t~                              | in     | Number of walkers                    |
    | ~charge~       | ~double[nucl_num]~                     | in     | charge of nucleus                    |
-    | ~en_distance~  | ~double[walk_num][nucl_num][elec_num]~ | in     | Electron-electron rescaled distances |
+    | ~en_distance~  | ~double[walk_num][elec_num][nucl_num]~ | in     | Electron-electron distances          |
    | ~en_potential~ | ~double[walk_num]~                     | out    | Electron-electron potential          |

    #+begin_src f90 :comments org :tangle (eval f) :noweb yes
@ -1526,7 +1541,7 @@ integer function qmckl_compute_en_potential_f(context, elec_num, nucl_num, walk_
  integer*8             , intent(in)  :: nucl_num
  integer*8             , intent(in)  :: walk_num
  double precision      , intent(in)  :: charge(nucl_num)
-  double precision      , intent(in)  :: en_distance(elec_num,nucl_num,walk_num)
+  double precision      , intent(in)  :: en_distance(nucl_num,elec_num,walk_num)
  double precision      , intent(out) :: en_potential(walk_num)

  integer*8 :: nw, i, j
@ -1550,10 +1565,10 @@ integer function qmckl_compute_en_potential_f(context, elec_num, nucl_num, walk_

  en_potential = 0.0d0
  do nw=1,walk_num
-    do j=1,nucl_num
-      do i=1,elec_num
-        if (dabs(en_distance(i,j,nw)) > 1e-5) then
-          en_potential(nw) = en_potential(nw) - charge(j)/(en_distance(i,j,nw))
+    do i=1,elec_num
+      do j=1,nucl_num
+        if (dabs(en_distance(j,i,nw)) > 1.d-6) then
+          en_potential(nw) = en_potential(nw) - charge(j)/(en_distance(j,i,nw))
        endif
      end do
    end do
@ -1592,7 +1607,7 @@ end function qmckl_compute_en_potential_f
      integer (c_int64_t) , intent(in)  , value :: nucl_num
      integer (c_int64_t) , intent(in)  , value :: walk_num
      real    (c_double ) , intent(in)          :: charge(nucl_num)
-      real    (c_double ) , intent(in)          :: en_distance(elec_num,nucl_num,walk_num)
+      real    (c_double ) , intent(in)          :: en_distance(nucl_num,elec_num,walk_num)
      real    (c_double ) , intent(out)         :: en_potential(walk_num)

      integer(c_int32_t), external :: qmckl_compute_en_potential_f
--- a/org/qmckl_error.org
+++ b/org/qmckl_error.org
@ -362,7 +362,7 @@ void qmckl_string_of_error_f(const qmckl_exit_code error, char result[<<MAX_STRI
       import
       implicit none
       integer (qmckl_exit_code), intent(in), value :: error
-       character, intent(out) :: string(<<MAX_STRING_LENGTH()>>)
+       character(c_char), intent(out) :: string(<<MAX_STRING_LENGTH()>>)
     end subroutine qmckl_string_of_error
  end interface
  #+end_src
@ -440,7 +440,7 @@ qmckl_set_error(qmckl_context context,

  Upon error, the error type and message can be obtained from the
  context using ~qmckl_get_error~. The message and function name
-  is returned in the variables provided. Therefore, passing a 
+  is returned in the variables provided. Therefore, passing a
  function name and message is mandatory.

  # Header
@ -494,7 +494,7 @@ qmckl_get_error(qmckl_context context,
    #+end_src

 * Failing
-  
+
   To make a function fail, the ~qmckl_failwith~ function should be
   called, such that information about the failure is stored in
   the context. The desired exit code is given as an argument, as
@ -511,7 +511,7 @@ qmckl_failwith(qmckl_context context,
               const char* function,
               const char* message) ;
   #+end_src
-   
+
   #+begin_src c :comments org :tangle (eval c) :exports none
 qmckl_exit_code
 qmckl_failwith(qmckl_context context,
@ -593,7 +593,7 @@ qmckl_last_error(qmckl_context context, char* buffer) {
  return QMCKL_SUCCESS;
 }
    #+end_src
-    
+
 ** Fortran inteface

   #+begin_src f90 :tangle (eval fh_func) :exports none :noweb yes
@ -603,7 +603,7 @@ qmckl_last_error(qmckl_context context, char* buffer) {
       import
       implicit none
       integer (c_int64_t) , intent(in), value :: context
-       character, intent(out) :: string(*)
+       character(c_char), intent(out) :: string(*)
     end subroutine qmckl_last_error
  end interface
   #+end_src
@ -625,7 +625,7 @@ qmckl_check(qmckl_context context, qmckl_exit_code rc);

 qmckl_exit_code
 qmckl_check(qmckl_context context, qmckl_exit_code rc)
-{  
+{

  if (rc != QMCKL_SUCCESS) {
    char fname[QMCKL_MAX_FUN_LEN];
@ -639,7 +639,7 @@ qmckl_check(qmckl_context context, qmckl_exit_code rc)
  return rc;
 }
    #+end_src
-    
+
  It should be used as:

    #+begin_src c
@ -648,7 +648,7 @@ rc = qmckl_check(context,
                 );
 assert (rc == QMCKL_SUCCESS);
     #+end_src
-    
+
 ** Fortran inteface

   #+begin_src f90 :tangle (eval fh_func) :exports none :noweb yes
--- a/org/qmckl_examples.org
+++ b/org/qmckl_examples.org
@ -1,28 +1,72 @@
 #+TITLE: Code examples
 #+SETUPFILE: ../tools/theme.setup
 #+INCLUDE: ../tools/lib.org
-  
-In this section, we present examples of usage of QMCkl.
-For simplicity, we assume that the wave function parameters are stored
-in a [[https://github.com/TREX-CoE/trexio][TREXIO]] file.

-* Python
-** Check numerically that MOs are orthonormal
+In this section, we provide hands-on examples to demonstrate the capabilities
+of the QMCkl library.  We furnish code samples in C, Fortran, and Python,
+serving as exhaustive tutorials for effectively leveraging QMCkl.
+For simplicity, we assume that the wave function parameters are stored in a
+[[https://github.com/TREX-CoE/trexio][TREXIO]] file.
+
+* Overlap matrix in the MO basis
+
+The focal point of this example is the numerical evaluation of the overlap
+matrix in the MO basis. Utilizing QMCkl, we approximate the MOs at
+discrete grid points to compute the overlap matrix \( S_{ij} \) as follows:
+\[
+S_{ij} = \int \phi_i(\mathbf{r})\, \phi_j(\mathbf{r}) \text{d}\mathbf{r} \approx
+         \sum_k \phi_i(\mathbf{r}_k)\, \phi_j(\mathbf{r}_k) \delta\mathbf{r}
+\]
+
+
+The code starts by reading a wave function from a TREXIO file.  This is
+accomplished using the ~qmckl_trexio_read~ function, which populates a
+~qmckl_context~ with the necessary wave function parameters.  The context
+serves as the primary interface for interacting with the QMCkl library,
+encapsulating the state and configurations of the system.
+Subsequently, the code retrieves various attributes such as the number of
+nuclei ~nucl_num~ and coordinates ~nucl_coord~.
+These attributes are essential for setting up the integration grid.
+
+The core of the example lies in the numerical computation of the overlap
+matrix. To achieve this, the code employs a regular grid in three-dimensional
+space, and the grid points are then populated into the ~qmckl_context~ using
+the ~qmckl_set_point~ function.
+
+The MO values at these grid points are computed using the
+~qmckl_get_mo_basis_mo_value~ function.  These values are then used to
+calculate the overlap matrix through a matrix multiplication operation
+facilitated by the ~qmckl_dgemm~ function.
+
+The code is also instrumented to measure the execution time for the MO
+value computation, providing an empirical assessment of the computational
+efficiency. Error handling is robustly implemented at each stage to ensure the
+reliability of the simulation.
+
+In summary, this example serves as a holistic guide for leveraging the QMCkl
+library, demonstrating its ease of use. It provides a concrete starting point
+for researchers and developers interested in integrating QMCkl into their QMC
+code.
+
+** Python
+   :PROPERTIES:
+   :header-args: :tangle mo_ortho.py
+   :END:

   In this example, we will compute numerically the overlap
   between the molecular orbitals:
-   
+
   \[
      S_{ij} = \int \phi_i(\mathbf{r}) \phi_j(\mathbf{r})
   \text{d}\mathbf{r} \sim \sum_{k=1}^{N} \phi_i(\mathbf{r}_k)
-   \phi_j(\mathbf{r}_k) \delta \mathbf{r} 
+   \phi_j(\mathbf{r}_k) \delta \mathbf{r}
   \]
   \[
-      S_{ij} = \langle \phi_i | \phi_j \rangle 
+      S_{ij} = \langle \phi_i | \phi_j \rangle
   \sim \sum_{k=1}^{N} \langle \phi_i | \mathbf{r}_k \rangle
   \langle \mathbf{r}_k | \phi_j \rangle
   \]
-   
+

   #+begin_src python :exports code
 import numpy as np
@ -30,11 +74,11 @@ import qmckl
   #+end_src

   #+RESULTS:
-   
+
   First, we create a context for the QMCkl calculation, and load the
   wave function stored in =h2o_5z.h5= inside it. It is a Hartree-Fock
   determinant for the water molecule in the cc-pV5Z basis set.
-   
+
   #+begin_src python :exports code
 trexio_filename = "..//share/qmckl/test_data/h2o_5z.h5"

@ -49,7 +93,7 @@ qmckl.trexio_read(context, trexio_filename)
   molecule.

   We fetch the nuclear coordinates from the context,
-   
+
   #+begin_src python :exports code
 nucl_num = qmckl.get_nucleus_num(context)

@ -72,7 +116,7 @@ for i in range(nucl_num):
   #+end_example

   and compute the coordinates of the grid points:
-   
+
   #+begin_src python :exports code
 nx = ( 120, 120, 120 )
 shift = np.array([5.,5.,5.])
@ -94,7 +138,7 @@ dr = step[0] * step[1] * step[2]
   #+end_src

   #+RESULTS:
-   
+
   Now the grid is ready, we can create the list of grid points
   $\mathbf{r}_k$ on which the MOs $\phi_i$ will be evaluated, and
   transfer them to the QMCkl context:
@ -113,7 +157,7 @@ qmckl.set_point(context, 'N', point_num, np.reshape(point, (point_num*3)))

   #+RESULTS:
   : None
-   
+
   Then, we evaluate all the MOs at the grid points (and time the execution),
   and thus obtain the matrix $M_{ki} = \langle \mathbf{r}_k | \phi_i \rangle =
   \phi_i(\mathbf{r}_k)$.
@ -127,7 +171,7 @@ before   = time.time()
 mo_value = qmckl.get_mo_basis_mo_value(context, point_num*mo_num)
 after    = time.time()

-mo_value = np.reshape( mo_value, (point_num, mo_num) )
+mo_value = np.reshape( mo_value, (point_num, mo_num) ).T   # Transpose to get mo_num x point_num

 print("Number of MOs: ", mo_num)
 print("Number of grid points: ", point_num)
@ -138,14 +182,14 @@ print("Execution time : ", (after - before), "seconds")
   #+begin_example
 Number of MOs:  201
 Number of grid points:  1728000
-Execution time :  3.511528968811035 seconds
+Execution time :  5.577778577804565 seconds
   #+end_example

   and finally we compute the overlap between all the MOs as
-   $M^\dagger M$.
+   $M.M^\dagger$.

   #+begin_src python :exports code
-overlap = mo_value.T @ mo_value * dr
+overlap = mo_value @ mo_value.T * dr
 print (overlap)
   #+end_src

@ -165,6 +209,477 @@ print (overlap)
      1.18264754e-09  8.97215950e-01]]
   #+end_example

+** C
+   In this example, electron-nucleus cusp fitting is added.
+
+   :PROPERTIES:
+   :header-args: :tangle mo_ortho.c
+   :END:
+
+   In this example, we will compute numerically the overlap
+   between the molecular orbitals:
+
+   \[
+      S_{ij} = \int \phi_i(\mathbf{r}) \phi_j(\mathbf{r})
+   \text{d}\mathbf{r} \sim \sum_{k=1}^{N} \phi_i(\mathbf{r}_k)
+   \phi_j(\mathbf{r}_k) \delta \mathbf{r}
+   \]
+   \[
+      S_{ij} = \langle \phi_i | \phi_j \rangle
+   \sim \sum_{k=1}^{N} \langle \phi_i | \mathbf{r}_k \rangle
+   \langle \mathbf{r}_k | \phi_j \rangle
+   \]
+
+   We apply the cusp fitting procedure, so the MOs might deviate
+   slightly from orthonormality.
+
+   #+begin_src c :exports code
+#include <qmckl.h>
+#include <stdio.h>
+#include <string.h>
+#include <sys/time.h>
+
+int main(int argc, char** argv)
+{
+  const char* trexio_filename = "..//share/qmckl/test_data/h2o_5z.h5";
+  qmckl_exit_code rc = QMCKL_SUCCESS;
+   #+end_src
+
+   First, we create a context for the QMCkl calculation, and load the
+   wave function stored in =h2o_5z.h5= inside it. It is a Hartree-Fock
+   determinant for the water molecule in the cc-pV5Z basis set.
+
+   #+begin_src c :exports code
+  qmckl_context context = qmckl_context_create();
+
+  rc = qmckl_trexio_read(context, trexio_filename, strlen(trexio_filename));
+
+  if (rc != QMCKL_SUCCESS) {
+    fprintf(stderr, "Error reading TREXIO file:\n%s\n", qmckl_string_of_error(rc));
+    exit(1);
+  }
+
+   #+end_src
+
+   We impose the electron-nucleus cusp fitting to occur when the
+   electrons are close to the nuclei. The critical distance
+   is 0.5 atomic units for hydrogens and 0.1 for the oxygen.
+   To identify which atom is an oxygen and which are hydrogens, we
+   fetch the nuclear charges from the context.
+
+   #+begin_src c :exports code
+  int64_t nucl_num;
+
+  rc = qmckl_get_nucleus_num(context, &nucl_num);
+
+  if (rc != QMCKL_SUCCESS) {
+    fprintf(stderr, "Error getting nucl_num:\n%s\n", qmckl_string_of_error(rc));
+    exit(1);
+  }
+
+
+  double nucl_charge[nucl_num];
+
+  rc = qmckl_get_nucleus_charge(context, &(nucl_charge[0]), nucl_num);
+
+  if (rc != QMCKL_SUCCESS) {
+    fprintf(stderr, "Error getting nucl_charge:\n%s\n", qmckl_string_of_error(rc));
+    exit(1);
+  }
+
+
+  double r_cusp[nucl_num];
+
+  for (size_t i=0 ; i<nucl_num ; ++i) {
+
+    switch ((int) nucl_charge[i]) {
+
+    case 1:
+      r_cusp[i] = 0.5;
+      break;
+
+    case 8:
+      r_cusp[i] = 0.1;
+      break;
+    }
+
+  }
+
+
+  rc = qmckl_set_mo_basis_r_cusp(context, &(r_cusp[0]), nucl_num);
+
+  if (rc != QMCKL_SUCCESS) {
+    fprintf(stderr, "Error setting r_cusp:\n%s\n", qmckl_string_of_error(rc));
+    exit(1);
+  }
+
+
+   #+end_src
+
+
+   We now define the grid points $\mathbf{r}_k$ as a regular grid around the
+   molecule.
+   We fetch the nuclear coordinates from the context,
+
+   #+begin_src c :exports code
+  double nucl_coord[nucl_num][3];
+
+  rc = qmckl_get_nucleus_coord(context, 'N', &(nucl_coord[0][0]), nucl_num*3);
+
+  if (rc != QMCKL_SUCCESS) {
+    fprintf(stderr, "Error getting nucl_coord:\n%s\n", qmckl_string_of_error(rc));
+    exit(1);
+  }
+
+  for (size_t i=0 ; i<nucl_num ; ++i) {
+    printf("%d  %+f %+f %+f\n",
+           (int32_t) nucl_charge[i],
+           nucl_coord[i][0],
+           nucl_coord[i][1],
+           nucl_coord[i][2]);
+  }
+   #+end_src
+
+   #+begin_example
+8  +0.000000 +0.000000 +0.000000
+1  -1.430429 +0.000000 -1.107157
+1  +1.430429 +0.000000 -1.107157
+   #+end_example
+
+   and compute the coordinates of the grid points:
+
+   #+begin_src c :exports code
+  size_t nx[3] = { 120, 120, 120 };
+  double shift[3] = {5.,5.,5.};
+  int64_t point_num = nx[0] * nx[1] * nx[2];
+
+  double rmin[3] = { nucl_coord[0][0], nucl_coord[0][1], nucl_coord[0][2] } ;
+  double rmax[3] = { nucl_coord[0][0], nucl_coord[0][1], nucl_coord[0][2] } ;
+
+  for (size_t i=0 ; i<nucl_num ; ++i) {
+    for (int j=0 ; j<3 ; ++j) {
+      rmin[j] = nucl_coord[i][j] < rmin[j] ? nucl_coord[i][j] : rmin[j];
+      rmax[j] = nucl_coord[i][j] > rmax[j] ? nucl_coord[i][j] : rmax[j];
+    }
+  }
+
+  rmin[0] -= shift[0]; rmin[1] -= shift[1]; rmin[2] -= shift[2];
+  rmax[0] += shift[0]; rmax[1] += shift[1]; rmax[2] += shift[2];
+
+  double step[3];
+
+  double* linspace[3];
+  for (int i=0 ; i<3 ; ++i) {
+
+    linspace[i] = (double*) calloc( sizeof(double), nx[i] );
+
+    if (linspace[i] == NULL) {
+      fprintf(stderr, "Allocation failed (linspace)\n");
+      exit(1);
+    }
+
+    step[i] = (rmax[i] - rmin[i]) / ((double) (nx[i]-1));
+
+    for (size_t j=0 ; j<nx[i] ; ++j) {
+      linspace[i][j] = rmin[i] + j*step[i];
+    }
+
+  }
+
+  double dr = step[0] * step[1] * step[2];
+   #+end_src
+
+   Now the grid is ready, we can create the list of grid points
+   $\mathbf{r}_k$ on which the MOs $\phi_i$ will be evaluated, and
+   transfer them to the QMCkl context:
+
+   #+begin_src c :exports code
+  double* point = (double*) calloc(sizeof(double), 3*point_num);
+
+  if (point == NULL) {
+    fprintf(stderr, "Allocation failed (point)\n");
+    exit(1);
+  }
+
+  size_t m = 0;
+  for (size_t i=0 ; i<nx[0] ; ++i) {
+    for (size_t j=0 ; j<nx[1] ; ++j) {
+      for (size_t k=0 ; k<nx[2] ; ++k) {
+
+        point[m] = linspace[0][i];
+        m++;
+
+        point[m] = linspace[1][j];
+        m++;
+
+        point[m] = linspace[2][k];
+        m++;
+
+      }
+    }
+  }
+
+  rc = qmckl_set_point(context, 'N', point_num, point, (point_num*3));
+
+  if (rc != QMCKL_SUCCESS) {
+    fprintf(stderr, "Error setting points:\n%s\n", qmckl_string_of_error(rc));
+    exit(1);
+  }
+   #+end_src
+
+   #+RESULTS:
+   : None
+
+   Then, we evaluate all the MOs at the grid points (and time the execution),
+   and thus obtain the matrix $M_{ki} = \langle \mathbf{r}_k | \phi_i
+   \rangle = \phi_i(\mathbf{r}_k)$.
+
+   #+begin_src c :exports code
+
+  int64_t mo_num;
+  rc = qmckl_get_mo_basis_mo_num(context, &mo_num);
+
+  long before, after;
+  struct timeval timecheck;
+
+  double* mo_value = (double*) calloc(sizeof(double), point_num*mo_num);
+  if (mo_value == NULL) {
+    fprintf(stderr, "Allocation failed (mo_value)\n");
+    exit(1);
+  }
+
+  gettimeofday(&timecheck, NULL);
+  before = (long)timecheck.tv_sec * 1000 + (long)timecheck.tv_usec / 1000;
+
+  rc = qmckl_get_mo_basis_mo_value(context, mo_value, point_num*mo_num);
+  if (rc != QMCKL_SUCCESS) {
+    fprintf(stderr, "Error getting mo_value)\n");
+    exit(1);
+  }
+
+  gettimeofday(&timecheck, NULL);
+  after = (long)timecheck.tv_sec * 1000 + (long)timecheck.tv_usec / 1000;
+
+  printf("Number of MOs: %ld\n", (long) mo_num);
+  printf("Number of grid points: %ld\n", (long) point_num);
+  printf("Execution time : %f seconds\n", (after - before)*1.e-3);
+
+   #+end_src
+
+   #+begin_example
+Number of MOs:  201
+Number of grid points:  1728000
+Execution time :  5.608000 seconds
+   #+end_example
+
+   and finally we compute the overlap between all the MOs as
+   $M.M^\dagger$.
+
+   #+begin_src c :exports code
+  double* overlap = (double*) malloc (mo_num*mo_num*sizeof(double));
+
+  rc = qmckl_dgemm(context, 'N', 'T', mo_num, mo_num, point_num, dr,
+                   mo_value, mo_num, mo_value, mo_num, 0.0,
+                   overlap, mo_num);
+
+  for (size_t i=0 ; i<mo_num ; ++i) {
+    printf("%4ld", i);
+    for (size_t j=0 ; j<mo_num ; ++j) {
+      printf(" %f", overlap[i*mo_num+j]);
+    }
+    printf("\n");
+  }
+
+  // Clean-up and exit
+  free(mo_value);
+  free(overlap);
+
+  rc = qmckl_context_destroy(context);
+  if (rc != QMCKL_SUCCESS) {
+    fprintf(stderr, "Error destroying context)\n");
+    exit(1);
+  }
+
+  return 0;
+}
+   #+end_src
+
+   #+begin_example
+   0 0.988765 0.002336 0.000000 -0.000734 0.000000 0.000530 0.000000 0.000446 0.000000 -0.000000 -0.000511 -0.000000 -0.000267 0.000000 0.000000 0.001007 0.000000 0.000168 -0.000000 -0.000000 -0.000670 -0.000000 0.000000 -0.000251 -0.000261 -0.000000 -0.000000 -0.000000 -0.000397 -0.000000 -0.000810 0.000000 0.000231 -0.000000 -0.000000 0.000000 -0.000000
+   ...
+ 200 0.039017 -0.013066 -0.000000 -0.001935 -0.000000 -0.003829 -0.000000 0.000996 -0.000000 0.000000 -0.003733 0.000000 0.000065 -0.000000 -0.000000 -0.002220 -0.000000 -0.001961 0.000000 0.000000 -0.004182 0.000000 -0.000000 -0.000165 -0.002445 0.000000 -0.000000 0.000000 0.001985 0.000000 0.001685 -0.000000 -0.002899 0.000000 0.000000 0.000000 -0.000000 0.002591 0.000000 -0.000000 0.000000 0.002385 0.000000 -0.011086 0.000000 -0.003885 0.000000 -0.000000 0.003602 -0.000000 0.000000 -0.003241 0.000000 0.000000 0.002613 -0.007344 -0.000000 -0.000000 0.000000 0.000017 0.000000 0.000000 0.000000 -0.008719 0.000000 -0.001358 -0.003233 0.000000 -0.000000 -0.000000 -0.000000 0.000000 0.003778 0.000000 0.000000 -0.000000 0.000000 0.000000 -0.001190 0.000000 0.000000 -0.000000 0.005578 -0.000000 -0.001502 0.003899 -0.000000 -0.000000 0.000000 -0.003291 -0.001775 -0.000000 -0.002374 0.000000 -0.000000 -0.000000 -0.000000 -0.001290 -0.000000 0.002178 0.000000 0.000000 0.000000 -0.001252 0.000000 -0.000000 -0.000926 0.000000 -0.000000 -0.013130 -0.000000 0.012124 0.000000 -0.000000 -0.000000 -0.000000 0.000000 0.025194 0.000343 -0.000000 0.000000 -0.000000 -0.004421 0.000000 0.000000 -0.000599 -0.000000 0.005289 0.000000 -0.000000 0.012826 -0.000000 0.000000 0.006190 0.000000 0.000000 -0.000000 0.000000 -0.000321 0.000000 -0.000000 -0.000000 0.000000 -0.000000 0.001499 -0.006629 0.000000 0.000000 0.000000 -0.000000 0.008737 -0.000000 0.006880 0.000000 -0.001851 -0.000000 -0.000000 0.000000 -0.007464 0.000000 0.010298 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 0.000000 0.000540 0.000000 -0.006616 -0.000000 0.000000 -0.002927 -0.000000 0.000000 0.010352 0.000000 -0.003103 -0.000000 -0.007640 -0.000000 -0.000000 0.005302 0.000000 0.000000 -0.000000 -0.000000 -0.010181 0.000000 -0.001108 0.000000 0.000000 -0.000000 0.000000 0.000000 -0.000998 -0.009754 0.013562 0.000000 -0.000000 0.887510
+   #+end_example
+
+   
+** Fortran
+   Here is the same piece of code translated in Fortran
+   #+begin_src f90
+#include <qmckl_f.F90>
+
+program main
+  use iso_c_binding
+  use qmckl
+  implicit none
+
+  ! Declare variables
+  integer :: argc
+  character(128) :: trexio_filename, err_msg
+  integer(qmckl_exit_code) :: rc
+  integer(qmckl_context) :: context
+  integer*8 :: nucl_num, mo_num, point_num
+  double precision, allocatable :: nucl_coord(:,:)
+  integer*8  :: nx(3)
+  double precision, dimension(3) :: shift, step, rmin, rmax
+  double precision, allocatable :: mo_value(:,:), overlap(:,:), point(:), linspace(:,:)
+  double precision :: before, after, dr
+  integer*8 :: i,j,k,m
+
+  ! Initialize variables
+  err_msg = ' ' 
+  argc = command_argument_count()
+  if (argc /= 1) then
+    print *, "Usage: ./program <TREXIO filename>"
+    stop -1
+  end if
+  call get_command_argument(1, trexio_filename)
+  rc = QMCKL_SUCCESS
+
+  ! Create a QMCkl context
+  context = qmckl_context_create()
+
+  ! Read the TREXIO file into the context
+  rc = qmckl_trexio_read(context, trim(trexio_filename), len(trexio_filename)*1_8)
+  if (rc /= QMCKL_SUCCESS) then
+    call qmckl_string_of_error(rc, err_msg)
+    write(*,*) "Error reading TREXIO file:", err_msg
+    stop -1
+  end if
+
+  ! Retrieve the number of nuclei
+  rc = qmckl_get_nucleus_num(context, nucl_num)
+  if (rc /= QMCKL_SUCCESS) then
+    call qmckl_string_of_error(rc, err_msg)
+    write(*,*) "Error getting nucl_num:", err_msg
+    stop -1
+  end if
+
+  ! Retrieve the nuclear coordinates
+  allocate(nucl_coord(3, nucl_num))
+  rc = qmckl_get_nucleus_coord(context, 'N', nucl_coord, nucl_num * 3_8)
+  if (rc /= QMCKL_SUCCESS) then
+    call qmckl_string_of_error(rc, err_msg)
+    write(*,*) "Error getting nucl_coord:", err_msg
+    stop -1
+  end if
+
+  ! Retrieve the number of MOs
+  rc = qmckl_get_mo_basis_mo_num(context, mo_num)
+  if (rc /= QMCKL_SUCCESS) then
+    call qmckl_string_of_error(rc, err_msg)
+    write(*,*) "Error getting mo_num:", err_msg
+    stop -1
+  end if
+
+  ! Initialize grid points for the calculation
+  nx = (/ 120, 120, 120 /)
+  shift = (/ 5.d0, 5.d0, 5.d0 /)
+  point_num = nx(1) * nx(2) * nx(3)
+
+  ! Initialize rmin and rmax
+  rmin = nucl_coord(:,1)
+  rmax = nucl_coord(:,1)
+
+  ! Update rmin and rmax based on nucl_coord
+  do i = 1, 3
+    do j = 1, nucl_num
+      rmin(i) = min(nucl_coord(i,j), rmin(i))
+      rmax(i) = max(nucl_coord(i,j), rmax(i))
+    end do
+  end do
+
+  ! Apply shift
+  rmin = rmin - shift
+  rmax = rmax + shift
+
+  ! Initialize linspace and step
+  allocate(linspace(3, maxval(nx)))
+
+  do i = 1, 3
+    step(i) = (rmax(i) - rmin(i)) / real(nx(i) - 1, 8)
+    do j = 1, nx(i)
+      linspace(i, j) = rmin(i) + (j - 1) * step(i)
+    end do
+  end do
+
+  ! Initialize point array
+  allocate(point(3 * point_num))
+  m = 1
+  do i = 1, nx(1)
+    do j = 1, nx(2)
+      do k = 1, nx(3)
+        point(m) = linspace(1, i); m = m + 1
+        point(m) = linspace(2, j); m = m + 1
+        point(m) = linspace(3, k); m = m + 1
+      end do
+    end do
+  end do
+
+  deallocate(linspace)
+
+
+  ! Set points in QMCKL context
+  rc = qmckl_set_point(context, 'N', point_num, point, point_num * 3)
+  if (rc /= QMCKL_SUCCESS) then
+    call qmckl_string_of_error(rc, err_msg)
+    write(*,*) "Error setting point:", err_msg
+    stop -1
+  end if
+
+
+
+
+  ! Perform the actual calculation and measure the time taken
+  call cpu_time(before)
+  allocate(mo_value(point_num, mo_num))
+  rc = qmckl_get_mo_basis_mo_value(context, mo_value, point_num * mo_num)
+  if (rc /= QMCKL_SUCCESS) then
+    call qmckl_string_of_error(rc, err_msg)
+    write(*,*) "Error getting mo_value:", err_msg
+    stop
+  end if
+  call cpu_time(after)
+
+  write(*,*) "Number of MOs:", mo_num
+  write(*,*) "Number of grid points:", point_num
+  write(*,*) "Execution time:", (after - before), "seconds"
+
+  ! Compute the overlap matrix
+  dr = step(1) * step(2) * step(3)
+
+  allocate(overlap(mo_num, mo_num))
+  rc = qmckl_dgemm(context, 'N', 'T', mo_num, mo_num, point_num, dr, &
+                   mo_value, mo_num, mo_value, mo_num, 0.d0, overlap, mo_num)
+
+  ! Print the overlap matrix
+  do i = 1, mo_num
+    write(*,'(i4)', advance='no') i
+    do j = 1, mo_num
+      write(*,'(f8.4)', advance='no') overlap(i, j)
+    end do
+    write(*,*)
+  end do
+
+  ! Clean-up and exit
+  deallocate(mo_value, overlap)
+  rc = qmckl_context_destroy(context)
+  if (rc /= QMCKL_SUCCESS) then
+    call qmckl_string_of_error(rc, err_msg)
+    write(*,*) "Error destroying context:", err_msg
+    stop -1
+  end if
+
+end program main
+   #+end_src
 * Fortran
 ** Checking errors

@ -173,7 +688,7 @@ print (overlap)
   error in text format and exits the program.

   #+NAME: qmckl_check_error
-   #+begin_src f90 
+   #+begin_src f90
 subroutine qmckl_check_error(rc, message)
  use qmckl
  implicit none
@ -188,7 +703,7 @@ subroutine qmckl_check_error(rc, message)
  end if
 end subroutine qmckl_check_error
   #+end_src
-  
+
 ** Computing an atomic orbital on a grid
   :PROPERTIES:
   :header-args: :tangle ao_grid.f90
@ -200,14 +715,14 @@ end subroutine qmckl_check_error
   atomic units in the borders.

   This program uses the ~qmckl_check_error~ function defined above.
-  
+
   To use this program, run
-  
+
   #+begin_src bash :tangle no :exports code
 $ ao_grid <trexio_file> <AO_id> <point_num>
   #+end_src

-  
+
   #+begin_src f90  :noweb yes
 <<qmckl_check_error>>

@ -237,7 +752,7 @@ program ao_grid

   Start by fetching the command-line arguments:

-   #+begin_src f90 
+   #+begin_src f90
  if (iargc() /= 3) then
     print *, 'Syntax: ao_grid <trexio_file> <AO_id> <point_num>'
     call exit(-1)
@ -257,7 +772,7 @@ program ao_grid
   Create the QMCkl context and initialize it with the wave function
   present in the TREXIO file:

-   #+begin_src f90 
+   #+begin_src f90
  qmckl_ctx = qmckl_context_create()
  rc  = qmckl_trexio_read(qmckl_ctx, trexio_filename, 1_8*len(trim(trexio_filename)))
  call qmckl_check_error(rc, 'Read TREXIO')
@ -265,8 +780,8 @@ program ao_grid

   We need to check that ~ao_id~ is in the range, so we get the total
   number of AOs from QMCkl:
-  
-   #+begin_src f90 
+
+   #+begin_src f90
  rc = qmckl_get_ao_basis_ao_num(qmckl_ctx, ao_num)
  call qmckl_check_error(rc, 'Getting ao_num')

@ -279,7 +794,7 @@ program ao_grid
   Now we will compute the limits of the box in which the molecule fits.
   For that, we first need to ask QMCkl the coordinates of nuclei.

-   #+begin_src f90 
+   #+begin_src f90
  rc = qmckl_get_nucleus_num(qmckl_ctx, nucl_num)
  call qmckl_check_error(rc, 'Get nucleus num')

@ -291,11 +806,11 @@ program ao_grid
   We now compute the coordinates of opposite points of the box, and
   the distance between points along the 3 directions:

-   #+begin_src f90 
+   #+begin_src f90
  rmin(1) = minval( nucl_coord(1,:) ) - 5.d0
  rmin(2) = minval( nucl_coord(2,:) ) - 5.d0
  rmin(3) = minval( nucl_coord(3,:) ) - 5.d0
-     
+
  rmax(1) = maxval( nucl_coord(1,:) ) + 5.d0
  rmax(2) = maxval( nucl_coord(2,:) ) + 5.d0
  rmax(3) = maxval( nucl_coord(3,:) ) + 5.d0
@ -305,8 +820,8 @@ program ao_grid

   We now produce the list of point coordinates where the AO will be
   evaluated:
-  
-   #+begin_src f90 
+
+   #+begin_src f90
  point_num = point_num_x**3
  allocate( points(point_num, 3) )
  ipoint=0
@ -327,19 +842,19 @@ program ao_grid
     z = z + dr(3)
  end do
   #+end_src
-  
+
   We give the points to QMCkl:

-   #+begin_src f90 
+   #+begin_src f90
  rc = qmckl_set_point(qmckl_ctx, 'T', point_num, points, size(points)*1_8 )
  call qmckl_check_error(rc, 'Setting points')
   #+end_src

   We allocate the space required to retrieve the values, gradients and
   Laplacian of all AOs, and ask to retrieve the values of the
-   AOs computed at the point positions. 
-  
-   #+begin_src f90 
+   AOs computed at the point positions.
+
+   #+begin_src f90
  allocate( ao_vgl(ao_num, 5, point_num) )
  rc = qmckl_get_ao_basis_ao_vgl(qmckl_ctx, ao_vgl, ao_num*5_8*point_num)
  call qmckl_check_error(rc, 'Setting points')
@ -347,13 +862,13 @@ program ao_grid

   We finally print the value and Laplacian of the AO:

-   #+begin_src f90 
+   #+begin_src f90
  do ipoint=1, point_num
     print '(3(F10.6,X),2(E20.10,X))', points(ipoint, 1:3), ao_vgl(ao_id,1,ipoint), ao_vgl(ao_id,5,ipoint)
  end do
   #+end_src

-   #+begin_src f90 
+   #+begin_src f90
  deallocate( nucl_coord, points, ao_vgl )
 end program ao_grid
   #+end_src
--- a/org/qmckl_jastrow.org
+++ b/org/qmckl_jastrow.org
--- a/org/qmckl_jastrow_champ.org
+++ b/org/qmckl_jastrow_champ.org
--- a/org/qmckl_local_energy.org
+++ b/org/qmckl_local_energy.org
@ -697,7 +697,7 @@ assert (rc == QMCKL_SUCCESS);

 const double  * mo_coefficient          =  &(chbrclf_mo_coef[0]);

-rc = qmckl_set_mo_basis_coefficient(context, mo_coefficient);
+rc = qmckl_set_mo_basis_coefficient(context, mo_coefficient, chbrclf_ao_num*chbrclf_mo_num);
 assert (rc == QMCKL_SUCCESS);

 assert(qmckl_mo_basis_provided(context));
--- a/org/qmckl_memory.org
+++ b/org/qmckl_memory.org
@ -104,12 +104,46 @@ typedef struct qmckl_memory_struct {
  properly freed when the library is de-initialized.
  If the allocation failed, the ~NULL~ pointer is returned.

+  The allocated memory block is zeroed using ~memset~.
+
  # Header
  #+begin_src c :tangle (eval h_private_func) :noexport
 void* qmckl_malloc(qmckl_context context,
                   const qmckl_memory_info_struct info);
  #+end_src

+  Here's a step-by-step explanation of ~qmckl_malloc~:
+
+  1. The function takes two parameters: a ~qmckl_context~ and a
+ ~qmckl_memory_info_struct~ containing the desired size of the memory
+     block to allocate.
+
+  2. The function checks if the provided ~qmckl_context~ is valid, using the
+ ~qmckl_context_check~ function.
+
+  3. The ~qmckl_context_struct~ pointer is retrieved from the provided
+ ~qmckl_context~.
+
+  4. The function then allocates memory:
+     If the ~HAVE_HPC~ and ~HAVE_POSIX_MEMALIGN~ macros are defined, the memory
+     allocation is done using the ~aligned_alloc~ function with a 64-byte alignment,
+     rounding up the requested size to the nearest multiple of 64 bytes. Else, the
+     memory allocation is done using the standard ~malloc~ function.
+
+  5  If the allocation fails, the function returns ~NULL~.
+
+  6. The allocated memory block is zeroed using ~memset~.
+
+  7. The function acquires a lock on the ~qmckl_context~ using ~qmckl_lock~.
+
+  8. Inside the locked section, the function checks if the
+ ~qmckl_memory_struct~ is full. If it is, it reallocates a larger array
+     by doubling its size and updating the ~array_size~ member of the
+ ~qmckl_memory_struct~.
+
+  9. The function finds the first available ~qmckl_memory_info_struct~ slot
+     in the element array of the ~qmckl_memory_struct~.
+
  # Source
  #+begin_src c :tangle (eval c)
 void* qmckl_malloc(qmckl_context context, const qmckl_memory_info_struct info) {
@ -119,10 +153,11 @@ void* qmckl_malloc(qmckl_context context, const qmckl_memory_info_struct info) {
  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;

 /* Allocate memory and zero it */
-#ifdef HAVE_HPC
-  void * pointer = aligned_alloc(64, ((info.size+64) >> 6) << 6 );
+  void * pointer = NULL;
+#if defined(HAVE_HPC) && defined(HAVE_POSIX_MEMALIGN)
+  if (posix_memalign(&pointer, 64, info.size) != 0) pointer = NULL;
 #else
-  void * pointer = malloc(info.size);
+  pointer = malloc(info.size);
 #endif
  if (pointer == NULL) {
    return NULL;
@ -159,6 +194,7 @@ void* qmckl_malloc(qmckl_context context, const qmckl_memory_info_struct info) {
    memcpy(&(ctx->memory.element[pos]), &info, sizeof(qmckl_memory_info_struct));
    ctx->memory.element[pos].pointer = pointer;
    ctx->memory.n_allocated += (size_t) 1;
+//printf("MALLOC: %5ld  %p\n", ctx->memory.n_allocated, ctx->memory.element[pos].pointer);
  }
  qmckl_unlock(context);

@ -203,6 +239,47 @@ qmckl_exit_code qmckl_free(qmckl_context context,
                           void * const ptr);
  #+end_src

+  Here's a step-by-step explanation of the ~qmckl_free~ function:
+
+  1. The function takes two parameters: a ~qmckl_context~ and a pointer to
+     the memory block (~ptr~) that needs to be deallocated.
+
+  2. The function checks if the provided ~qmckl_context~ is valid, using the
+ ~qmckl_context_check~ function. If it is not valid, it returns an error
+     code ~QMCKL_INVALID_CONTEXT~ using the ~qmckl_failwith~ function.
+   
+  3. The function checks if the provided pointer is ~NULL~. If it is, it
+     returns an error code ~QMCKL_INVALID_ARG_2~ using the ~qmckl_failwith~
+     function.
+   
+  4. The ~qmckl_context_struct~ pointer is retrieved from the provided
+ ~qmckl_context~.
+   
+  5. The function acquires a lock on the ~qmckl_context~ using ~qmckl_lock~.
+     
+  6. Inside the locked section, the function searches for the pointer in
+     the element array of the ~qmckl_memory_struct~.
+   
+  7. If the pointer is not found in the array, it releases the lock and
+     returns an error code ~QMCKL_INVALID_ARG_2~ using the ~qmckl_failwith~
+     function.
+     
+  8. If the pointer is found, the memory block is deallocated using the
+     standard ~free~ function.
+
+  9. The ~qmckl_memory_info_struct~ at the found position is zeroed
+     using ~memset~. This marks the slot as available for future
+     allocations.
+     
+ 10. The ~n_allocated~ member of the ~qmckl_memory_struct~ is decremented
+     by one, as the memory block has been deallocated.
+
+ 11. The function releases the lock on the ~qmckl_context~ using ~qmckl_unlock~.
+
+ 12. Finally, the function returns ~QMCKL_SUCCESS~ to indicate
+     successful deallocation of the memory block.
+     
+     
   # Source
   #+begin_src c :tangle (eval c)
 qmckl_exit_code qmckl_free(qmckl_context context, void * const ptr) {
@ -240,10 +317,13 @@ qmckl_exit_code qmckl_free(qmckl_context context, void * const ptr) {
                            "Pointer not found in context");
    }

+ /* Found */
+
    free(ptr);

-    memset( &(ctx->memory.element[pos]), 0, sizeof(qmckl_memory_info_struct) );
    ctx->memory.n_allocated -= (size_t) 1;
+//printf("FREE  : %5ld  %p\n", ctx->memory.n_allocated, ctx->memory.element[pos].pointer);
+    ctx->memory.element[pos] = qmckl_memory_info_struct_zero;
  }
  qmckl_unlock(context);

--- a/org/qmckl_mo.org
+++ b/org/qmckl_mo.org
--- a/org/qmckl_numprec.org
+++ b/org/qmckl_numprec.org
@ -1,4 +1,4 @@
-#+TITLE: Numerical precision
+3+TITLE: Numerical precision
 #+SETUPFILE: ../tools/theme.setup
 #+INCLUDE: ../tools/lib.org

@ -40,6 +40,42 @@ int main() {
 #include <stdlib.h>
 #include <string.h>

+#ifdef HAVE_FPE
+#define _GNU_SOURCE
+#include <fenv.h>
+#include <signal.h>
+#include <stdio.h>
+#include <execinfo.h>
+
+
+#define MAX_BACKTRACE_SIZE 100
+
+void floatingPointExceptionHandler(int signal) {
+    void* backtraceArray[MAX_BACKTRACE_SIZE];
+    int backtraceSize = backtrace(backtraceArray, MAX_BACKTRACE_SIZE);
+    char** backtraceSymbols = backtrace_symbols(backtraceArray, backtraceSize);
+
+    // Print the backtrace
+    for (int i = 0; i < backtraceSize; ++i) {
+            printf("[%d] %s\n", i, backtraceSymbols[i]);
+    }
+
+    // Clean up the memory used by backtrace_symbols
+    free(backtraceSymbols);
+
+    exit(EXIT_FAILURE);
+}
+
+static void __attribute__ ((constructor))
+trapfpe ()
+{
+ /* Enable some exceptions.  At startup all exceptions are masked.  */
+
+  feenableexcept (FE_INVALID|FE_DIVBYZERO|FE_OVERFLOW);
+  signal(SIGFPE, floatingPointExceptionHandler);
+}
+#endif
+
 #include "qmckl.h"
 #include "qmckl_context_private_type.h"
  #+end_src
@ -50,8 +86,8 @@ int main() {
 default parameters determining the target numerical precision and
 range are defined. Following the IEEE Standard for Floating-Point
 Arithmetic (IEEE 754),
- /precision/ refers to the number of significand bits and /range/
- refers to the number of exponent bits.
+ /precision/ refers to the number of significand bits (including the
+ sign bit) and /range/ refers to the number of exponent bits.

 #+NAME: table-precision
 | ~QMCKL_DEFAULT_PRECISION~ | 53 |
@ -293,23 +329,28 @@ int qmckl_get_numprec_range(const qmckl_context context) {

 * Helper functions

- ~qmckl_get_numprec_epsilon~ returns $\epsilon = 2^{1-n}$ where ~n~ is the precision.
-  We need to remove the sign bit from the precision.
+** Epsilon

-  #+begin_src c :comments org :tangle (eval h_func) :exports none
+ ~qmckl_get_numprec_epsilon~ returns $\epsilon = 2^{1-n}$ where ~n~ is the precision.
+   We need to remove the sign bit from the precision.
+
+   #+begin_src c :comments org :tangle (eval h_func) :exports none
 double qmckl_get_numprec_epsilon(const qmckl_context context);
-  #+end_src
+   #+end_src

-  # Source
-  #+begin_src c :tangle (eval c)
+   # Source
+   #+begin_src c :tangle (eval c)
 double qmckl_get_numprec_epsilon(const qmckl_context context) {
-  const int precision = qmckl_get_numprec_precision(context);
-  return 1. /  (double) (1L << (precision-2));
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT)
+    return QMCKL_INVALID_CONTEXT;
+  const qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
+  const int precision = ctx->numprec.precision;
+  return 1. /  (double) ( ((uint64_t) 1) << (precision-2));
 }
-  #+end_src
+   #+end_src

-  # Fortran interface
-  #+begin_src f90 :tangle (eval fh_func) :exports none
+   # Fortran interface
+   #+begin_src f90 :tangle (eval fh_func) :exports none
  interface
     real (c_double) function qmckl_get_numprec_epsilon(context) bind(C)
       use, intrinsic :: iso_c_binding
@ -317,7 +358,137 @@ double qmckl_get_numprec_epsilon(const qmckl_context context) {
       integer (qmckl_context), intent(in), value :: context
     end function qmckl_get_numprec_epsilon
  end interface
-  #+end_src
+   #+end_src
+
+** Testing the number of unchanged bits
+
+   To test that a given approximation keeps a given number of bits
+   unchanged, we need a function that returns the number of unchanged
+   bits in the range, and in the precision.
+
+   For this, we first count by how many units in the last place (ulps) two
+   numbers differ.
+
+   #+begin_src c :tangle (eval c)
+int64_t countUlpDifference_64(double a, double b) {
+
+  union int_or_float {
+    int64_t i;
+    double f;
+  } x, y;
+
+  x.f = a;
+  y.f = b;
+
+  // Handle sign bit discontinuity: if the signs are different and either value is not zero
+  if ((x.i < 0) != (y.i < 0) && (x.f != 0.0) && (y.f != 0.0)) {
+    // Use the absolute values and add the distance to zero for both numbers
+    int64_t distanceToZeroForX = x.i < 0 ? INT64_MAX + x.i : INT64_MAX - x.i;
+    int64_t distanceToZeroForY = y.i < 0 ? INT64_MAX + y.i : INT64_MAX - y.i;
+    return distanceToZeroForX + distanceToZeroForY;
+  }
+
+  // Calculate the difference in their binary representations
+  int64_t result = x.i - y.i;
+  result = result > 0 ? result : -result;
+  return result;
+}
+   #+end_src
+
+   #+begin_src c :comments org :tangle (eval h_func) :exports none
+int32_t qmckl_test_precision_64(double a, double b);
+int32_t qmckl_test_precision_32(float a, float b);
+   #+end_src
+
+   #+begin_src c :tangle (eval c)
+int32_t qmckl_test_precision_64(double a, double b) {
+
+  int64_t diff = countUlpDifference_64(a,b);
+
+  if (diff == 0) return 53;
+
+  int32_t result = 53;
+
+  for (int i=0 ; i<53 && diff != 0 ; ++i) {
+    diff >>= 1;
+    result--;
+  }
+
+  return result;
+}
+   #+end_src
+
+
+   #+begin_src c :tangle (eval c)
+int32_t qmckl_test_precision_32(float a, float b) {
+  return qmckl_test_precision_64( (double) a, (double) b );
+}
+   #+end_src
+
+   #+begin_src f90 :tangle (eval fh_func) :exports none
+  interface
+     integer (c_int) function qmckl_test_precision_32(a,b) bind(C)
+       use, intrinsic :: iso_c_binding
+       import
+       real (c_float), intent(in), value :: a, b
+     end function qmckl_test_precision_32
+  end interface
+
+  interface
+     integer (c_int) function qmckl_test_precision_64(a,b) bind(C)
+       use, intrinsic :: iso_c_binding
+       import
+       real (c_double), intent(in), value :: a, b
+     end function qmckl_test_precision_64
+  end interface
+   #+end_src
+
+* Approximate functions
+
+** Exponential
+
+   Fast exponential function, adapted from Johan Rade's implementation
+   (https://gist.github.com/jrade/293a73f89dfef51da6522428c857802d). It
+   is based on Schraudolph's paper:
+
+        N. Schraudolph, "A Fast, Compact Approximation of the Exponential Function",
+ /Neural Computation/ *11*, 853–862 (1999).
+        (available at https://nic.schraudolph.org/pubs/Schraudolph99.pdf)
+
+   #+begin_src c :tangle (eval c)
+float fastExpf(float x)
+{
+  const float a = 12102203.0;
+  const float b = 1064986816.0;
+  x = a * x + b;
+
+  const float c = 8388608.0;
+  const float d = 2139095040.0;
+  if (x < c || x > d)
+    x = (x < c) ? 0.0f : d;
+
+  uint32_t n = (uint32_t) x;
+  memcpy(&x, &n, 4);
+  return x;
+}
+
+
+double fastExp(double x)
+{
+  const double a = 6497320848556798.0;
+  const double b = 4606985713057410560.0;
+  x = a * x + b;
+
+  const double c = 4503599627370496.0;
+  const double d = 9218868437227405312.0;
+  if (x < c || x > d)
+    x = (x < c) ? 0.0 : d;
+
+  uint64_t n = (uint64_t) x;
+  memcpy(&x, &n, 8);
+  return x;
+}
+   #+end_src

 * End of files                                                     :noexport:

--- a/org/qmckl_point.org
+++ b/org/qmckl_point.org
@ -19,7 +19,7 @@ coordinates of all the walkers.
 #include <stdbool.h>
 #include "qmckl_blas_private_type.h"
  #+end_src
-  
+
  #+begin_src c :tangle (eval h_private_func)
 #ifndef QMCKL_POINT_HPF
 #define QMCKL_POINT_HPF
@ -318,8 +318,8 @@ qmckl_set_point (qmckl_context context,
  if (num != ctx->point.num) {

    if (ctx->point.coord.data != NULL) {
-      rc = qmckl_matrix_free(context, &(ctx->point.coord));
-      assert (rc == QMCKL_SUCCESS);
+        rc = qmckl_matrix_free(context, &(ctx->point.coord));
+        assert (rc == QMCKL_SUCCESS);
    }

    ctx->point.coord = qmckl_matrix_alloc(context, num, 3);
--- a/org/qmckl_sherman_morrison_woodbury.org
+++ b/org/qmckl_sherman_morrison_woodbury.org
--- a/org/qmckl_tests.org
+++ b/org/qmckl_tests.org
@ -60114,8 +60114,8 @@ double n2_elec_coord[n2_walk_num][n2_elec_num][3] = { {
 #define n2_dim_c_vec        ((int64_t) 23)

 int64_t n2_type_nucl_vector[n2_nucl_num] = {
-  1,
-  1};
+  0,
+  0};

 double n2_a_vector[n2_aord_num + 1][n2_type_nucl_num] = {
  { 0.      },
--- a/org/qmckl_trexio.org
+++ b/org/qmckl_trexio.org
@ -267,7 +267,7 @@ qmckl_trexio_read_nucleus_X(qmckl_context context, trexio_t* const file)


   #+begin_src c :tangle (eval c)
-
+  assert ( qmckl_nucleus_provided(context) );
  return QMCKL_SUCCESS;
 }
 #endif
@ -368,7 +368,7 @@ qmckl_trexio_read_ao_X(qmckl_context context, trexio_t* const file)
 *** Number of atomic orbitals

   #+begin_src c :tangle (eval c)
-  int64_t ao_num = 0LL;
+  int64_t ao_num = 0;

  rcio = trexio_read_ao_num_64(file, &ao_num);
  if (rcio != TREXIO_SUCCESS) {
@ -515,7 +515,7 @@ qmckl_trexio_read_ao_X(qmckl_context context, trexio_t* const file)
    for (int i=0 ; i<nucleus_num ; ++i) {
      nucleus_shell_num[i] = 0;
    }
-    
+
    for (int i=0 ; i<shell_num ; ++i) {
      const int k = tmp_array[i];
      if (k < 0 || k >= nucleus_num) {
@ -601,7 +601,7 @@ qmckl_trexio_read_ao_X(qmckl_context context, trexio_t* const file)
    if (shell_prim_num == NULL) {
      return qmckl_failwith( context,
                             QMCKL_ALLOCATION_FAILED,
-                             "qmckl_trexio_read_basis_shell_prim_num_X",
+                             "qmckl_trexio_read_basis_shell_index",
                             NULL);
    }

@ -617,7 +617,7 @@ qmckl_trexio_read_ao_X(qmckl_context context, trexio_t* const file)
      shell_prim_num = NULL;
      return qmckl_failwith( context,
                             QMCKL_ALLOCATION_FAILED,
-                             "qmckl_trexio_read_basis_shell_prim_num_X",
+                             "qmckl_trexio_read_basis_shell_index",
                             NULL);
    }

@ -632,7 +632,7 @@ qmckl_trexio_read_ao_X(qmckl_context context, trexio_t* const file)
      tmp_array = NULL;
      return qmckl_failwith( context,
                             QMCKL_FAILURE,
-                             "trexio_read_basis_shell_prim_num",
+                             "qmckl_trexio_read_basis_shell_index",
                             trexio_string_of_error(rcio));
    }

@ -640,16 +640,18 @@ qmckl_trexio_read_ao_X(qmckl_context context, trexio_t* const file)
    for (int i=0 ; i<shell_num ; ++i) {
      shell_prim_num[i] = 0;
    }
-    
+
    for (int i=0 ; i<prim_num ; ++i) {
      const int k = tmp_array[i];
      if (k < 0 || k >= shell_num) {
        qmckl_free(context, tmp_array);
        qmckl_free(context, shell_prim_num);
+        char msg[128];
+        sprintf(&msg[0], "Irrelevant data in TREXIO file: k = %d", k);
        return qmckl_failwith( context,
                                QMCKL_FAILURE,
-                                "trexio_read_basis_shell_prim_num",
-                                "Irrelevant data in TREXIO file");
+                                "qmckl_trexio_read_basis_shell_index",
+                                &msg[0]);
      }
      shell_prim_num[k] += 1;
    }
@ -1030,7 +1032,7 @@ qmckl_trexio_read_mo_X(qmckl_context context, trexio_t* const file)
                             trexio_string_of_error(rcio));
    }

-    rc = qmckl_set_mo_basis_coefficient(context, mo_coef);
+    rc = qmckl_set_mo_basis_coefficient(context, mo_coef, ao_num*mo_num);

    qmckl_free(context, mo_coef);
    mo_coef = NULL;
@ -1332,7 +1334,7 @@ double * mo_coef = (double*) malloc (ao_num * mo_num * sizeof(double));
 rc = qmckl_get_mo_basis_coefficient(context, mo_coef, mo_num*ao_num);
 qmckl_check(context, rc);
 for (int i=0 ; i<ao_num * mo_num ; i++) {
-  printf("%d %e %e %e\n", i, mo_coef[i], chbrclf_mo_coef[i], 
+  printf("%d %e %e %e\n", i, mo_coef[i], chbrclf_mo_coef[i],
         ( fabs(mo_coef[i] - chbrclf_mo_coef[i])/fabs(mo_coef[i])) );
  assert ( fabs(mo_coef[i] - chbrclf_mo_coef[i])/fabs(mo_coef[i]) < 1.e-12 );
 }
--- a/org/table_of_contents
+++ b/org/table_of_contents
@ -12,7 +12,7 @@ qmckl_ao.org
 qmckl_mo.org
 qmckl_determinant.org
 qmckl_sherman_morrison_woodbury.org
-qmckl_jastrow.org
+qmckl_jastrow_champ.org
 qmckl_local_energy.org
 qmckl_utils.org
 qmckl_trexio.org
--- a/python/setup.py
+++ b/python/setup.py
@ -46,7 +46,7 @@ qmckl_module   =  Extension(name                 = "._" + MODULE_NAME,


 setup(name             = MODULE_NAME,
-      version          = "0.3.1",
+      version          = "1.0.0",
      author           = "TREX-CoE",
      author_email     = "posenitskiy@irsamc.ups-tlse.fr",
      description      = """Python API of the QMCkl library""",
--- a/python/src/qmckl.i
+++ b/python/src/qmckl.i
@ -58,6 +58,11 @@ import_array();
 %apply ( double* ARGOUT_ARRAY1 , int64_t DIM1 ) { ( double* const C, const int64_t size_max_C) };
 %apply ( double* ARGOUT_ARRAY1 , int64_t DIM1 ) { ( double* const B, const int64_t size_max_B) };

+%apply ( int64_t* IN_ARRAY1 , int64_t DIM1 ) { ( const int64_t* A, const int64_t size_max_A) };
+%apply ( int64_t* IN_ARRAY1 , int64_t DIM1 ) { ( const int64_t* B, const int64_t size_max_B) };
+%apply ( int64_t* ARGOUT_ARRAY1 , int64_t DIM1 ) { ( int64_t* const C, const int64_t size_max_C) };
+%apply ( int64_t* ARGOUT_ARRAY1 , int64_t DIM1 ) { ( int64_t* const B, const int64_t size_max_B) };
+
 /* Handle properly get_point */


--- a/share/qmckl/test_data/sm_test.h
+++ b/share/qmckl/test_data/sm_test.h
--- a/tools/build_makefile.py
+++ b/tools/build_makefile.py
@ -73,7 +73,7 @@ def main():
        TEXT[org]     = text
        HTML[org]     = html

-        grep = open(org, "r").read()
+        grep = open(org, "r", encoding="utf-8").read()

        if "(eval c)" in grep:
           C_FILES   += [c]
@ -135,14 +135,14 @@ def main():
            F_FILES += [F90]

            if F90 in DEPS:
-                DEPS[F90] += [tangled, "$(src_qmckl_fo)"]
+                DEPS[F90] += [tangled, "$(qmckl_fo)"]
            else:
-                DEPS[F90]  = [tangled, "$(src_qmckl_fo)"]
+                DEPS[F90]  = [tangled, "$(qmckl_fo)"]

            if fo in DEPS:
-                DEPS[fo] += [F90, "$(src_qmckl_fo)"]
+                DEPS[fo] += [F90, "$(qmckl_fo)"]
            else:
-                DEPS[fo]  = [F90, "$(src_qmckl_fo)"]
+                DEPS[fo]  = [F90, "$(qmckl_fo)"]

        if "(eval fh_func)" in grep:
            FH_FUNC_FILES += [fh_func]
@ -178,14 +178,14 @@ def main():
            F_TEST_FILES += [f_test]

            if f_test in DEPS:
-                DEPS_TEST[f_test] += [tangled, "$(test_qmckl_fo)"]
+                DEPS_TEST[f_test] += [tangled]
            else:
-                DEPS_TEST[f_test]  = [tangled, "$(test_qmckl_fo)"]
+                DEPS_TEST[f_test]  = [tangled]

            if c_test_x in TESTS:
-                TESTS[c_test_x] += [f_test, "$(test_qmckl_fo)"]
+                TESTS[c_test_x] += [f_test]
            else:
-                TESTS[c_test_x]  = [f_test, "$(test_qmckl_fo)"]
+                TESTS[c_test_x]  = [f_test]

    output = ["",
              "## Source files",
@ -232,7 +232,7 @@ def main():
        if DEPS[f][-1].endswith(".tangled"):
            output += [ f + ": " + " ".join(DEPS[f]) ]
    output += [ "endif",
-                "$(src_qmckl_fo): $(src_qmckl_f)" ]
+                "$(qmckl_fo): $(qmckl_f)" ]
    for f in sorted(DEPS.keys()):
        if not DEPS[f][-1].endswith(".tangled"):
            output += [ f + ": include/config.h " + " ".join(DEPS[f]) ]
@ -240,8 +240,7 @@ def main():
    output+= ["",
              "## Test files",
              "",
-              "$(header_tests): $(TANGLED_FILES)",
-              "$(test_qmckl_fo): $(test_qmckl_f)"]
+              "$(header_tests): $(TANGLED_FILES)" ]
    output += ["",
                "check_PROGRAMS = $(TESTS)" ]
    for f in sorted(TESTS.keys()):
--- a/tools/build_qmckl_f.sh
+++ b/tools/build_qmckl_f.sh
@ -76,7 +76,7 @@ cat << EOF > ${qmckl_f}
 !    
 !    
 !
-module qmckl
+module qmckl_constants
  use, intrinsic :: iso_c_binding
 EOF

@ -85,6 +85,13 @@ do
    cat $i >> ${qmckl_f}
 done

+cat << EOF >> ${qmckl_f}
+end module qmckl_constants
+
+module qmckl
+  use qmckl_constants
+EOF
+
 for i in ${HEADERS}
 do
    cat $i >> ${qmckl_f}
--- a/tools/init.el
+++ b/tools/init.el
@ -58,8 +58,8 @@
 (add-hook 'org-babel-after-execute-hook 'org-display-inline-images)
 '(indent-tabs-mode nil)

-(require 'evil)
 (setq evil-want-C-i-jump nil)
+(require 'evil)
 (evil-mode 1)
 (global-font-lock-mode t)
 (global-superword-mode 1)
--- a/tools/lib.org
+++ b/tools/lib.org
@ -32,14 +32,13 @@
     | ~C~       | ~double[n][ldc]~ | out    | Array containing the $m \times n$ matrix $C$  |
     | ~ldc~     | ~int64_t~        | in     | Leading dimension of array ~C~                |

-
 *** Fortran-C type conversions

    #+NAME:f_of_c
-    #+BEGIN_SRC python :var table=test :var rettyp="integer" :var fname=[] :results value :noweb yes :wrap "src f90 :tangle (eval f) :comments org :exports none"
+    #+BEGIN_SRC python :var table=test :var rettyp="qmckl_exit_code" :var fname=[] :results value :noweb yes :wrap "src f90 :tangle (eval f) :comments org :exports none"
 f_of_c_d = { '' : ''
-     , 'qmckl_context'   : 'integer (c_int64_t)'
-     , 'qmckl_exit_code' : 'integer (c_int32_t)'
+     , 'qmckl_context'   : 'integer (qmckl_context)'
+     , 'qmckl_exit_code' : 'integer (qmckl_exit_code)'
     , 'bool'            : 'logical*8'
     , 'int32_t'         : 'integer (c_int32_t)'
     , 'int64_t'         : 'integer (c_int64_t)'
@ -47,15 +46,15 @@ f_of_c_d = { '' : ''
     , 'uint64_t'        : 'integer (c_int64_t)'
     , 'float'           : 'real    (c_float  )'
     , 'double'          : 'real    (c_double )'
-     , 'char'            : 'character'
+     , 'char'            : 'character(c_char  )'
 }
    #+END_SRC

    #+NAME:c_of_f
    #+BEGIN_SRC python :var table=test :var rettyp="integer" :var fname=[] :results value :noweb yes :wrap "src f90 :tangle (eval f) :comments org :exports none"
 ctypeid_d = { '' : ''
-     , 'qmckl_context'   : 'integer(c_int64_t)'
-     , 'qmckl_exit_code' : 'integer(c_int32_t)'
+     , 'qmckl_context'   : 'integer(qmckl_context)'
+     , 'qmckl_exit_code' : 'integer(qmckl_exit_code)'
     , 'integer'         : 'integer(c_int32_t)'
     , 'integer*8'       : 'integer(c_int64_t)'
     , 'real'            : 'real(c_float)'
@ -132,26 +131,40 @@ return template

    #+END_SRC

-    #+RESULTS: generate_c_header
-    #+begin_src c :tangle (eval h_func) :comments org
-    qmckl_exit_code [] (
-          const qmckl_context context,
-          const char transa,
-          const char transb,
-          const int64_t m,
-          const int64_t n,
-          const double* A,
-          const int64_t lda,
-          const double* B,
-          const int64_t ldb,
-          double* const C,
-          const int64_t ldc );
-    #+end_src
+    #+NAME: generate_private_c_header
+    #+BEGIN_SRC python :var table=test :var rettyp="qmckl_exit_code" :var fname=[] :results drawer :noweb yes :wrap "src c :tangle (eval h_private_func) :comments org"
+<<parse_table>>

+results = []
+for d in parse_table(table):
+    name = d["name"]
+    c_type = d["c_type"]
+
+    # Add star for arrays
+    if d["rank"] > 0 or d["inout"] in ["out", "inout"]:
+        c_type += "*"
+
+    if d["inout"] == "out":
+        c_type += " const"
+
+    # Only inputs are const
+    if d["inout"] == "in":
+        const = "const "
+    else:
+        const = ""
+
+    results += [ f"      {const}{c_type} {name}" ]
+
+results=',\n'.join(results)
+template = f"""{rettyp} {fname} (\n{results} ); """
+return template
+
+    #+END_SRC
+    
 *** Generates a C interface to the Fortran function

    #+NAME: generate_c_interface
-    #+BEGIN_SRC python :var table=[]   :var rettyp="integer" :var fname=[] :results value :noweb yes :wrap "src f90 :tangle (eval f) :comments org :exports none"
+    #+BEGIN_SRC python :var table=[]   :var rettyp="qmckl_exit_code" :var fname=[] :results value :noweb yes :wrap "src f90 :tangle (eval f) :comments org :exports none"
 <<c_of_f>>
 <<f_of_c>>
 <<parse_table>>
@ -168,7 +181,7 @@ results = [ f"{rettyp_c} function {fname} &"
 , f"    ({args}) &"
 ,  "    bind(C) result(info)"
 ,  ""
-,  "  use, intrinsic :: iso_c_binding"
+,  "  use qmckl_constants"
 ,  "  implicit none"
 ,  ""
 ]
@ -207,7 +220,7 @@ return results
 *** Generates a Fortran interface to the C function

    #+NAME: generate_f_interface
-    #+BEGIN_SRC python :var table=test :var rettyp="integer" :var fname=[] :results value :noweb yes :wrap "src f90 :tangle (eval fh_func) :comments org :exports none"
+    #+BEGIN_SRC python :var table=test :var rettyp="qmckl_exit_code" :var fname=[] :results value :noweb yes :wrap "src f90 :tangle (eval fh_func) :comments org :exports none"
 <<c_of_f>>
 <<f_of_c>>
 <<parse_table>>
@ -221,7 +234,7 @@ results = [ f"interface"
 , f"  {rettyp_c} function {fname} &"
 , f"      ({args}) &"
 ,  "      bind(C)"
-,  "    use, intrinsic :: iso_c_binding"
+,  "    use qmckl_constants"
 ,  "    import"
 ,  "    implicit none"
 ,  ""
@ -255,7 +268,54 @@ results='\n'.join(results)
 return results
    #+END_SRC

+    #+NAME: generate_private_f_interface
+    #+BEGIN_SRC python :var table=test :var rettyp="qmckl_exit_code" :var fname=[] :results value :noweb yes :wrap "src f90 :tangle (eval fh_private_func) :comments org :exports none"
+<<c_of_f>>
+<<f_of_c>>
+<<parse_table>>
+d = parse_table(table)

+args = ", ".join([ x["name"] for x in d ])
+
+rettyp_c = ctypeid_d[rettyp.lower()]
+
+results = [ f"interface"
+, f"  {rettyp_c} function {fname} &"
+, f"      ({args}) &"
+,  "      bind(C)"
+,  "    use qmckl_constants"
+,  "    import"
+,  "    implicit none"
+,  ""
+]
+
+for d in parse_table(table):
+    f_type = f_of_c_d[d["c_type"]]
+    inout  = "intent("+d["inout"]+")"
+    name = d["name"]
+
+    # Input scalars are passed by value
+    if d["rank"] == 0 and d["inout"] == "in":
+        value = ", value"
+    else:
+        value = "       "
+
+    # Append dimensions to the name
+    if d["rank"] == 0:
+        dims = ""
+    else:
+        d["dims"].reverse()
+        dims = "(" + ",".join(d["dims"]) + ")"
+
+    results += [ f"    {f_type:20}, {inout:12}{value} :: {name}{dims}" ]
+
+results += [ ""
+, f"  end function {fname}"
+, f"end interface"
+]
+results='\n'.join(results)
+return results
+    #+END_SRC

 ** Creating provide functions

@ -421,3 +481,4 @@ return msg
      return QMCKL_SUCCESS;
    }
    #+end_src
+    
--- a/tools/tangle.sh
+++ b/tools/tangle.sh
@ -18,8 +18,13 @@ if [[ -z ${top_builddir} ]] ; then
   exit 1
 fi

+EMACS="${VARIABLE:=emacs}" 
+
+EXTENSIONS="_f.F90 _fh_func.F90 _fh_type.F90 .c _func.h _type.h _private_type.h _private_func.h"
+
 function tangle()
 {
+    local backup_dir=$(mktemp -d)
    local org_file=$1
    local c_file=${org_file%.org}.c
    local f_file=${org_file%.org}.F90
@ -29,10 +34,32 @@ function tangle()
    elif [[ ${org_file} -ot ${f_file} ]] ; then
        return
    fi
+
+    local prefix=${top_builddir}/src/$(basename ${org_file})
+    prefix=${prefix%.org}
+    for ext in $EXTENSIONS ; do
+      if [[ -f ${prefix}${ext} ]] ; then
+         mv ${prefix}${ext} ${backup_dir}
+      fi
+    done
+
    ${srcdir}/tools/missing \
-        emacs --no-init-file --no-site-lisp --quick --batch ${org_file} \
+        $EMACS --no-init-file --no-site-lisp --quick --batch ${org_file} \
         --load=${srcdir}/tools/config_tangle.el \
        -f org-babel-tangle
+    
+    for ext in $EXTENSIONS ; do
+      local new_file=${prefix}${ext}
+      local old_file=${backup_dir}/$(basename ${new_file})
+      diff $new_file $old_file &> /dev/null
+      if [[ $? -eq 0 ]] ; then
+         echo "${old_file} unchanged" 
+         mv $old_file $new_file
+      else
+         echo "${old_file}   changed" 
+      fi
+    done
+    rm -rf ${backup_dir}
 }

 for i in $@
Author	SHA1	Message	Date
Anthony Scemama	9553059bbe	Fix wrong gradient at nodes of AOs	2024-05-08 13:56:32 +02:00
Anthony Scemama	8c5ec872ed	CHange version	2024-04-10 16:38:10 +02:00
Anthony Scemama	5ee297a1c9	Removed scorecards	2024-03-29 09:27:40 +01:00
Anthony Scemama	5040ce1087	Removed redundant iso_c_binding	2024-02-26 09:23:23 +01:00
Anthony Scemama	fd9ce7ed5e	Deactivated hpc versions woodburry	2024-02-24 01:18:06 +02:00
Anthony Scemama	574cde88e5	Removed redundant iso_c_binding	2024-02-24 00:58:11 +02:00
Anthony Scemama	e0abd84059	Fix Cray fortran errors	2024-02-23 16:37:31 +01:00
Anthony Scemama	b2395ece87	Ordering problem in Fortran interface fixed	2024-02-23 16:15:56 +01:00
Anthony Scemama	f745899f3a	Replaced iso_c_binding by qmckl_constants	2024-02-23 12:06:13 +01:00
Anthony Scemama	83dea2b773	Fix warnings with Cray ftn	2024-02-23 11:56:28 +01:00
Anthony Scemama	21f40b3a13	Fixed gradient of Jastrow in HPC version	2024-02-21 16:33:49 +01:00
Anthony Scemama	e07b8bfa55	Avoid memset in Jastrow	2024-02-20 23:59:28 +01:00
Anthony Scemama	41615ba14b	Avoid memset in Jastrow	2024-02-20 23:38:34 +01:00
Anthony Scemama	2f0ca9f674	Improved Jastrow	2024-02-14 11:26:10 +01:00
Anthony Scemama	be2a7199c2	simd	2024-02-14 11:11:50 +01:00
Anthony Scemama	2228ab23c5	Vectorization	2024-02-14 10:59:31 +01:00
Anthony Scemama	48b80f68f1	HPC implementation in Jastrow	2024-02-14 09:31:06 +01:00
Anthony Scemama	24e3f8dd11	OpenMP in Fortran function	2024-02-13 17:21:53 +01:00
Anthony Scemama	949cfb6f82	Accelerated Jastrow (OpenMP)	2024-02-13 17:07:59 +01:00
Anthony Scemama	6ce1d2cbb9	Merge branch 'master' of github.com:TREX-CoE/qmckl	2024-02-07 17:06:32 +01:00
Anthony Scemama	6caf3521a4	commented loops	2024-02-07 17:06:27 +01:00
Anthony Scemama	c6ea9c13a4	Removed ivdep	2024-02-07 17:04:18 +01:00
Anthony Scemama	e20be734cc	Merge branch 'master' of github.com:TREX-CoE/qmckl	2024-02-07 14:58:05 +01:00
Anthony Scemama	023c9cda85	Moved ivdep	2024-02-07 14:57:24 +01:00
Anthony Scemama	f08ed5da6d	Fixing memset	2024-02-06 22:39:52 +01:00
Anthony Scemama	3a73e5722b	Merge branch 'master' of github.com:TREX-CoE/qmckl	2024-02-06 22:27:27 +01:00
Anthony Scemama	a0e1843963	Fix memset	2024-02-06 22:27:17 +01:00
Anthony Scemama	5060bde30f	Moved ivdep after omp simd	2024-01-30 23:46:06 +01:00
Anthony Scemama	dd3db966b0	Replace += by = ... + for better FMA	2024-01-30 11:31:07 +01:00
Anthony Scemama	ffbeb97df4	Improve BLAS detection for ARM	2024-01-30 11:24:17 +01:00
Anthony Scemama	43ebd409a8	Improved vectorization of mo_value	2024-01-29 11:59:39 +01:00
Anthony Scemama	098b6deec3	Merge branch 'master' of github.com:TREX-CoE/qmckl	2024-01-11 14:33:48 +01:00
Anthony Scemama	d257e28b92	Fix bug in HPC laplacian AO	2024-01-11 14:33:40 +01:00
Anthony Scemama	7e1dad0e4e	Update README.md	2024-01-09 08:19:36 +01:00
Anthony Scemama	7fc10a47a1	Fixed memory leak	2023-11-30 19:14:32 +01:00
Anthony Scemama	43b4aa81bd	Revert to old way of computing ranges	2023-11-30 12:56:06 +01:00
Anthony Scemama	b1891b267e	New way to compute the nucleus range	2023-11-30 12:50:06 +01:00
Anthony Scemama	141a0a866e	Fixed tests	2023-11-30 01:22:08 +01:00
Anthony Scemama	dba15f6b84	Merge branch 'master' of github.com:TREX-CoE/qmckl	2023-11-30 01:17:23 +01:00
Anthony Scemama	27b1134a4c	Precision adjusted for MOs	2023-11-30 01:17:18 +01:00
Anthony Scemama	034f2e81e8	Merge pull request #112 from TREX-CoE/addman2-patch-1 Update qmckl_ao.org	2023-11-29 21:46:04 +01:00
addman2	2f69d2af21	Update qmckl_ao.org typo in documentation	2023-11-29 19:53:03 +01:00
Anthony Scemama	f150eb1610	Added functions to test the number of bits of precision	2023-11-29 11:41:16 +01:00
Anthony Scemama	4df54f21c7	Removed calloc	2023-11-29 02:10:20 +01:00
Anthony Scemama	c6d193887a	Fixed single precision	2023-11-29 01:18:15 +01:00
Anthony Scemama	dbb49a2df5	Introduced SP in ao->mo	2023-11-28 23:37:15 +01:00
Anthony Scemama	6bf9388a4e	Added control of precision in AOs	2023-11-28 17:00:39 +01:00
Anthony Scemama	063aada9e4	Added --with-icx	2023-11-28 12:44:35 +01:00
Anthony Scemama	952ca05bf0	Fix previous commit	2023-11-15 13:20:30 +01:00
Anthony Scemama	f1764a5717	Better error message in trexio read	2023-11-15 13:03:07 +01:00
Anthony Scemama	5d8dfacffe	Bump version to 0.5.4	2023-10-06 11:35:02 +02:00
Anthony Scemama	a7523fbf77	Now using posix_memalign	2023-10-06 11:33:33 +02:00
Anthony Scemama	eaa44b45c4	Fixed bug in Jastrow en HPC	2023-09-27 22:12:26 +02:00
Anthony Scemama	c70b7b246b	Fixed Jastrow GL for spin-independent Jastrow	2023-09-27 15:56:37 +02:00
Anthony Scemama	5118359099	Bump version	2023-09-26 17:37:17 +02:00
Anthony Scemama	0ddaf0cd29	cast	2023-09-26 17:34:57 +02:00
Anthony Scemama	89a4a57c32	Added spin-independent Jastrow	2023-09-26 17:32:51 +02:00
Anthony Scemama	ad378103a5	Fixed cast	2023-09-26 16:25:24 +02:00
Anthony Scemama	1dc1c0f192	Fixed warning with Clang	2023-09-26 16:10:37 +02:00
Anthony Scemama	233edeeae2	Fixed warning with Clang	2023-09-26 16:02:10 +02:00
Anthony Scemama	47c4ee7d01	Fixed parallel build with Fortran	2023-09-26 16:00:49 +02:00
Anthony Scemama	ab596fe408	Fix CI	2023-09-26 15:18:22 +02:00
Anthony Scemama	de98045fe4	Simplifying Fortran	2023-09-22 16:56:48 +02:00
Anthony Scemama	0d9af3c497	Cleaning Fortran	2023-09-22 16:41:43 +02:00
Anthony Scemama	50fa3aa754	Introduced qmckl_constants module	2023-09-22 09:33:54 +02:00
Anthony Scemama	7a995a0f6b	Simplify Fortran call	2023-09-21 13:02:13 +02:00
Anthony Scemama	0d2327cae3	Accelerated 1-body Jastrow	2023-09-21 12:37:57 +02:00
Anthony Scemama	46785ddb7e	Fixing IVDEP CI	2023-09-15 10:23:29 +02:00
Anthony Scemama	6939891ac3	Bug in J_gl	2023-09-14 17:41:45 +02:00
Anthony Scemama	42222f73a5	Merge olympe2:qmckl	2023-09-14 11:01:13 +02:00
Anthony Scemama	6919e16f9a	Update version to 0.5.2	2023-09-14 11:01:08 +02:00
Anthony Scemama	561373fe4f	Improved configure for nvc	2023-09-14 11:00:24 +02:00
Anthony Scemama	c66188e641	Cleaning in SM	2023-09-14 09:56:28 +02:00
Anthony Scemama	932263d22f	Cleaning in SM	2023-09-14 09:54:50 +02:00
Anthony Scemama	10ee050050	Removed IVPDEP from SM	2023-09-14 09:53:17 +02:00
Anthony Scemama	2291103a9b	Merge branch 'master' of github.com:TREX-CoE/qmckl Conflicts: org/qmckl_sherman_morrison_woodbury.org	2023-09-14 09:43:59 +02:00
Anthony Scemama	fd2addb370	Added IVDEP and ALIGNED in configure.ac	2023-09-14 09:41:15 +02:00
Anthony Scemama	d77dc26e52	Fix bugs for Python and Jastrow	2023-09-13 16:22:23 +02:00
Anthony Scemama	3db1765cdb	Reduced OMP stack size	2023-09-13 13:23:16 +02:00
Anthony Scemama	5303bf88b3	Fixed bug in Jastrow gradient and Laplacian	2023-09-12 17:50:13 +02:00
Anthony Scemama	4133a6ba0e	Fixed bug in 3-body Jastrow. Tested in CHAMP->correct	2023-09-12 16:41:26 +02:00
Anthony Scemama	fc7d8b2af5	Testing CI	2023-09-12 11:28:14 +02:00
Anthony Scemama	1b5ce50b78	Testing CI	2023-09-12 11:25:53 +02:00
Anthony Scemama	a5139dd66d	Update ax_openmp.m4	2023-09-12 11:22:27 +02:00
Anthony Scemama	85e7592a6a	Workflows	2023-09-12 11:12:55 +02:00
Anthony Scemama	cbc8b9bd58	Fixes in Fortran interface	2023-09-11 17:05:41 +02:00
Anthony Scemama	c176cd86c3	Update version in python	2023-08-31 18:13:37 +02:00
Anthony Scemama	ecc19af2ba	Remove GPU from configure	2023-08-31 17:31:27 +02:00
Anthony Scemama	a14d8abd52	Example in C	2023-08-31 12:05:37 +02:00
Anthony Scemama	5d1373a2fb	Version 0.4.1	2023-08-29 11:26:16 +02:00
Anthony Scemama	bbf596bb4c	Fix AO bug in HPC	2023-08-24 10:14:37 +02:00
Anthony Scemama	803f914fb3	Initialize ao_vgl to zero	2023-08-24 09:59:57 +02:00
Anthony Scemama	1f3a08fa30	Avoid FPE in Jastrow GL	2023-08-23 14:32:33 +02:00
Anthony Scemama	2a38543ba0	Avoid FPE in Jastrow	2023-08-23 13:58:51 +02:00
Anthony Scemama	71d271572e	Check values of type_nucl in Jastrow	2023-08-23 13:40:32 +02:00
Anthony Scemama	5b64b44b1f	Merge pull request #111 from fmgjcoppens/wb3x3 Wb3x3	2023-07-13 09:44:50 +02:00
Anthony Scemama	38a7c3ba6f	Merge pull request #109 from fmgjcoppens/wb2x2_doc Wb2x2 doc kernel	2023-07-13 09:44:31 +02:00
Francois Coppens	44323468e6	Added woodbury 3x3 HPC and DOC kernels back. They pass the tests.	2023-07-12 15:48:43 +02:00
Francois Coppens	81be8263a0	Fixed mismached array size in splitting doc kernel. Test passes now.	2023-07-12 15:48:05 +02:00
Francois Coppens	5f888abe5b	Added missing headers and interfaces.	2023-07-12 15:48:04 +02:00
Francois Coppens	17398059d5	Woodbury 2x2 doc kernel passed test.	2023-07-12 15:48:04 +02:00
Francois Coppens	bb8377edc8	Made documentation vectors, matrices, conversion functions square for easier code comprehenisibility.	2023-07-12 15:48:04 +02:00
Francois Coppens	41745409a0	Doc kernel.	2023-07-12 15:48:04 +02:00
Francois Coppens	883416d84c	Added WB2 kernel back.	2023-07-12 15:48:04 +02:00
Anthony Scemama	db13db8afa	Enable FPE tracking	2023-06-29 15:56:11 +02:00
Anthony Scemama	5228c287ad	Fix #110	2023-06-21 14:42:40 +02:00
Anthony Scemama	9211bf576f	Fix typo	2023-06-14 11:45:39 +02:00
Anthony Scemama	53df240ca3	Fix typo	2023-06-14 11:44:00 +02:00
Anthony Scemama	15e3c7a4c8	Trick to make MKL efficient on AMD	2023-06-14 09:11:58 +02:00
Anthony Scemama	cf14f1c32b	Merge branch 'master' of github.com:TREX-CoE/qmckl	2023-05-26 09:51:52 +02:00
Anthony Scemama	dee0054c34	More compact error checking in Jastrow	2023-05-26 09:51:15 +02:00
Anthony Scemama	42d89b6e2f	Delete devskim.yml	2023-05-25 03:00:33 +02:00
Anthony Scemama	e3f99d0030	Jastrow ee and en OK in QMC=Chem	2023-05-25 01:38:50 +02:00
Anthony Scemama	edbe33f40f	Updated tests for Jastrow with kappa=0.6	2023-05-25 01:12:05 +02:00
Anthony Scemama	1e4bf9631f	Tests for qmckl_distance_rescaled	2023-05-24 23:56:19 +02:00
Anthony Scemama	5c019b06e3	Fixed sign error in jastrow gradient	2023-05-24 12:41:45 +02:00
Anthony Scemama	04d599649b	Rewrote qmckl_compute_jastrow_champ_factor_ee_gl_hpc	2023-05-24 11:32:23 +02:00
Anthony Scemama	7987d6a18a	Renamed deriv_e in to gl	2023-05-24 11:12:15 +02:00
Anthony Scemama	19a0a4a675	Rewrote Jastrow ee	2023-05-23 13:49:26 +02:00
Anthony Scemama	cfda515885	Improved HPC jastrow fee	2023-05-23 09:51:55 +02:00
Anthony Scemama	92705b7c87	Debugging Jastrow	2023-05-22 19:15:17 +02:00
Anthony Scemama	95b579dfc8	More flexibility in setting Jastrow	2023-05-19 16:35:05 +02:00
Anthony Scemama	252baa4721	ao_num=30 in H2 example	2023-05-02 14:24:30 +02:00
Anthony Scemama	7bec8b7984	Improved HPC of jastrow deriv	2023-04-11 19:16:14 +02:00
Anthony Scemama	b0bec3bc6c	Improved HPC of jastrow deriv	2023-04-11 18:51:05 +02:00
Anthony Scemama	5093b2c35c	HPC version of qmckl_compute_jastrow_champ_factor_een_deriv_e	2023-04-11 16:20:56 +02:00
Anthony Scemama	2153cfccf6	Optimized Jasrow	2023-03-31 19:58:30 +02:00
Anthony Scemama	e4023b426e	Fixed address sanitizer	2023-03-31 19:20:25 +02:00
Anthony Scemama	78cf825219	Caching of tangled files for better behvior of make	2023-03-31 15:10:30 +02:00
Anthony Scemama	5ae8828684	Jastrow OK	2023-03-31 14:41:32 +02:00
Anthony Scemama	daddb57200	Total jastrow value	2023-03-31 13:37:35 +02:00
Anthony Scemama	0c35d11165	Total value of the Jastrow	2023-03-30 18:17:33 +02:00
Anthony Scemama	0c136ab950	Renamed Jastrow into Jastrow_champ	2023-03-30 17:07:11 +02:00
Anthony Scemama	fdf6b905bb	Working on Jastrow	2023-03-30 12:34:17 +02:00
Anthony Scemama	2ecfc55dbc	Added eN cusp fitting in MOs	2023-03-17 14:54:58 +01:00
Anthony Scemama	e10c7584ff	Improved tensors in qmcalk_blas:	2023-03-16 16:52:01 +01:00
Anthony Scemama	c0131d5da4	Improved en_distance	2023-03-14 19:13:49 +01:00
Anthony Scemama	21336e0178	Working on e-n cusp	2023-03-14 14:59:14 +01:00
Anthony Scemama	4241461a20	Added ~ao_ang_mom~ and ~ao_nucl~	2023-03-13 17:06:41 +01:00
Anthony Scemama	3f33db6887	Transposed en_distance	2023-03-13 15:32:35 +01:00
Anthony Scemama	44c4c6c6d5	qmckl_mo_basis_select_mo is qmckl_exit_code, not bool	2023-03-09 11:01:55 +01:00
Anthony Scemama	c715f3e31f	Cleaning in Jastrow	2023-03-07 14:42:54 +01:00
Anthony Scemama	b269cd7403	Documentation	2023-03-05 23:38:34 +01:00
Anthony Scemama	da70f7a032	Fixed Fortran dependencies in parallel builds	2023-03-02 11:10:44 +01:00
Anthony Scemama	3ebb304218	Merge branch 'master' of github.com:TREX-CoE/qmckl	2023-03-02 10:10:49 +01:00
Anthony Scemama	71ea32ef2e	Cleaning Jastrow	2023-03-01 15:36:49 +01:00
Anthony Scemama	ea21ec2ef7	Removed GPU from Jastrow	2023-03-01 14:47:32 +01:00
Anthony Scemama	a051a1dd7d	Merge pull request #108 from fmgjcoppens/sm_splitting_fortran SM Splitting Fortran doc-kernel	2023-02-27 15:37:43 +01:00
Francois Coppens	b01c7c306b	Done with splitting doc.	2023-02-27 12:04:04 +01:00
Francois Coppens	1640eb60f9	Changed argument order.	2023-02-27 11:29:01 +01:00
Francois Coppens	656d268187	qmckl_sm_splitting_doc kernel works.	2023-02-26 12:34:10 +01:00
Francois Coppens	8e2674a3b2	reorder	2023-02-26 12:34:10 +01:00
Francois Coppens	7a97aa4a77	Fixed Fortran function call bug.	2023-02-26 12:34:10 +01:00
Francois Coppens	8216f682b3	Strange Fortran type error...	2023-02-26 12:34:10 +01:00
Francois Coppens	8ba882675e	Renamed function prefixes.	2023-02-26 12:34:10 +01:00
Francois Coppens	5e5c15a09d	- Added qmckl_context to Slagel Splitting kernel - Renamed it to Sherman-Morrison Splitting Core. - Sherman-Morrison Splitting Core now callable on its own. User is responsible for what to do with the output data. - Added default switch cases with asserts to generate crash with message if a template for a specific size is missing. - Added switch breaks to prevent the default case to always execute and make the kernel crash at the assert. - Reorganisded the Sherman-Morrison Splitting kernel so that the HPC variant always calls the Core HPC variant and not the generec variant and make duplicate decisions.	2023-02-26 12:34:10 +01:00
Anthony Scemama	0834c77946	Merge pull request #107 from fmgjcoppens/generate_private Generate private C-header / Fortran interface	2023-02-23 18:40:14 +01:00
Francois Coppens	7c57fe2b6f	Added fucntion that generates private fortran interfaces to C functions.	2023-02-23 17:30:12 +01:00
Francois Coppens	216fcebf70	Added fucntion that generates private c headers.	2023-02-23 16:17:39 +01:00
Francois Coppens	60a2d2c986	Fixes tab-key not working to open/close Org-mode sections in Emacs Evil mode. Fix works on Linux and macOS.	2023-02-23 15:03:33 +01:00
Anthony Scemama	1a0ea315ad	Merge pull request #102 from fmgjcoppens/master Fixed minor SIMD bug in tests.	2023-02-16 16:16:43 +01:00
Francois Coppens	b0f05b7c25	Commented out unused testdata.	2023-02-16 15:58:57 +01:00
Francois Coppens	1ee9635590	Added SM Splitting with doc version in Fortran skelleton plus Fortran/C interface.	2023-02-16 14:54:59 +01:00
Francois Coppens	181f662c68	Added macro HPC/DOC switch	2023-02-15 19:03:11 +01:00
Francois Coppens	4f0bdda4ff	...and the Fortran interfaces to the C-functions.	2023-02-15 18:49:12 +01:00
Francois Coppens	54a51b6ecc	Added Slagel splitting back + pedagogical skeleton function and interface.	2023-02-15 18:41:53 +01:00
Francois Coppens	c07553480c	Pedagogical Naive kernel works.	2023-02-15 11:46:48 +01:00
Anthony Scemama	6a637c394e	Merge pull request #103 from haampie/patch-1 bash -> sh	2023-02-14 13:46:31 +01:00
Anthony Scemama	4c27eb0078	Merge pull request #104 from haampie/patch-2 Use encoding in open	2023-02-14 13:33:39 +01:00
Harmen Stoppels	13f208165c	Use encoding in open Python's open is locale dependent, LC_ALL=C may open it in ascii mode and fail	2023-02-14 11:18:51 +01:00
Harmen Stoppels	eb8d8bf34e	bash -> sh	2023-02-14 11:04:04 +01:00
Francois Coppens	87d6acb49a	Adding documentation to ORG file.	2023-02-13 17:50:20 +01:00
Francois Coppens	42f4556fa3	Adding documentation to ORG file.	2023-02-13 17:49:18 +01:00
Francois Coppens	3482c832ac	Adding documentation to ORG file.	2023-02-13 17:48:31 +01:00
Francois Coppens	707fa17e09	Adding documentation to ORG file.	2023-02-13 17:44:11 +01:00
Francois Coppens	c0d4f766b1	Reorganising ORG file.	2023-02-13 15:08:37 +01:00
Francois Coppens	6ad4aabdfa	Still working	2023-02-10 17:16:08 +01:00
Francois Coppens	cc17b79316	Still working	2023-02-10 16:45:22 +01:00
Francois Coppens	06127f24cb	added return value to fortran interface.	2023-02-02 17:34:33 +01:00
Francois Coppens	8a89003bf2	Commented call to `_doc` kernel.	2023-02-02 17:23:14 +01:00
Francois Coppens	d3aebe52ff	Started adding the pedagogical kernels for the HAVE_DOC builds.	2023-02-02 17:04:34 +01:00
Francois Coppens	2e45927e04	Added AVX2 detection to autoconfig script. Fixed minor SIMD bug in tests.	2023-01-30 17:35:11 +01:00
Anthony Scemama	9a779f2a94	Avoid SIMD length=127	2023-01-30 16:25:02 +01:00
Anthony Scemama	30c3e48d91	Merge pull request #101 from fmgjcoppens/master Woodbury 2x2, 3x3 & Slagel Splitting kernel template generator	2023-01-30 08:56:53 +01:00
Francois Coppens	31ea30cdc3	Added Slagel Splitting kernel template generator.	2023-01-27 19:33:27 +01:00
Francois Coppens	6c0430a509	Added Woodbury 3x3 kernel template generator.	2023-01-27 17:41:32 +01:00
Francois Coppens	549413abca	Pulled out kernel template range so it can be set at the top, instead of at 8 different places throughout the code.	2023-01-27 15:24:52 +01:00
Francois Coppens	c58cf3c7f6	Added Woodbury 2x2 kernel template generator.	2023-01-27 14:31:25 +01:00
Francois Coppens	2d02b8cd63	Trivial rename	2023-01-27 14:31:25 +01:00
Anthony Scemama	345cf8525b	Merge pull request #100 from fmgjcoppens/master Return QMCKL_FAILURE if return code from qmckl_slagel_splitting equal…	2023-01-27 13:16:24 +01:00
Francois Coppens	5c0024f3f2	Return QMCKL_FAILURE if return code from qmckl_slagel_splitting equals QMCKL_FAILURE	2023-01-27 11:13:57 +01:00