Merge pull request #2 from TREX-CoE/main

Main

.github/workflows/gh-pages.yml

@@ -0,0 +1,34 @@
+name: github pages
+
+on:
+  push:
+    branches:
+      - master
+
+jobs:
+  deploy:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v2
+
+      - name: install extra repository
+        run: sudo add-apt-repository ppa:kelleyk/emacs
+
+      - name: refresh apt
+        run: sudo apt-get update
+
+      - name: install dependencies
+        run: sudo apt-get install emacs26
+
+      - name: install htmlize
+        run: git clone https://github.com/hniksic/emacs-htmlize && cp emacs-htmlize/htmlize.el docs/
+
+      - name: make
+        run: make -C src/ doc
+
+      - name: Deploy
+        uses: peaceiris/actions-gh-pages@v3
+        with:
+          github_token: ${{ secrets.GITHUB_TOKEN }}
+          publish_dir: ./docs
+

.github/workflows/test-build.yml

@@ -2,9 +2,9 @@ name: test-build
 
 on:
   push:
-    branches: [ main ]
+    branches: [ master ]
   pull_request:
-    branches: [ main ]
+    branches: [ master ]
 
 jobs:
   build:
@@ -13,7 +13,28 @@ jobs:
 
     steps:
     - uses: actions/checkout@v2
+
     - name: install dependencies
       run: sudo apt-get install emacs
+
     - name: make
       run: make -C src/
+
+
+  test:
+
+    runs-on: ubuntu-latest
+
+    steps:
+    - uses: actions/checkout@v2
+    - name: install dependencies
+      run: sudo apt-get install emacs
+    - name: Checkout submodules using a PAT
+      run: |
+           git config --file .gitmodules --get-regexp url | while read url; do
+             git config --file=.gitmodules $(echo "$url" | sed -E "s/git@github.com:|https:\/\/github.com\//https:\/\/${{ secrets.CI_PAT }}:${{ secrets.CI_PAT }}@github.com\//")
+           done
+           git submodule sync
+           git submodule update --init --recursive
+    - name: make
+      run: make -C src/ test

.gitmodules

@@ -0,0 +1,3 @@
+[submodule "munit"]
+	path = munit
+	url = https://github.com/nemequ/munit/

README.md

@@ -1,7 +1,28 @@
-# qmckl
+# QMCkl: Quantum Monte Carlo Kernel Library
 
 ![Build Status](https://github.com/TREX-CoE/qmckl/workflows/test-build/badge.svg?branch=main)
 
-Quantum Monte Carlo Kernel Library.
+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
+simple language and provide a standard API and tests to enable the
+development of high-performance QMCkl implementations taking
+advantage of modern hardware.
+
+See the [source code](https://github.com/TREX-CoE/qmckl/tree/main/src)
+to read the documentation.
+
+
+To clone the repository, use:
+```
+git clone --recursive https://github.com/TREX-CoE/qmckl.git
+```
+to dowload also the [munit](https://github.com/nemequ/munit) unit testing
+framework.
+
+
+------------------------------
+
+![European flag](https://trex-coe.eu/sites/default/files/inline-images/euflag.jpg)
+[TREX: Targeting Real Chemical Accuracy at the Exascale](https://trex-coe.eu) project has received funding from the European Union’s Horizon 2020 - Research and Innovation program - under grant agreement no. 952165. The content of this document does not represent the opinion of the European Union, and the European Union is not responsible for any use that might be made of such content.
 
-See the [Wiki](https://github.com/TREX-CoE/qmckl/wiki) for more information.

TODO.org

@@ -10,3 +10,9 @@ qmckl_malloc, where the domain id is something obtained from the
 context.
 
 
+* TRANSA, TRANSB
+* Performance info
+* Benchmark interpolation of basis functions
+* Complex numbers
+* Adjustable number for derivatives (1,2,3)
+

docs/config.el

@@ -0,0 +1,71 @@
+;; Thanks to Tobias's answer on Emacs Stack Exchange:
+;; https://emacs.stackexchange.com/questions/38437/org-mode-batch-export-missing-syntax-highlighting
+
+(package-initialize)
+(require 'htmlize)
+(require 'font-lock)
+(require 'subr-x) ;; for `when-let'
+
+(unless (boundp 'maximal-integer)
+  (defconst maximal-integer (lsh -1 -1)
+    "Maximal integer value representable natively in emacs lisp."))
+
+(defun face-spec-default (spec)
+  "Get list containing at most the default entry of face SPEC.
+Return nil if SPEC has no default entry."
+  (let* ((first (car-safe spec))
+     (display (car-safe first)))
+    (when (eq display 'default)
+      (list (car-safe spec)))))
+
+(defun face-spec-min-color (display-atts)
+  "Get min-color entry of DISPLAY-ATTS pair from face spec."
+  (let* ((display (car-safe display-atts)))
+    (or (car-safe (cdr (assoc 'min-colors display)))
+    maximal-integer)))
+
+(defun face-spec-highest-color (spec)
+  "Search face SPEC for highest color.
+That means the DISPLAY entry of SPEC
+with class 'color and highest min-color value."
+  (let ((color-list (cl-remove-if-not
+             (lambda (display-atts)
+               (when-let ((display (car-safe display-atts))
+                  (class (and (listp display)
+                          (assoc 'class display)))
+                  (background (assoc 'background display)))
+             (and (member 'light (cdr background))
+                  (member 'color (cdr class)))))
+             spec)))
+    (cl-reduce (lambda (display-atts1 display-atts2)
+         (if (> (face-spec-min-color display-atts1)
+            (face-spec-min-color display-atts2))
+             display-atts1
+           display-atts2))
+           (cdr color-list)
+           :initial-value (car color-list))))
+
+(defun face-spec-t (spec)
+  "Search face SPEC for fall back."
+  (cl-find-if (lambda (display-atts)
+        (eq (car-safe display-atts) t))
+          spec))
+
+(defun my-face-attribute (face attribute &optional frame inherit)
+  "Get FACE ATTRIBUTE from `face-user-default-spec' and not from `face-attribute'."
+  (let* ((face-spec (face-user-default-spec face))
+     (display-attr (or (face-spec-highest-color face-spec)
+               (face-spec-t face-spec)))
+     (attr (cdr display-attr))
+     (val (or (plist-get attr attribute) (car-safe (cdr (assoc attribute attr))))))
+    ;; (message "attribute: %S" attribute) ;; for debugging
+    (when (and (null (eq attribute :inherit))
+           (null val))
+      (let ((inherited-face (my-face-attribute face :inherit)))
+    (when (and inherited-face
+           (null (eq inherited-face 'unspecified)))
+      (setq val (my-face-attribute inherited-face attribute)))))
+    ;; (message "face: %S attribute: %S display-attr: %S, val: %S" face attribute display-attr val) ;; for debugging
+    (or val 'unspecified)))
+
+(advice-add 'face-attribute :override #'my-face-attribute)

munit

@@ -0,0 +1 @@
+Subproject commit fbbdf1467eb0d04a6ee465def2e529e4c87f2118

src/.gitignore

@@ -1,3 +1,11 @@
 *.o
 *.c
+*.f90
 *.h
+*.fh
+*.html
+*~
+*.so
+Makefile.generated
+test_qmckl
+merged_qmckl.org

src/Makefile

@@ -1,30 +1,66 @@
-CC=gcc
+COMPILER=GNU
-CFLAGS=-fexceptions -Wall -Werror -Wpedantic -Wextra
+#COMPILER=INTEL
+#COMPILER=LLVM
 
-FC=gfortran
+ifeq ($(COMPILER),GNU)
-FFLAGS=-fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant -Wuninitialized  -fbacktrace -ffpe-trap=zero,overflow,underflow -finit-real=nan
+CC=gcc -g
+CFLAGS=-fPIC -fexceptions -Wall -Werror -Wpedantic -Wextra
+
+FC=gfortran -g
+FFLAGS=-fPIC  -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant -Wuninitialized  -fbacktrace -ffpe-trap=zero,overflow,underflow -finit-real=nan
+
+LIBS=-lgfortran -lm
+endif
+
+ifeq ($(COMPILER),INTEL)
+CC=icc -xHost
+CFLAGS=-fPIC -g -O2
+
+FC=ifort -xHost
+FFLAGS=-fPIC -g -O2
+
+LIBS=-lm -lifcore -lirc
+endif
+
+#TODO
+ifeq ($(COMPILER),LLVM)
+CC=clang
+CFLAGS=-fPIC -g -O2
+
+FC=flang
+FFLAGS=fPIC -g -O2
+
+LIBS=-lm
+endif
 
 
-ORG_SOURCE_FILES=qmckl_context.org
+export CC CFLAGS FC FFLAGS LIBS
-OBJECT_FILES=$(patsubst %.org,%.o,$(ORG_SOURCE_FILES))
+
+MERGED_ORG=merged_qmckl.org
+ORG_SOURCE_FILES=$(wildcard *.org)
+OBJECT_FILES=$(filter-out $(EXCLUDED_OBJECTS), $(patsubst %.org,%.o,$(ORG_SOURCE_FILES)))
 
 .PHONY: clean
+.SECONDARY: # Needed to keep the produced C and Fortran files
 
-all: $(OBJECT_FILES)
+libqmckl.so: Makefile.generated 
+	$(MAKE) -f Makefile.generated
 
-%.c %.h: %.org
+test: Makefile.generated
-	emacs --quick --no-init-file --batch --eval "(require 'org)" --eval '(org-babel-tangle-file "$^")'
+	$(MAKE) -f Makefile.generated test
 
-%.c %.h %_f.f90: %.org
-	emacs --quick --no-init-file --batch --eval "(require 'org)" --eval '(org-babel-tangle-file "$^")'
 
-%.o: %.c
+doc: $(ORG_SOURCE_FILES)
-	$(CC) $(CFLAGS) -c $*.c -o $*.o
+	./merge_org.sh
+	./create_doc.sh $(MERGED_ORG)
+	rm $(MERGED_ORG)
 
-%.o: %.f90
-	$(FC) $(FFLAGS) -c $*.f90 -o $*.o
 
 clean:
-	rm -f qmckl_*.f90 qmckl_*.c qmckl_*.o qmckl_*.h
+	rm -f qmckl.h test_qmckl_* test_qmckl.c test_qmckl qmckl_*.f90 qmckl_*.c qmckl_*.o qmckl_*.h Makefile.generated libqmckl.so *.html *.fh *.mod
 
+Makefile.generated: Makefile create_makefile.sh  $(ORG_SOURCE_FILES)
+	./merge_org.sh
+	./create_makefile.sh $(MERGED_ORG) 
+	rm $(MERGED_ORG)
 

src/README.org

@@ -1,26 +1,51 @@
-* QMCkl source code
+#+TITLE: QMCkl source code documentation
+#+EXPORT_FILE_NAME: index.html
 
-** Introduction
+#+SETUPFILE: https://fniessen.github.io/org-html-themes/setup/theme-readtheorg.setup 
 
-   The main objective of present library is documentation. Therefore,
+* Introduction
-   literate programming is particularly adapted in this context.
+
-   Source files are written in org-mode format, to provide useful
+  The  ultimate  goal  of  QMCkl  is  to  provide  a  high-performance
+  implementation  of  the main  kernels  of  QMC. In  this  particular
+  repository, we focus on the definition of the API and the tests, and
+  on a /pedagogical/ presentation of  the algorithms.  We  expect the
+  HPC experts  to use  this repository as  a reference  for re-writing
+  optimized libraries.
+
+  Literate  programming  is  particularly  adapted  in  this  context.
+  Source  files are  written  in [[https://karl-voit.at/2017/09/23/orgmode-as-markup-only/][org-mode]]  format,  to provide  useful
   comments and LaTex formulas close to the code. There exists multiple
-   possibilities to convert org-mode files into different formats such as
+  possibilities to convert org-mode  files into different formats such
-   HTML or pdf.
+  as  HTML  or pdf.   For  a  tutorial  on literate  programming  with
-   For a tutorial on literate programming with org-mode, follow 
+  org-mode, follow [[http://www.howardism.org/Technical/Emacs/literate-programming-tutorial.html][this link]].
-   [[http://www.howardism.org/Technical/Emacs/literate-programming-tutorial.html][this link]].
 
-   The code is extracted from the org files using Emacs as a command-line
+  The  code  is  extracted  from  the  org  files  using  Emacs  as  a
-   tool in the =Makefile=, and then the produced files are compiled.
+  command-line tool in the =Makefile=, and then the produced files are
+  compiled.
 
-*** Source code editing
+** Language used
+
+   Fortran is one of the most  common languages used by the community,
+   and  is simple  enough to  make the  algorithms readable.  Hence we
+   propose in this pedagogical implementation  of QMCkl to use Fortran
+   to express  the algorithms.  For specific internal  functions where
+   the C language is more natural, C is used.
+
+   As Fortran  modules generate  compiler-dependent files, the  use of
+   modules is restricted to the internal use of the library, otherwise
+   the compliance with C is violated.
+
+   The external dependencies  should be kept as small  as possible, so
+   external libraries should be used /only/ if their used is strongly
+   justified.
+
+** Source code editing
 
    Any text  editor can be used  to edit org-mode files.  For a better
-    user experience Emacs is recommended.
+   user experience Emacs  is recommended.  For users  hating Emacs, it
-    For users hating Emacs, it is good to know that Emacs can behave
+   is good to  know that Emacs can behave like  Vim when switched into
-    like Vim when switched into ``Evil'' mode. There also exists
+   ``Evil''  mode.  There  also   exists  [[https://www.spacemacs.org][Spacemacs]]  which  helps  the
-    [[https://www.spacemacs.org][Spacemacs]] which is particularly well adapted to Vim users.
+   transition for Vim users.
 
    For users  with a preference  for Jupyter notebooks,  the following
    script can convert jupyter notebooks to org-mode files:
@@ -40,14 +65,143 @@ rm ${nb}.md
 
    And pandoc can convert multiple markdown formats into org-mode.
 
-*** Writing in Fortran
+** Writing in Fortran
 
    The  Fortran  source  files  should provide  a  C  interface  using
-    iso-c-binding. The name of the Fortran source files should end
+ =iso_c_binding=. The  name of the  Fortran source files  should end
-    with =_f.f90= to be properly handled by the Makefile.
+   with =_f.f90= to be properly handled by the Makefile.  The names of
+   the  functions defined  in  fortran  should be  the  same as  those
+   exposed  in the  API  suffixed by =_f=.   Fortran interface  files
+   should also be written in the =qmckl_f.f90= file.
    
-** Documentation
+   For more guidelines on using Fortran to generate a C interface, see
+   [[http://fortranwiki.org/fortran/show/Generating+C+Interfaces][this link]].
 
-- [[qmckl_context.org][Context]]
+** Coding style
+   # TODO: decide on a coding style
+
+   To improve  readability, we maintain  a consistent coding  style in
+   the library.
+
+   - For C source files, we will use __(decide on a coding style)__
+   - For  Fortran source  files, we  will  use __(decide  on a  coding
+     style)__
+
+   Coding style can be automatically checked with [[https://clang.llvm.org/docs/ClangFormat.html][clang-format]].
+
+** Design of the library
+
+   The proposed API should allow the library to:
+   - deal with memory transfers between CPU and accelerators
+   - use different levels of floating-point precision
+
+   We  chose  a multi-layered  design  with  low-level and  high-level
+   functions (see below).
+
+*** Naming conventions
+
+    Use =qmckl_= as a prefix for all exported functions and variables.
+    All exported header  files should have a filename  with the prefix
+ =qmckl_=.
+
+    If the  name of the  org-mode file is =xxx.org=, the name  of the
+    produced C files should be =xxx.c= and =xxx.h= and the name of the
+    produced Fortran files should be =xxx.f90=
+   
+    Arrays are in uppercase and scalars are in lowercase.
+    
+    In the  names of  the variables and  functions, only  the singular
+    form is allowed.
+
+*** Application programming interface
+
+    The  application programming  interface  (API) is  designed to  be
+    compatible with  the C programming  language (not C++),  to ensure
+    that the  library will be  easily usable in /any/  language.  This
+    implies that only the following data types are allowed in the API:
+
+    - 32-bit and 64-bit floats and arrays (=real= and =double=)
+    - 32-bit and 64-bit integers and arrays (=int32_t= and =int64_t=)
+    - Pointers  should  be represented  as  64-bit  integers (even  on
+      32-bit architectures)
+    - ASCII  strings are  represented  as a  pointers  to a  character
+      arrays and terminated by a zero character (C convention).
+   
+    Complex numbers can be represented by an array of 2 floats.
+
+    # TODO : Link to repositories for bindings
+    To facilitate the  use in other languages than C,  we provide some
+    bindings in other languages in other repositories.
+
+*** Global state
+
+    Global variables should  be avoided in the library,  because it is
+    possible that one  single program needs to  use multiple instances
+    of the library. To solve this  problem we propose to use a pointer
+    to  a =context=  variable,  built   by  the  library   with  the
+ =qmckl_context_create= function. The =context= contains the global
+    state of  the library, and is  used as the first  argument of many
+    QMCkl functions.
+
+    The internal structure of the context  is not specified, to give a
+    maximum of  freedom to  the different  implementations.  Modifying
+    the  state   is  done   by  setters   and  getters,   prefixed  by
+ =qmckl_context_set_=  an =qmckl_context_get_=.   When  a  context
+    variable is modified by a setter, a copy of the old data structure
+    is made and updated, and the  pointer to the new data structure is
+    returned, such that the old contexts can still be accessed.  It is
+    also possible to modify the state  in an impure fashion, using the
+ =qmckl_context_update_=  functions.   The   context  and  its  old
+    versions can be destroyed with =qmckl_context_destroy=.
+
+*** Low-level functions
+
+    Low-level functions are very simple  functions which are leaves of
+    the function call tree (they don't call any other QMCkl function).
+
+    These  functions   are /pure/,   and  unaware   of   the   QMCkl
+ =context=. They are not allowed to allocate/deallocate memory, and
+    if they need temporary memory it should be provided in input.
+
+*** High-level functions
+
+    High-level functions  are at  the top of  the function  call tree.
+    They  are  able  to  choose which  lower-level  function  to  call
+    depending on the required precision, and do the corresponding type
+    conversions.  These functions are  also responsible for allocating
+    temporary storage, to simplify the use of accelerators.
+
+    The high-level  functions should be pure,  unless the introduction
+    of non-purity is justified. All the side effects should be made in
+    the =context= variable.
+
+    # TODO : We need an identifier for impure functions
+
+*** Numerical precision
+
+    The number of bits of precision  required for a function should be
+    given as an input of low-level computational functions. This input
+    will be used to define the values of the different thresholds that
+    might be  used to  avoid computing unnecessary  noise.  High-level
+    functions  will  use  the  precision specified  in  the =context=
+    variable.
+
+** Algorithms
+  
+   Reducing the scaling of an  algorithm usually implies also reducing
+   its arithmetic  complexity (number  of flops per  byte). Therefore,
+   for  small  sizes   \(\mathcal{O}(N^3)\)  and  \(\mathcal{O}(N^2)\)
+   algorithms are  better adapted than linear  scaling algorithms.  As
+   QMCkl is a  general purpose library, multiple  algorithms should be
+   implemented adapted to different problem sizes.
+
+** Rules for the API
+   
+   - =stdint= should be used for integers (=int32_t=, =int64_t=)
+   - integers used for counting should always be =int64_t=
+   - floats should be by default =double=, unless explicitly mentioned 
+   - pointers are converted to =int64_t= to increase portability
+
+* Documentation
 
   

src/create_doc.sh

@@ -0,0 +1,13 @@
+#!/bin/bash
+INPUT=$1
+
+if [[ -f ../docs/htmlize.el ]]
+then
+    emacs --batch --load ../docs/htmlize.el --load ../docs/config.el $INPUT -f org-html-export-to-html
+else
+    emacs --batch --load ../docs/config.el $INPUT -f org-html-export-to-html
+fi
+
+mv index.html ../docs
+
+

src/create_makefile.sh

@@ -0,0 +1,95 @@
+#!/bin/bash
+
+INPUT=$1
+OUTPUT=Makefile.generated
+
+# Tangle org files
+emacs \
+     $INPUT \
+     --batch \
+     -f org-babel-tangle \
+     --kill
+
+
+
+# Create the list of *.o files to be created
+
+OBJECTS=""
+for i in $(ls qmckl_*.c) ; do
+    FILE=${i%.c}
+    OBJECTS="${OBJECTS} ${FILE}.o"
+done >> $OUTPUT
+
+for i in $(ls qmckl_*.f90) ; do
+    FILE=${i%.f90}
+    OBJECTS="${OBJECTS} ${FILE}.o"
+done >> $OUTPUT
+
+TESTS=""
+for i in $(ls test_qmckl_*.c) ; do
+    FILE=${i%.c}.o
+    TESTS="${TESTS} ${FILE}"
+done >> $OUTPUT
+
+TESTS_F=""
+for i in $(ls test_qmckl_*.f90) ; do
+    FILE=${i%.f90}.o
+    TESTS_F="${TESTS_F} ${FILE}"
+done >> $OUTPUT
+
+
+# Write the Makefile
+
+cat << EOF > $OUTPUT
+CC=$CC
+CFLAGS=$CFLAGS -I../munit/ 
+
+FC=$FC
+FFLAGS=$FFLAGS
+OBJECT_FILES=$OBJECTS
+TESTS=$TESTS
+TESTS_F=$TESTS_F
+
+LIBS=$LIBS
+
+libqmckl.so: \$(OBJECT_FILES)
+	\$(CC) -shared \$(OBJECT_FILES) -o libqmckl.so
+
+%.o: %.c 
+	\$(CC) \$(CFLAGS) -c \$*.c -o \$*.o
+
+%.o: %.f90 qmckl_f.o
+	\$(FC) \$(FFLAGS) -c \$*.f90 -o \$*.o
+
+test_qmckl: test_qmckl.c libqmckl.so \$(TESTS) \$(TESTS_F)
+	\$(CC) \$(CFLAGS) -Wl,-rpath,$PWD -L. \
+	../munit/munit.c \$(TESTS) \$(TESTS_F) -lqmckl \$(LIBS) test_qmckl.c -o test_qmckl
+
+test: test_qmckl
+	./test_qmckl
+
+.PHONY: test
+EOF
+
+for i in $(ls qmckl_*.c) ; do
+    FILE=${i%.c}
+    echo "${FILE}.o: ${FILE}.c " *.h
+done >> $OUTPUT
+
+for i in $(ls qmckl_*.f90) ; do
+    FILE=${i%.f90}
+    echo "${FILE}.o: ${FILE}.f90"
+done >> $OUTPUT
+
+for i in $(ls test_qmckl_*.c) ; do
+    FILE=${i%.c}
+    echo "${FILE}.o: ${FILE}.c qmckl.h" 
+done >> $OUTPUT
+
+
+for i in $(ls test_qmckl*.f90) ; do
+    FILE=${i%.f90}
+    echo "${FILE}.o: ${FILE}.f90"
+done >> $OUTPUT
+
+

src/merge_org.sh

@@ -0,0 +1,13 @@
+#!/bin/bash
+
+for i in README.org \
+             qmckl.org \
+             qmckl_memory.org \
+             qmckl_context.org \
+             qmckl_distance.org \
+             qmckl_ao.org \
+             qmckl_footer.org \
+             test_qmckl.org
+do
+    cat $i >> merged_qmckl.org
+done

src/qmckl.org

@@ -0,0 +1,64 @@
+** =qmckl.h= header file
+
+ This file produces the =qmckl.h= header file, which is to be included
+ when qmckl functions are used.
+
+ We also create here the =qmckl_f.f90= which is the Fortran interface file.
+
+*** Top of header files                                            :noexport:
+  
+    #+BEGIN_SRC C :tangle qmckl.h
+#ifndef QMCKL_H
+#define QMCKL_H
+#include <stdlib.h>
+#include <stdint.h>
+#include <math.h>
+    #+END_SRC
+
+    #+BEGIN_SRC f90 :tangle qmckl_f.f90
+module qmckl
+  use, intrinsic :: iso_c_binding
+    #+END_SRC
+
+    The bottoms of the files are located in the [[qmckl_footer.org]] file.
+   
+*** Constants
+
+**** Success/failure
+
+     These are the codes returned by the functions to indicate success
+     or failure. All such functions should have as a return type =qmckl_exit_code=.
+
+      #+BEGIN_SRC C :comments org :tangle qmckl.h
+#define QMCKL_SUCCESS 0
+#define QMCKL_FAILURE 1
+
+typedef int32_t qmckl_exit_code;
+typedef int64_t qmckl_context ;
+
+      #+END_SRC
+
+  #+BEGIN_SRC f90 :comments org :tangle qmckl_f.f90
+integer, parameter :: QMCKL_SUCCESS = 0
+integer, parameter :: QMCKL_FAILURE = 0
+  #+END_SRC
+
+**** Precision-related constants
+
+     Controlling numerical precision enables optimizations. Here, the
+     default parameters determining the target numerical precision and
+     range are defined.
+
+      #+BEGIN_SRC C :comments org :tangle qmckl.h
+#define QMCKL_DEFAULT_PRECISION 53
+#define QMCKL_DEFAULT_RANGE     11
+      #+END_SRC
+
+  #+BEGIN_SRC f90 :comments org :tangle qmckl_f.f90
+integer, parameter :: QMCKL_DEFAULT_PRECISION = 53
+integer, parameter :: QMCKL_DEFAULT_RANGE = 11
+  #+END_SRC
+ 
+
+  # -*- mode: org -*-
+  # vim: syntax=c

src/qmckl_ao.org

@@ -0,0 +1,742 @@
+** Atomic Orbitals
+
+
+ This  files contains  all the  routines  for the  computation of  the
+ values, gradients and Laplacian of the atomic basis functions.
+
+ 3 files are produced:
+ - a source file : =qmckl_ao.f90=
+ - a C test file : =test_qmckl_ao.c=
+ - a Fortran test file : =test_qmckl_ao_f.f90=
+
+*** Test                                                           :noexport:
+      #+BEGIN_SRC C :tangle test_qmckl_ao.c
+#include "qmckl.h"
+#include "munit.h"
+MunitResult test_qmckl_ao() {
+  qmckl_context context;
+  context = qmckl_context_create();
+      #+END_SRC
+
+*** Polynomials
+
+    \[
+    P_l(\mathbf{r},\mathbf{R}_i)  =   (x-X_i)^a (y-Y_i)^b (z-Z_i)^c 
+    \]
+    \begin{eqnarray*} 
+    \frac{\partial }{\partial x} P_l\left(\mathbf{r},\mathbf{R}_i \right) &
+                   = & a (x-X_i)^{a-1} (y-Y_i)^b (z-Z_i)^c \\
+    \frac{\partial }{\partial y} P_l\left(\mathbf{r},\mathbf{R}_i \right) &
+                   = & b (x-X_i)^a (y-Y_i)^{b-1} (z-Z_i)^c \\
+    \frac{\partial }{\partial z} P_l\left(\mathbf{r},\mathbf{R}_i \right) &
+                   = & c (x-X_i)^a (y-Y_i)^b (z-Z_i)^{c-1} \\
+    \end{eqnarray*} 
+    \begin{eqnarray*} 
+    \left( \frac{\partial }{\partial x^2} + 
+               \frac{\partial }{\partial y^2} + 
+               \frac{\partial }{\partial z^2} \right) P_l
+               \left(\mathbf{r},\mathbf{R}_i \right) &  = &  
+             a(a-1) (x-X_i)^{a-2} (y-Y_i)^b (z-Z_i)^c + \\
+          && b(b-1) (x-X_i)^a (y-Y_i)^{b-1} (z-Z_i)^c + \\
+          && c(c-1) (x-X_i)^a (y-Y_i)^b (z-Z_i)^{c-1} 
+    \end{eqnarray*}
+
+**** =qmckl_ao_power=
+
+     Computes all  the powers of  the =n= input  data up to  the given
+     maximum value given in input for each of the $n$ points:
+   
+     \[ P_{ij} = X_j^i \]
+
+***** Arguments
+
+      | =context=  | input  | Global state                                      |
+      | =n=        | input  | Number of values                                  |
+      | =X(n)=     | input  | Array containing the input values                 |
+      | =LMAX(n)=  | input  | Array containing the maximum power for each value |
+      | =P(LDP,n)= | output | Array containing all the powers of =X=            |
+      | =LDP=      | input  | Leading dimension of array =P=                    |
+    
+***** Requirements
+
+      - =context= is not 0
+      - =n= > 0
+      - =X= is allocated with at least $n \times 8$ bytes
+      - =LMAX= is allocated with at least $n \times 4$ bytes
+      - =P= is allocated with at least $n \times \max_i \text{LMAX}_i \times 8$ bytes
+      - =LDP= >= $\max_i$ =LMAX[i]=
+
+***** Header
+      #+BEGIN_SRC C :tangle qmckl.h
+qmckl_exit_code qmckl_ao_power(const qmckl_context context,
+                const int64_t n, 
+                const double *X, const int32_t *LMAX,
+                const double *P, const int64_t LDP);
+      #+END_SRC
+    
+***** Source
+      #+BEGIN_SRC f90 :tangle qmckl_ao.f90
+integer function qmckl_ao_power_f(context, n, X, LMAX, P, ldp) result(info)
+  implicit none
+  integer*8 , intent(in)  :: context
+  integer*8 , intent(in)  :: n
+  real*8    , intent(in)  :: X(n)
+  integer   , intent(in)  :: LMAX(n)
+  real*8    , intent(out) :: P(ldp,n)
+  integer*8 , intent(in)  :: ldp
+
+  integer*8  :: i,j
+
+  info = 0
+
+  if (context == 0_8) then
+     info = -1
+     return
+  endif
+  
+  if (LDP < MAXVAL(LMAX)) then
+     info = -2
+     return
+  endif
+  
+  do j=1,n
+    P(1,j) = X(j)
+    do i=2,LMAX(j)
+       P(i,j) = P(i-1,j) * X(j) 
+    end do
+  end do
+
+end function qmckl_ao_power_f
+      #+END_SRC 
+
+***** C interface                                                  :noexport:
+      #+BEGIN_SRC f90 :tangle qmckl_ao.f90
+integer(c_int32_t) function qmckl_ao_power(context, n, X, LMAX, P, ldp) &
+     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 :: n
+  real    (c_double)  , intent(in)         :: X(n)
+  integer (c_int32_t) , intent(in)         :: LMAX(n)
+  real    (c_double)  , intent(out)        :: P(ldp,n)
+  integer (c_int64_t) , intent(in) , value :: ldp
+  
+  integer, external :: qmckl_ao_power_f
+  info = qmckl_ao_power_f(context, n, X, LMAX, P, ldp)
+end function qmckl_ao_power
+      #+END_SRC
+
+      #+BEGIN_SRC f90 :tangle qmckl_f.f90
+  interface
+     integer(c_int32_t) function qmckl_ao_power(context, n, X, LMAX, P, ldp) bind(C)
+       use, intrinsic :: iso_c_binding
+       integer (c_int64_t) , intent(in) , value :: context
+       integer (c_int64_t) , intent(in) , value :: n
+       integer (c_int64_t) , intent(in) , value :: ldp
+       real    (c_double)  , intent(in)         :: X(n)
+       integer (c_int32_t) , intent(in)         :: LMAX(n)
+       real    (c_double)  , intent(out)        :: P(ldp,n)
+     end function qmckl_ao_power
+  end interface
+      #+END_SRC
+  
+***** Test                                                         :noexport:
+      #+BEGIN_SRC f90 :tangle test_qmckl_ao_f.f90
+integer(c_int32_t) function test_qmckl_ao_power(context) bind(C)
+  use qmckl
+  implicit none
+
+  integer(c_int64_t), intent(in), value :: context
+  
+  integer*8                     :: n, LDP 
+  integer, allocatable          :: LMAX(:) 
+  double precision, allocatable :: X(:), P(:,:)
+  integer*8                     :: i,j
+  double precision              :: epsilon
+
+  epsilon = qmckl_context_get_epsilon(context)
+
+  n = 100;
+  LDP = 10;
+  
+  allocate(X(n), P(LDP,n), LMAX(n))
+  
+  do j=1,n
+     X(j) = -5.d0 + 0.1d0 * dble(j)
+     LMAX(j) = 1 + int(mod(j, 5),4)
+  end do
+  
+  test_qmckl_ao_power = qmckl_ao_power(context, n, X, LMAX, P, LDP) 
+  if (test_qmckl_ao_power /= 0) return
+  
+  test_qmckl_ao_power = -1
+  
+  do j=1,n
+     do i=1,LMAX(j)
+        if ( X(j)**i == 0.d0 ) then
+           if ( P(i,j) /= 0.d0) return
+        else
+           if ( dabs(1.d0 - P(i,j) / (X(j)**i)) > epsilon ) return
+        end if
+     end do
+  end do
+
+  test_qmckl_ao_power = 0
+  deallocate(X,P,LMAX)
+end function test_qmckl_ao_power
+      #+END_SRC
+
+      #+BEGIN_SRC C :tangle test_qmckl_ao.c
+int  test_qmckl_ao_power(qmckl_context context);
+munit_assert_int(0, ==, test_qmckl_ao_power(context));
+      #+END_SRC
+  
+
+**** =qmckl_ao_polynomial_vgl=
+   
+     Computes the values, gradients and Laplacians at a given point of
+     all polynomials with an angular momentum up to =lmax=.
+
+***** Arguments
+
+      | =context=    | input  | Global state                                         |
+      | =X(3)=       | input  | Array containing the coordinates of the points       |
+      | =R(3)=       | input  | Array containing the x,y,z coordinates of the center |
+      | =lmax=       | input  | Maximum angular momentum                             |
+      | =n=          | output | Number of computed polynomials                       |
+      | =L(ldl,n)=   | output | Contains a,b,c for all =n= results                   |
+      | =ldl=        | input  | Leading dimension of =L=                             |
+      | =VGL(ldv,n)= | output | Value, gradients and Laplacian of the polynomials    |
+      | =ldv=        | input  | Leading dimension of array =VGL=                     |
+    
+***** Requirements
+
+      - =context= is not 0
+      - =n= > 0
+      - =lmax= >= 0
+      - =ldl= >= 3
+      - =ldv= >= 5
+      - =X= is allocated with at least $3 \times 8$ bytes
+      - =R= is allocated with at least $3 \times 8$ bytes
+      - =n= >= =(lmax+1)(lmax+2)(lmax+3)/6=
+      - =L= is allocated with at least $3 \times n \times 4$ bytes
+      - =VGL= is allocated with at least $5 \times n \times 8$ bytes
+      - On output, =n= should be equal to =(lmax+1)(lmax+2)(lmax+3)/6=
+      - On output, the powers are given in the following order (l=a+b+c):
+        - Increase values of =l=
+        - Within a given value of =l=, alphabetical order of the
+          string made by a*"x" + b*"y" + c*"z" (in Python notation).
+          For example, with a=0, b=2 and c=1 the string is "yyz"
+        
+***** Error codes
+
+      | -1 | Null context        |
+      | -2 | Inconsistent =ldl=  |
+      | -3 | Inconsistent =ldv=  |
+      | -4 | Inconsistent =lmax= |
+
+***** Header
+      #+BEGIN_SRC C :tangle qmckl.h
+qmckl_exit_code qmckl_ao_polynomial_vgl(const qmckl_context context,
+                const double *X, const double *R,
+                const int32_t lmax, const int64_t *n,
+                const int32_t *L,   const int64_t ldl,
+                const double *VGL,  const int64_t ldv);
+      #+END_SRC
+    
+***** Source
+      #+BEGIN_SRC f90 :tangle qmckl_ao.f90
+integer function qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL, ldv) result(info)
+  implicit none
+  integer*8 , intent(in)  :: context
+  real*8    , intent(in)  :: X(3), R(3)
+  integer   , intent(in)  :: lmax
+  integer*8 , intent(out) :: n
+  integer   , intent(out) :: L(ldl,(lmax+1)*(lmax+2)*(lmax+3)/6)
+  integer*8 , intent(in)  :: ldl
+  real*8    , intent(out) :: VGL(ldv,(lmax+1)*(lmax+2)*(lmax+3)/6)
+  integer*8 , intent(in)  :: ldv
+
+  integer*8         :: i,j
+  integer           :: a,b,c,d
+  real*8            :: Y(3)
+  integer           :: lmax_array(3)
+  real*8            :: pows(-2:lmax,3)
+  integer, external :: qmckl_ao_power_f
+  double precision  :: xy, yz, xz
+  double precision  :: da, db, dc, dd
+  
+  info = 0
+  
+  if (context == 0_8) then
+     info = -1
+     return
+  endif
+  
+  if (ldl < 3) then
+     info = -2
+     return
+  endif
+  
+  if (ldv < 5) then
+     info = -3
+     return
+  endif
+  
+  if (lmax <= 0) then
+     info = -4
+     return
+  endif
+  
+  
+  do i=1,3
+     Y(i) = X(i) - R(i)
+  end do
+
+  lmax_array(1:3) = lmax
+  if (lmax == 0) then
+    VGL(1,1) = 1.d0
+    vgL(2:5,1) = 0.d0
+    l(1:3,1) = 0
+    n=1
+  else if (lmax > 0) then
+    pows(-2:0,1:3) = 1.d0
+    do i=1,lmax
+        pows(i,1) = pows(i-1,1) * Y(1) 
+        pows(i,2) = pows(i-1,2) * Y(2) 
+        pows(i,3) = pows(i-1,3) * Y(3) 
+    end do
+
+    VGL(1:5,1:4) = 0.d0
+    l(1:3,1:4) = 0
+
+    VGL(1,1) = 1.d0
+    vgl(1:5,2:4) = 0.d0
+
+    l(1,2) = 1
+    vgl(1,2) = pows(1,1)
+    vgL(2,2) = 1.d0
+
+    l(2,3) = 1
+    vgl(1,3) = pows(1,2)
+    vgL(3,3) = 1.d0
+
+    l(3,4) = 1
+    vgl(1,4) = pows(1,3)
+    vgL(4,4) = 1.d0
+
+    n=4
+  endif
+        
+  ! l>=2
+  dd = 2.d0
+  do d=2,lmax
+     da = dd
+     do a=d,0,-1
+        db = dd-da
+        do b=d-a,0,-1
+           c  = d  - a  - b
+           dc = dd - da - db
+           n = n+1
+
+           l(1,n) = a
+           l(2,n) = b
+           l(3,n) = c
+           
+           xy = pows(a,1) * pows(b,2)
+           yz = pows(b,2) * pows(c,3)
+           xz = pows(a,1) * pows(c,3)
+           
+           vgl(1,n) = xy * pows(c,3)
+           
+           xy = dc * xy
+           xz = db * xz
+           yz = da * yz
+           
+           vgl(2,n) = pows(a-1,1) * yz
+           vgl(3,n) = pows(b-1,2) * xz
+           vgl(4,n) = pows(c-1,3) * xy
+           
+           vgl(5,n) = &
+                (da-1.d0) * pows(a-2,1) * yz + &
+                (db-1.d0) * pows(b-2,2) * xz + &
+                (dc-1.d0) * pows(c-2,3) * xy
+
+           db = db - 1.d0
+        end do
+        da = da - 1.d0
+     end do
+     dd = dd + 1.d0
+  end do
+
+  info = 0
+
+end function qmckl_ao_polynomial_vgl_f
+      #+END_SRC 
+
+***** C interface                                                  :noexport:
+      #+BEGIN_SRC f90 :tangle qmckl_ao.f90
+integer(c_int32_t) function qmckl_ao_polynomial_vgl(context, X, R, lmax, n, L, ldl, VGL, ldv) &
+     bind(C) result(info)
+  use, intrinsic :: iso_c_binding
+  implicit none
+  integer (c_int64_t) , intent(in) , value :: context
+  real    (c_double)  , intent(in)         :: X(3), R(3)
+  integer (c_int32_t) , intent(in) , value :: lmax
+  integer (c_int64_t) , intent(out)        :: n
+  integer (c_int32_t) , intent(out)        :: L(ldl,(lmax+1)*(lmax+2)*(lmax+3)/6)
+  integer (c_int64_t) , intent(in) , value :: ldl
+  real    (c_double)  , intent(out)        :: VGL(ldv,(lmax+1)*(lmax+2)*(lmax+3)/6)
+  integer (c_int64_t) , intent(in) , value :: ldv
+
+  integer, external :: qmckl_ao_polynomial_vgl_f
+  info = qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL, ldv) 
+end function qmckl_ao_polynomial_vgl
+      #+END_SRC
+
+***** Fortran interface                                            :noexport:
+      #+BEGIN_SRC f90 :tangle qmckl_f.f90
+  interface
+     integer(c_int32_t) function qmckl_ao_polynomial_vgl(context, X, R, lmax, n, L, ldl, VGL, ldv) &
+          bind(C) 
+       use, intrinsic :: iso_c_binding
+       integer (c_int64_t) , intent(in) , value :: context
+       integer (c_int32_t) , intent(in) , value :: lmax
+       integer (c_int64_t) , intent(in) , value :: ldl
+       integer (c_int64_t) , intent(in) , value :: ldv
+       real    (c_double)  , intent(in)         :: X(3), R(3)
+       integer (c_int64_t) , intent(out)        :: n
+       integer (c_int32_t) , intent(out)        :: L(ldl,(lmax+1)*(lmax+2)*(lmax+3)/6)
+       real    (c_double)  , intent(out)        :: VGL(ldv,(lmax+1)*(lmax+2)*(lmax+3)/6)
+     end function qmckl_ao_polynomial_vgl
+  end interface
+      #+END_SRC
+***** Test                                                         :noexport:
+      #+BEGIN_SRC f90 :tangle test_qmckl_ao_f.f90
+integer(c_int32_t) function test_qmckl_ao_polynomial_vgl(context) bind(C)
+  use qmckl
+  implicit none
+
+  integer(c_int64_t), intent(in), value :: context
+  
+  integer                       :: lmax, d, i
+  integer, allocatable          :: L(:,:)
+  integer*8                     :: n, ldl, ldv, j
+  double precision              :: X(3), R(3), Y(3)
+  double precision, allocatable :: VGL(:,:)
+  double precision              :: w
+  double precision              :: epsilon
+
+  epsilon = qmckl_context_get_epsilon(context)
+
+  X = (/ 1.1 , 2.2 ,  3.3 /)
+  R = (/ 0.1 , 1.2 , -2.3 /)
+  Y(:) = X(:) - R(:)
+
+  lmax = 4;
+  n = 0;
+  ldl = 3;
+  ldv = 100;
+  
+  d = (lmax+1)*(lmax+2)*(lmax+3)/6
+
+  allocate (L(ldl,d), VGL(ldv,d))
+
+  test_qmckl_ao_polynomial_vgl = &
+       qmckl_ao_polynomial_vgl(context, X, R, lmax, n, L, ldl, VGL, ldv)
+  if (test_qmckl_ao_polynomial_vgl /= 0) return
+
+  test_qmckl_ao_polynomial_vgl = -1
+
+  if (n /= d) return 
+
+  do j=1,n
+     test_qmckl_ao_polynomial_vgl = -11
+     do i=1,3
+        if (L(i,j) < 0) return
+     end do
+     test_qmckl_ao_polynomial_vgl = -12
+     if (dabs(1.d0 - VGL(1,j) / (&
+          Y(1)**L(1,j) * Y(2)**L(2,j) * Y(3)**L(3,j)  &
+          )) > epsilon ) return
+
+     test_qmckl_ao_polynomial_vgl = -13
+     if (L(1,j) < 1) then
+        if (VGL(2,j) /= 0.d0) return
+     else
+        if (dabs(1.d0 - VGL(2,j) / (&
+             L(1,j) * Y(1)**(L(1,j)-1) * Y(2)**L(2,j) * Y(3)**L(3,j) &
+             )) > epsilon ) return
+     end if
+
+     test_qmckl_ao_polynomial_vgl = -14
+     if (L(2,j) < 1) then
+        if (VGL(3,j) /= 0.d0) return
+     else
+        if (dabs(1.d0 - VGL(3,j) / (&
+             L(2,j) * Y(1)**L(1,j) * Y(2)**(L(2,j)-1) * Y(3)**L(3,j) &
+             )) > epsilon ) return
+     end if
+
+     test_qmckl_ao_polynomial_vgl = -15
+     if (L(3,j) < 1) then
+        if (VGL(4,j) /= 0.d0) return
+     else
+        if (dabs(1.d0 - VGL(4,j) / (&
+             L(3,j) * Y(1)**L(1,j) * Y(2)**L(2,j) * Y(3)**(L(3,j)-1) &
+             )) > epsilon ) return
+     end if
+     
+     test_qmckl_ao_polynomial_vgl = -16
+     w = 0.d0
+     if (L(1,j) > 1) then
+        w = w + L(1,j) * (L(1,j)-1) * Y(1)**(L(1,j)-2) * Y(2)**L(2,j) * Y(3)**L(3,j) 
+     end if
+     if (L(2,j) > 1) then
+        w = w + L(2,j) * (L(2,j)-1) * Y(1)**L(1,j) * Y(2)**(L(2,j)-2) * Y(3)**L(3,j) 
+     end if
+     if (L(3,j) > 1) then
+        w = w + L(3,j) * (L(3,j)-1) * Y(1)**L(1,j) * Y(2)**L(2,j) * Y(3)**(L(3,j)-2) 
+     end if
+     if (dabs(1.d0 - VGL(5,j) / w) > epsilon ) return
+  end do
+
+  test_qmckl_ao_polynomial_vgl = 0
+     
+  deallocate(L,VGL)
+end function test_qmckl_ao_polynomial_vgl
+      #+END_SRC
+
+      #+BEGIN_SRC C :tangle test_qmckl_ao.c
+int  test_qmckl_ao_polynomial_vgl(qmckl_context context);
+munit_assert_int(0, ==, test_qmckl_ao_polynomial_vgl(context));
+      #+END_SRC
+      #+END_SRC
+   
+*** Gaussian basis functions
+
+**** =qmckl_ao_gaussian_vgl=
+   
+     Computes the values, gradients and Laplacians at a given point of
+ =n= Gaussian functions centered at the same point:
+   
+     \[ v_i = exp(-a_i |X-R|^2) \]
+     \[ \nabla_x v_i = -2 a_i (X_x -  R_x) v_i \]
+     \[ \nabla_y v_i = -2 a_i (X_y -  R_y) v_i \]
+     \[ \nabla_z v_i = -2 a_i (X_z -  R_z) v_i \]
+     \[ \Delta v_i = a_i (4 |X-R|^2 a_i - 6) v_i \]
+
+***** Arguments
+
+      | =context=    | input  | Global state                                         |
+      | =X(3)=       | input  | Array containing the coordinates of the points       |
+      | =R(3)=       | input  | Array containing the x,y,z coordinates of the center |
+      | =n=          | input  | Number of computed gaussians                         |
+      | =A(n)=       | input  | Exponents of the Gaussians                           |
+      | =VGL(ldv,5)= | output | Value, gradients and Laplacian of the Gaussians      |
+      | =ldv=        | input  | Leading dimension of array =VGL=                     |
+    
+***** Requirements
+
+      - =context= is not 0
+      - =n= > 0
+      - =ldv= >= 5
+      - =A(i)= > 0 for all =i=
+      - =X= is allocated with at least $3 \times 8$ bytes
+      - =R= is allocated with at least $3 \times 8$ bytes
+      - =A= is allocated with at least $n \times 8$ bytes
+      - =VGL= is allocated with at least $n \times 5 \times 8$ bytes
+
+***** Header
+      #+BEGIN_SRC C :tangle qmckl.h
+qmckl_exit_code qmckl_ao_gaussian_vgl(const qmckl_context context,
+                const double *X, const double *R,
+                const int64_t *n, const int64_t *A,
+                const double *VGL,  const int64_t ldv);
+      #+END_SRC
+    
+***** Source
+      #+BEGIN_SRC f90 :tangle qmckl_ao.f90
+integer function qmckl_ao_gaussian_vgl_f(context, X, R, n, A, VGL, ldv) result(info)
+  implicit none
+  integer*8 , intent(in)  :: context
+  real*8    , intent(in)  :: X(3), R(3)
+  integer*8 , intent(in)  :: n
+  real*8    , intent(in)  :: A(n)
+  real*8    , intent(out) :: VGL(ldv,5)
+  integer*8 , intent(in)  :: ldv
+
+  integer*8         :: i,j
+  real*8            :: Y(3), r2, t, u, v
+  
+  info = 0
+  
+  if (context == 0_8) then
+     info = -1
+     return
+  endif
+  
+  if (n <= 0) then
+     info = -2
+     return
+  endif
+  
+  if (ldv < n) then
+     info = -3
+     return
+  endif
+  
+  
+  do i=1,3
+     Y(i) = X(i) - R(i)
+  end do
+  r2 = Y(1)*Y(1) + Y(2)*Y(2) + Y(3)*Y(3)
+  
+  do i=1,n
+     VGL(i,1) = dexp(-A(i) * r2)
+  end do
+
+  do i=1,n
+     VGL(i,5) = A(i) * VGL(i,1)
+  end do
+
+  t = -2.d0 * ( X(1) - R(1) )
+  u = -2.d0 * ( X(2) - R(2) )
+  v = -2.d0 * ( X(3) - R(3) )
+
+  do i=1,n
+     VGL(i,2) = t * VGL(i,5)
+     VGL(i,3) = u * VGL(i,5)
+     VGL(i,4) = v * VGL(i,5)
+  end do
+
+  t = 4.d0 * r2
+  do i=1,n
+     VGL(i,5) = (t * A(i) - 6.d0) *  VGL(i,5)
+  end do
+
+end function qmckl_ao_gaussian_vgl_f
+      #+END_SRC 
+
+***** C interface                                                  :noexport:
+      #+BEGIN_SRC f90 :tangle qmckl_ao.f90
+integer(c_int32_t) function qmckl_ao_gaussian_vgl(context, X, R, n, A, VGL, ldv) &
+     bind(C) result(info)
+  use, intrinsic :: iso_c_binding
+  implicit none
+  integer (c_int64_t) , intent(in) , value :: context
+  real    (c_double)  , intent(in)         :: X(3), R(3)
+  integer (c_int64_t) , intent(in) , value :: n
+  real    (c_double)  , intent(in)         :: A(n)
+  real    (c_double)  , intent(out)        :: VGL(ldv,5)
+  integer (c_int64_t) , intent(in) , value :: ldv
+
+  integer, external :: qmckl_ao_gaussian_vgl_f
+  info = qmckl_ao_gaussian_vgl_f(context, X, R, n, A, VGL, ldv) 
+end function qmckl_ao_gaussian_vgl
+      #+END_SRC
+
+      #+BEGIN_SRC f90 :tangle qmckl_f.f90
+  interface
+     integer(c_int32_t) function qmckl_ao_gaussian_vgl(context, X, R, n, A, VGL, ldv) &
+          bind(C) 
+       use, intrinsic :: iso_c_binding
+       integer (c_int64_t) , intent(in) , value :: context
+       integer (c_int64_t) , intent(in) , value :: ldv
+       integer (c_int64_t) , intent(in) , value :: n
+       real    (c_double)  , intent(in)         :: X(3), R(3), A(n)
+       real    (c_double)  , intent(out)        :: VGL(ldv,5)
+     end function qmckl_ao_gaussian_vgl
+  end interface
+      #+END_SRC
+***** Test                                                         :noexport:
+      #+BEGIN_SRC f90 :tangle test_qmckl_ao_f.f90
+integer(c_int32_t) function test_qmckl_ao_gaussian_vgl(context) bind(C)
+  use qmckl
+  implicit none
+
+  integer(c_int64_t), intent(in), value :: context
+  
+  integer*8                     :: n, ldv, j, i
+  double precision              :: X(3), R(3), Y(3), r2
+  double precision, allocatable :: VGL(:,:), A(:)
+  double precision              :: epsilon
+
+  epsilon = qmckl_context_get_epsilon(context)
+
+  X = (/ 1.1 , 2.2 ,  3.3 /)
+  R = (/ 0.1 , 1.2 , -2.3 /)
+  Y(:) = X(:) - R(:)
+  r2 = Y(1)**2 + Y(2)**2 + Y(3)**2
+
+  n = 10;
+  ldv = 100;
+
+  allocate (A(n), VGL(ldv,5))
+  do i=1,n
+     A(i) = 0.0013 * dble(ishft(1,i))
+  end do
+
+
+  test_qmckl_ao_gaussian_vgl = &
+       qmckl_ao_gaussian_vgl(context, X, R, n, A, VGL, ldv)
+  if (test_qmckl_ao_gaussian_vgl /= 0) return
+
+  test_qmckl_ao_gaussian_vgl = -1
+
+  do i=1,n
+     test_qmckl_ao_gaussian_vgl = -11
+     if (dabs(1.d0 - VGL(i,1) / (&
+          dexp(-A(i) * r2) &
+          )) > epsilon ) return
+     
+     test_qmckl_ao_gaussian_vgl = -12
+     if (dabs(1.d0 - VGL(i,2) / (&
+          -2.d0 * A(i) * Y(1) * dexp(-A(i) * r2) &
+          )) > epsilon ) return
+     
+     test_qmckl_ao_gaussian_vgl = -13
+     if (dabs(1.d0 - VGL(i,3) / (&
+          -2.d0 * A(i) * Y(2) * dexp(-A(i) * r2) &
+          )) > epsilon ) return
+     
+     test_qmckl_ao_gaussian_vgl = -14
+     if (dabs(1.d0 - VGL(i,4) / (&
+          -2.d0 * A(i) * Y(3) * dexp(-A(i) * r2) &
+          )) > epsilon ) return
+     
+     test_qmckl_ao_gaussian_vgl = -15
+     if (dabs(1.d0 - VGL(i,5) / (&
+          A(i) * (4.d0*r2*A(i) - 6.d0) * dexp(-A(i) * r2) &
+          )) > epsilon ) return
+  end do
+
+  test_qmckl_ao_gaussian_vgl = 0
+     
+  deallocate(VGL)
+end function test_qmckl_ao_gaussian_vgl
+      #+END_SRC
+
+      #+BEGIN_SRC C :tangle test_qmckl_ao.c
+int  test_qmckl_ao_gaussian_vgl(qmckl_context context);
+munit_assert_int(0, ==, test_qmckl_ao_gaussian_vgl(context));
+      #+END_SRC
+   
+    
+*** TODO Slater basis functions
+  
+*** End of files                                                   :noexport:
+  
+***** Test
+    #+BEGIN_SRC C :tangle test_qmckl_ao.c
+  if (qmckl_context_destroy(context) != QMCKL_SUCCESS)
+    return QMCKL_FAILURE;
+  return MUNIT_OK;
+}
+
+    #+END_SRC
+  
+
+  # -*- mode: org -*-
+  # vim: syntax=c

src/qmckl_context.org

@@ -1,193 +1,820 @@
-# -*- mode: org -*-
+** Context
 
-#+TITLE: Context
+   This file is written in C because it is more natural to express the
-
+   context in C than in Fortran.
-This file is written in C because it is more natural to express the context in
-C than in Fortran.
 
+   2 files are produced:
+   - a source file : =qmckl_context.c=
+   - a test   file : =test_qmckl_context.c=
 
+*** Headers                                                        :noexport:
     #+BEGIN_SRC C :tangle qmckl_context.c
-#include <stdlib.h> /* malloc */
+#include "qmckl.h"
-#include "qmckl_context.h"
     #+END_SRC
 
-* Context
+    #+BEGIN_SRC C :tangle test_qmckl_context.c
+#include "qmckl.h"
+#include "munit.h"
+MunitResult test_qmckl_context() {
+    #+END_SRC
+
+*** Context 
 
     The context variable is a handle for the state of the library, and
     is stored  in the  following data structure,  which can't  be seen
-  outside of the library.
+    outside  of  the library.  To  simplify  compatibility with  other
+    languages, the pointer to the internal data structure is converted
+    into a 64-bit signed integer, defined in the =qmckl_context= type.
+    A value of 0 for the context is equivalent to a =NULL= pointer.
 
+    # The following code block should be kept to insert comments into
+    # the qmckl.h file
 
-  #+BEGIN_SRC C :tangle qmckl_context.h
+    #+BEGIN_SRC C :comments org :tangle qmckl.h :export none
-#define QMCKL_DEFAULT_PRECISION 53
-#define QMCKL_DEFAULT_RANGE     2
-
-/* 64-bit integer */
-typedef long long int qmckl_context ;
     #+END_SRC
 
+**** Basis set data structure
+
+     Data  structure for  the info  related to  the atomic  orbitals
+     basis set.
+
+     #+BEGIN_SRC C :comments org :tangle qmckl_context.c
+typedef struct qmckl_ao_basis_struct {
+
+  int64_t   shell_num;      
+  int64_t   prim_num;      
+  int64_t * shell_center;  
+  int32_t * shell_ang_mom; 
+  double  * shell_factor;  
+  double  * exponent    ;  
+  double  * coefficient ;  
+  int64_t * shell_prim_num;
+  char      type;          
+
+} qmckl_ao_basis_struct;
+     #+END_SRC
+
+**** Source
+      
+     The tag is used internally to check if the memory domain pointed
+     by a pointer is a valid context.
+
+     #+BEGIN_SRC C :comments org :tangle qmckl_context.c
+typedef struct qmckl_context_struct {
+
+  struct qmckl_context_struct * prev;
+
+ /* Molecular system */
+  //  struct qmckl_nucleus_struct     * nucleus;  
+  //  struct qmckl_electron_struct    * electron;
+  struct qmckl_ao_basis_struct      * ao_basis;
+  //  struct qmckl_mo_struct          * mo;
+  //  struct qmckl_determinant_struct * det;
+
+ /* Numerical precision */
+  uint32_t tag;
+  int32_t precision;
+  int32_t range;
 
-  #+BEGIN_SRC C :tangle qmckl_context.c
-typedef struct qmckl_context_struct_ {
-  struct qmckl_context_struct_ * prev;
-  int precision;
-  int range;
 } qmckl_context_struct;
+
+#define VALID_TAG   0xBEEFFACE
+#define INVALID_TAG 0xDEADBEEF
      #+END_SRC
 
-** =qmckl_context_create=
+**** Test                                                          :noexport:
+     #+BEGIN_SRC C :tangle test_qmckl_context.c
+qmckl_context  context;
+qmckl_context  new_context;
+     #+END_SRC
 
-   To create a new context, use =qmckl_context_create()=. If the creation
-   failed, the function returns =0=. On success, a pointer to a context
-   is returned as a 64-bit integer.
 
-   #+BEGIN_SRC C :tangle qmckl_context.h
+**** =qmckl_context_check=
+
+     Checks if the  domain pointed by the pointer is  a valid context.
+     Returns  the input =qmckl_context= if  the context  is valid,  0
+     otherwise.
+
+     #+BEGIN_SRC C :comments org :tangle qmckl.h
+qmckl_context qmckl_context_check(const qmckl_context context) ;
+     #+END_SRC
+
+***** Source
+      #+BEGIN_SRC C :tangle qmckl_context.c
+qmckl_context qmckl_context_check(const qmckl_context context) {
+
+  if (context == (qmckl_context) 0) return (qmckl_context) 0;
+
+  const qmckl_context_struct * ctx = (qmckl_context_struct*) context;
+
+  if (ctx->tag != VALID_TAG) return (qmckl_context) 0;
+
+  return context;
+}
+      #+END_SRC
+
+**** =qmckl_context_create=
+
+     To create a new context, use =qmckl_context_create()=.
+     - On success, returns a pointer to a context using the =qmckl_context= type
+     - Returns 0 upon failure to allocate the internal data structure
+
+       #+BEGIN_SRC C :comments org :tangle qmckl.h
 qmckl_context qmckl_context_create();
        #+END_SRC
 
+***** Source
       #+BEGIN_SRC C :tangle qmckl_context.c
 qmckl_context qmckl_context_create() {
 
-  qmckl_context_struct*  context;
+  qmckl_context_struct* context =
-
+    (qmckl_context_struct*) qmckl_malloc ((qmckl_context) 0, sizeof(qmckl_context_struct));
-  context = (qmckl_context_struct*) malloc (sizeof(qmckl_context_struct));
   if (context == NULL) {
     return (qmckl_context) 0;
   }
 
   context->prev      = NULL;
+  context->ao_basis  = NULL;
   context->precision = QMCKL_DEFAULT_PRECISION;
   context->range     = QMCKL_DEFAULT_RANGE;
+  context->tag       = VALID_TAG;
 
   return (qmckl_context) context;
 }
       #+END_SRC
 
-** =qmckl_context_copy=
+***** Fortran interface
+      #+BEGIN_SRC f90 :tangle qmckl_f.f90
+  interface
+     integer (c_int64_t) function qmckl_context_create() bind(C)
+       use, intrinsic :: iso_c_binding
+     end function qmckl_context_create
+  end interface
+      #+END_SRC
 
-   #+BEGIN_SRC C :tangle qmckl_context.h
+***** Test                                                         :noexport:
+      #+BEGIN_SRC C :comments link :tangle test_qmckl_context.c
+context = qmckl_context_create();
+munit_assert_int64( context, !=, (qmckl_context) 0);
+munit_assert_int64( qmckl_context_check(context), ==,  context);
+      #+END_SRC
+
+**** =qmckl_context_copy=
+
+     This function makes a shallow copy of the current context.
+     - Copying the 0-valued context returns 0
+     - On success, returns a pointer to the new context using the =qmckl_context= type
+     - Returns 0 upon failure to allocate the internal data structure
+       for the new context
+
+       #+BEGIN_SRC C :comments org :tangle qmckl.h
 qmckl_context qmckl_context_copy(const qmckl_context context);
        #+END_SRC
 
+***** Source
       #+BEGIN_SRC C :tangle qmckl_context.c
 qmckl_context qmckl_context_copy(const qmckl_context context) {
 
-  qmckl_context_struct* old_context;
+  const qmckl_context checked_context = qmckl_context_check(context);
-  qmckl_context_struct* new_context;
 
-  new_context = (qmckl_context_struct*) malloc (sizeof(qmckl_context_struct));
+  if (checked_context == (qmckl_context) 0) {
+    return (qmckl_context) 0;
+  }
+
+  qmckl_context_struct* old_context = (qmckl_context_struct*) checked_context;
+
+  qmckl_context_struct* new_context = 
+    (qmckl_context_struct*) qmckl_malloc (context, sizeof(qmckl_context_struct));
   if (new_context == NULL) {
     return (qmckl_context) 0;
   }
 
-  old_context = (qmckl_context_struct*) context;
-
   new_context->prev      = old_context;
+  new_context->ao_basis  = old_context->ao_basis;
   new_context->precision = old_context->precision;
   new_context->range     = old_context->range;
+  new_context->tag       = VALID_TAG;
 
   return (qmckl_context) new_context;
 }
+
       #+END_SRC
 
-** =qmckl_context_destroy=
+***** Fortran interface
-
+      #+BEGIN_SRC f90 :tangle qmckl_f.f90
-   To delete a new context, use =qmckl_context_destroy()=. If the deletion
+  interface
-   failed, the function returns =0=. On success, the function returns =1=
+     integer (c_int64_t) function qmckl_context_copy(context) bind(C)
-   implying that the context has been freed.
+       use, intrinsic :: iso_c_binding
-
+       integer (c_int64_t), intent(in), value :: context
-   #+BEGIN_SRC C :tangle qmckl_context.h
+     end function qmckl_context_copy
-int qmckl_context_destroy(qmckl_context context);
+  end interface
       #+END_SRC
 
+***** Test                                                         :noexport:
+      #+BEGIN_SRC C :comments link :tangle test_qmckl_context.c
+new_context = qmckl_context_copy(context);
+munit_assert_int64(new_context, !=, (qmckl_context) 0);
+munit_assert_int64(new_context, !=, context);
+munit_assert_int64(qmckl_context_check(new_context), ==, new_context);
+      #+END_SRC
+
+**** =qmckl_context_previous=
+
+     Returns the previous context
+     - On success, returns the ancestor of the current context
+     - Returns 0 for the initial context
+     - Returns 0 for the 0-valued context
+
+       #+BEGIN_SRC C :comments org :tangle qmckl.h
+qmckl_context qmckl_context_previous(const qmckl_context context);
+       #+END_SRC
+
+***** Source
       #+BEGIN_SRC C :tangle qmckl_context.c
-int qmckl_context_destroy(qmckl_context context) {
+qmckl_context qmckl_context_previous(const qmckl_context context) {
 
-  qmckl_context_struct* ctx;
+  const qmckl_context checked_context = qmckl_context_check(context);
+  if (checked_context == (qmckl_context) 0) {
+    return (qmckl_context) 0;
+  }
 
-  ctx = (qmckl_context_struct*) context;
+  const qmckl_context_struct* ctx = (qmckl_context_struct*) checked_context;
+  return qmckl_context_check((qmckl_context) ctx->prev);
+}
+      #+END_SRC
 
-  if (ctx == NULL) {
+***** Fortran interface
+      #+BEGIN_SRC f90 :tangle qmckl_f.f90
+  interface
+     integer (c_int64_t) function qmckl_context_previous(context) bind(C)
+       use, intrinsic :: iso_c_binding
+       integer (c_int64_t), intent(in), value :: context
+     end function qmckl_context_previous
+  end interface
+      #+END_SRC
+
+***** Test                                                         :noexport:
+      #+BEGIN_SRC C :comments link :tangle test_qmckl_context.c
+munit_assert_int64(qmckl_context_previous(new_context), !=, (qmckl_context) 0);
+munit_assert_int64(qmckl_context_previous(new_context), ==, context);
+munit_assert_int64(qmckl_context_previous(context), ==, (qmckl_context) 0);
+munit_assert_int64(qmckl_context_previous((qmckl_context) 0), ==, (qmckl_context) 0);
+      #+END_SRC
+
+**** =qmckl_context_destroy=
+
+     Destroys the current context, leaving the ancestors untouched.
+     - Succeeds if the current context is properly destroyed
+     - Fails otherwise
+     - Fails if the 0-valued context is given in argument
+     - Fails if the the pointer is not a valid context
+
+     #+BEGIN_SRC C :comments org :tangle qmckl.h
+qmckl_exit_code qmckl_context_destroy(qmckl_context context);
+     #+END_SRC
+
+***** Source
+      #+BEGIN_SRC C :tangle qmckl_context.c
+qmckl_exit_code qmckl_context_destroy(const qmckl_context context) {
+
+  const qmckl_context checked_context = qmckl_context_check(context);
+  if (checked_context == (qmckl_context) 0) return QMCKL_FAILURE;
+
+  qmckl_context_struct* ctx = (qmckl_context_struct*) context;
+  if (ctx == NULL) return QMCKL_FAILURE;
+
+  ctx->tag = INVALID_TAG;
+  qmckl_free(ctx);
+  return QMCKL_SUCCESS;
+}
+      #+END_SRC
+
+***** Fortran interface
+      #+BEGIN_SRC f90 :tangle qmckl_f.f90
+  interface
+     integer (c_int32_t) function qmckl_context_destroy(context) bind(C)
+       use, intrinsic :: iso_c_binding
+       integer (c_int64_t), intent(in), value :: context
+     end function qmckl_context_destroy
+  end interface
+      #+END_SRC
+
+***** Test                                                         :noexport:
+      #+BEGIN_SRC C :tangle test_qmckl_context.c
+munit_assert_int64(qmckl_context_check(new_context), ==, new_context);
+munit_assert_int64(new_context, !=, (qmckl_context) 0);
+munit_assert_int32(qmckl_context_destroy(new_context), ==, QMCKL_SUCCESS);
+munit_assert_int64(qmckl_context_check(new_context), !=, new_context);
+munit_assert_int64(qmckl_context_check(new_context), ==, (qmckl_context) 0);
+munit_assert_int64(qmckl_context_destroy((qmckl_context) 0), ==, QMCKL_FAILURE);
+      #+END_SRC
+
+*** Basis set
+
+     For H_2 with the following basis set, 
+
+     #+BEGIN_EXAMPLE
+HYDROGEN
+S   5
+1         3.387000E+01           6.068000E-03
+2         5.095000E+00           4.530800E-02
+3         1.159000E+00           2.028220E-01                                                 
+4         3.258000E-01           5.039030E-01
+5         1.027000E-01           3.834210E-01
+S   1
+1         3.258000E-01           1.000000E+00
+S   1
+1         1.027000E-01           1.000000E+00
+P   1
+1         1.407000E+00           1.000000E+00
+P   1
+1         3.880000E-01           1.000000E+00
+D   1
+1         1.057000E+00           1.0000000
+     #+END_EXAMPLE
+
+     we have:
+
+     #+BEGIN_EXAMPLE
+type = 'G'
+shell_num = 12
+prim_num = 20
+SHELL_CENTER = [1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2]
+SHELL_ANG_MOM = ['S', 'S', 'S', 'P', 'P', 'D', 'S', 'S', 'S', 'P', 'P', 'D']
+SHELL_PRIM_NUM = [5, 1, 1, 1, 1, 1, 5, 1, 1, 1, 1, 1]
+prim_index = [1, 6, 7, 8, 9, 10, 11, 16, 17, 18, 19, 20]
+EXPONENT = [ 33.87, 5.095, 1.159, 0.3258, 0.1027, 0.3258, 0.1027,
+             1.407, 0.388, 1.057, 33.87, 5.095, 1.159, 0.3258, 0.1027,
+             0.3258, 0.1027, 1.407, 0.388, 1.057]
+COEFFICIENT = [ 0.006068, 0.045308, 0.202822, 0.503903, 0.383421,
+                1.0, 1.0, 1.0, 1.0, 1.0, 0.006068, 0.045308, 0.202822,
+                0.503903, 0.383421, 1.0, 1.0, 1.0, 1.0, 1.0]
+     #+END_EXAMPLE
+
+**** =qmckl_context_update_ao_basis=
+
+     Updates the data describing the AO basis set into the context.
+
+     | =type=                        | Gaussian or Slater                                                   |
+     | =shell_num=                   | Number of shells                                                     |
+     | =prim_num=                    | Total number of primitives                                           |
+     | =SHELL_CENTER(shell_num)=     | Id of the nucleus on which the shell is centered                     |
+     | =SHELL_ANG_MOM(shell_num)=    | Id of the nucleus on which the shell is centered                     |
+     | =SHELL_FACTOR(shell_num)=     | Normalization factor for the shell                                   |
+     | =SHELL_PRIM_NUM(shell_num)=   | Number of primitives in the shell                                    |
+     | =SHELL_PRIM_INDEX(shell_num)= | Address of the first primitive of the shelll in the =EXPONENT= array |
+     | =EXPONENT(prim_num)=          | Array of exponents                                                   |
+     | =COEFFICIENT(prim_num)=       | Array of coefficients                                                |
+
+     #+BEGIN_SRC C :comments org :tangle qmckl.h
+qmckl_exit_code
+qmckl_context_update_ao_basis(qmckl_context   context     , const char      type,
+                              const int64_t   shell_num   , const int64_t   prim_num, 
+                              const int64_t * SHELL_CENTER, const int32_t * SHELL_ANG_MOM,
+                              const double  * SHELL_FACTOR, const int64_t * SHELL_PRIM_NUM,
+                              const int64_t * SHELL_PRIM_INDEX,
+                              const double  * EXPONENT    , const double  * COEFFICIENT);
+     #+END_SRC
+
+***** Source
+      #+BEGIN_SRC C :tangle qmckl_context.c
+qmckl_exit_code
+qmckl_context_update_ao_basis(qmckl_context   context     , const char      type,
+                              const int64_t   shell_num   , const int64_t   prim_num, 
+                              const int64_t * SHELL_CENTER, const int32_t * SHELL_ANG_MOM,
+                              const double  * SHELL_FACTOR, const int64_t * SHELL_PRIM_NUM,
+                              const int64_t * SHELL_PRIM_INDEX,
+                              const double  * EXPONENT    , const double  * COEFFICIENT)
+{
+
+  int64_t i;
+
+ /* Check input */
+
+  if (type != 'G' && type != 'S') return QMCKL_FAILURE;
+  if (shell_num <= 0) return QMCKL_FAILURE;
+  if (prim_num  <= 0) return QMCKL_FAILURE;
+  if (prim_num  < shell_num) return QMCKL_FAILURE;
+  
+  for (i=0 ; i<shell_num ; i++) {
+    if (SHELL_CENTER[i] <= 0) return QMCKL_FAILURE;
+    if (SHELL_PRIM_NUM[i] <= 0) return QMCKL_FAILURE;
+    if (SHELL_ANG_MOM[i] < 0) return QMCKL_FAILURE;
+    if (SHELL_PRIM_INDEX[i] < 0) return QMCKL_FAILURE;
+  }
+  
+  for (i=0 ; i<prim_num ; i++) {
+    if (EXPONENT[i] <= 0) return QMCKL_FAILURE;
+  }
+  
+  qmckl_context_struct* ctx = (qmckl_context_struct*) context;
+  if (ctx == NULL) return QMCKL_FAILURE;
+  
+  qmckl_ao_basis_struct* basis = (qmckl_ao_basis_struct*) malloc (sizeof(qmckl_ao_basis_struct));
+  if (basis == NULL) return QMCKL_FAILURE;
+  
+
+ /* Memory allocations */
+
+  basis->shell_center  = (int64_t*) malloc (shell_num * sizeof(int64_t));
+  if (basis->shell_center == NULL) {
+    free(basis);
+    return QMCKL_FAILURE;
+  }
+  
+  basis->shell_ang_mom = (int32_t*) malloc (shell_num * sizeof(int32_t));
+  if (basis->shell_ang_mom == NULL) {
+    free(basis->shell_center);
+    free(basis);
+    return QMCKL_FAILURE;
+  }
+  
+  basis->shell_prim_num= (int64_t*) malloc (shell_num * sizeof(int64_t));
+  if (basis->shell_prim_num == NULL) {
+    free(basis->shell_ang_mom);
+    free(basis->shell_center);
+    free(basis);
+    return QMCKL_FAILURE;
+  }
+  
+  basis->shell_factor  = (double *) malloc (shell_num * sizeof(double ));
+  if (basis->shell_factor == NULL) {
+    free(basis->shell_prim_num);
+    free(basis->shell_ang_mom);
+    free(basis->shell_center);
+    free(basis);
+    return QMCKL_FAILURE;
+  }
+
+  basis->exponent      = (double *) malloc (prim_num  * sizeof(double ));
+  if (basis->exponent == NULL) {
+    free(basis->shell_factor);
+    free(basis->shell_prim_num);
+    free(basis->shell_ang_mom);
+    free(basis->shell_center);
+    free(basis);
+    return QMCKL_FAILURE;
+  }
+
+  basis->coefficient   = (double *) malloc (prim_num  * sizeof(double ));
+  if (basis->coefficient == NULL) {
+    free(basis->exponent);
+    free(basis->shell_factor);
+    free(basis->shell_prim_num);
+    free(basis->shell_ang_mom);
+    free(basis->shell_center);
+    free(basis);
+    return QMCKL_FAILURE;
+  }
+
+  
+ /* Assign data */
+
+  basis->type      = type;
+  basis->shell_num = shell_num;
+  basis->prim_num  = prim_num;      
+
+  for (i=0 ; i<shell_num ; i++) {
+    basis->shell_center  [i] = SHELL_CENTER  [i];
+    basis->shell_ang_mom [i] = SHELL_ANG_MOM [i];
+    basis->shell_prim_num[i] = SHELL_PRIM_NUM[i];
+    basis->shell_factor  [i] = SHELL_FACTOR  [i];
+  }
+
+  for (i=0 ; i<prim_num ; i++) {
+    basis->exponent   [i] = EXPONENT[i];
+    basis->coefficient[i] = COEFFICIENT[i];
+  }
+
+  ctx->ao_basis = basis;
+  return QMCKL_SUCCESS;
+}
+      #+END_SRC
+      
+***** Fortran interface
+      #+BEGIN_SRC f90 :tangle qmckl_f.f90
+  interface
+     integer (c_int32_t) function qmckl_context_update_ao_basis(context, &
+          typ, shell_num, prim_num, SHELL_CENTER, SHELL_ANG_MOM, SHELL_FACTOR, &
+          SHELL_PRIM_NUM, SHELL_PRIM_INDEX, EXPONENT, COEFFICIENT) bind(C)
+       use, intrinsic :: iso_c_binding
+       integer (c_int64_t), intent(in), value :: context
+       character(c_char)  , intent(in), value :: typ
+       integer (c_int64_t), intent(in), value :: shell_num
+       integer (c_int64_t), intent(in), value :: prim_num
+       integer (c_int64_t), intent(in)        :: SHELL_CENTER(shell_num)
+       integer (c_int32_t), intent(in)        :: SHELL_ANG_MOM(shell_num)
+       double precision   , intent(in)        :: SHELL_FACTOR(shell_num)
+       integer (c_int64_t), intent(in)        :: SHELL_PRIM_NUM(shell_num)
+       integer (c_int64_t), intent(in)        :: SHELL_PRIM_INDEX(shell_num)
+       double precision   , intent(in)        :: EXPONENT(prim_num)
+       double precision   , intent(in)        :: COEFFICIENT(prim_num)
+     end function qmckl_context_update_ao_basis
+  end interface
+      #+END_SRC
+
+***** TODO Test
+
+**** =qmckl_context_set_ao_basis=
+
+     Sets the data describing the AO basis set into the context.
+
+     | =type=                        | Gaussian or Slater                                                   |
+     | =shell_num=                   | Number of shells                                                     |
+     | =prim_num=                    | Total number of primitives                                           |
+     | =SHELL_CENTER(shell_num)=     | Id of the nucleus on which the shell is centered                     |
+     | =SHELL_ANG_MOM(shell_num)=    | Id of the nucleus on which the shell is centered                     |
+     | =SHELL_FACTOR(shell_num)=     | Normalization factor for the shell                                   |
+     | =SHELL_PRIM_NUM(shell_num)=   | Number of primitives in the shell                                    |
+     | =SHELL_PRIM_INDEX(shell_num)= | Address of the first primitive of the shelll in the =EXPONENT= array |
+     | =EXPONENT(prim_num)=          | Array of exponents                                                   |
+     | =COEFFICIENT(prim_num)=       | Array of coefficients                                                |
+
+     #+BEGIN_SRC C :comments org :tangle qmckl.h
+qmckl_context
+qmckl_context_set_ao_basis(const qmckl_context context     , const char      type,
+                           const int64_t       shell_num   , const int64_t   prim_num, 
+                           const int64_t *     SHELL_CENTER, const int32_t * SHELL_ANG_MOM,
+                           const double  *     SHELL_FACTOR, const int64_t * SHELL_PRIM_NUM,
+                           const int64_t *     SHELL_PRIM_INDEX,
+                           const double  *     EXPONENT    , const double  * COEFFICIENT);
+     #+END_SRC
+
+***** Source
+      #+BEGIN_SRC C :tangle qmckl_context.c
+qmckl_context
+qmckl_context_set_ao_basis(const qmckl_context context     , const char      type,
+                           const int64_t       shell_num   , const int64_t   prim_num, 
+                           const int64_t *     SHELL_CENTER, const int32_t * SHELL_ANG_MOM,
+                           const double  *     SHELL_FACTOR, const int64_t * SHELL_PRIM_NUM,
+                           const int64_t *     SHELL_PRIM_INDEX,
+                           const double  *     EXPONENT    , const double  * COEFFICIENT)
+{
+
+  qmckl_context new_context = qmckl_context_copy(context);
+  if (new_context == 0) return 0;
+
+  if (qmckl_context_update_ao_basis(context, type, shell_num, prim_num, 
+                                    SHELL_CENTER, SHELL_ANG_MOM, SHELL_FACTOR, 
+                                    SHELL_PRIM_NUM, SHELL_PRIM_INDEX, EXPONENT,
+                                    COEFFICIENT
+                                    ) == QMCKL_FAILURE)
     return 0;
-  }
 
-  free(ctx);
+  return new_context;
-  return 1;
 }
       #+END_SRC
       
-* Precision
+***** Fortran interface
-
+      #+BEGIN_SRC f90 :tangle qmckl_f.f90
-  The following functions set and get the expected required precision
+  interface
-  and range. =precision= should be an integer between 2 and 53, and
+     integer (c_int64_t) function qmckl_context_set_ao_basis(context, &
- =range= should be an integer between 2 and 11.
+          typ, shell_num, prim_num, SHELL_CENTER, SHELL_ANG_MOM, SHELL_FACTOR, &
-  The setter functions functions return a new context as a 64-bit integer.
+          SHELL_PRIM_NUM, SHELL_PRIM_INDEX, EXPONENT, COEFFICIENT) bind(C)
-  The getter functions return the value, as a 32-bit integer.
+       use, intrinsic :: iso_c_binding
-
+       integer (c_int64_t), intent(in), value :: context
-** =qmckl_context_set_precision=
+       character(c_char)  , intent(in), value :: typ
-
+       integer (c_int64_t), intent(in), value :: shell_num
-   #+BEGIN_SRC C :tangle qmckl_context.h
+       integer (c_int64_t), intent(in), value :: prim_num
-qmckl_context qmckl_context_set_precision(const qmckl_context context, int precision);
+       integer (c_int64_t), intent(in)        :: SHELL_CENTER(shell_num)
+       integer (c_int32_t), intent(in)        :: SHELL_ANG_MOM(shell_num)
+       double precision   , intent(in)        :: SHELL_FACTOR(shell_num)
+       integer (c_int64_t), intent(in)        :: SHELL_PRIM_NUM(shell_num)
+       integer (c_int64_t), intent(in)        :: SHELL_PRIM_INDEX(shell_num)
+       double precision   , intent(in)        :: EXPONENT(prim_num)
+       double precision   , intent(in)        :: COEFFICIENT(prim_num)
+     end function qmckl_context_set_ao_basis
+  end interface
       #+END_SRC
 
+***** TODO Test
+
+*** Precision
+
+    The  following  functions  set   and  get  the  expected  required
+    precision and  range. =precision= should  be an integer  between 2
+    and 53, and =range= should be an integer between 2 and 11.
+
+    The setter  functions functions return  a new context as  a 64-bit
+    integer.   The getter  functions  return the  value,  as a  32-bit
+    integer.    The  update   functions   return =QMCKL_SUCCESS=   or
+ =QMCKL_FAILURE=.
+
+**** =qmckl_context_update_precision=
+     Modifies the parameter for the numerical precision in a given context.
+     #+BEGIN_SRC C :comments org :tangle qmckl.h
+qmckl_exit_code qmckl_context_update_precision(const qmckl_context context, const int precision);
+     #+END_SRC
+
+***** Source
       #+BEGIN_SRC C :tangle qmckl_context.c
-qmckl_context qmckl_context_set_precision(const qmckl_context context, int precision) {
+qmckl_exit_code qmckl_context_update_precision(const qmckl_context context, const int precision) {
-  qmckl_context_struct* ctx;
 
-  if (precision <  2) return (qmckl_context) 0;
+  if (precision <  2) return QMCKL_FAILURE;
-  if (precision > 53) return (qmckl_context) 0;
+  if (precision > 53) return QMCKL_FAILURE;
+
+  qmckl_context_struct* ctx = (qmckl_context_struct*) context;
+  if (ctx == NULL) return QMCKL_FAILURE;
 
-  ctx = (qmckl_context_struct*) qmckl_context_copy(context);
   ctx->precision = precision;
-  return (qmckl_context) ctx;
+  return QMCKL_SUCCESS;
 }
       #+END_SRC
 
-** =qmckl_context_set_range=
+***** Fortran interface
-   #+BEGIN_SRC C :tangle qmckl_context.h
+      #+BEGIN_SRC f90 :tangle qmckl_f.f90
-qmckl_context qmckl_context_set_range(const qmckl_context context, int range);
+  interface
+     integer (c_int32_t) function qmckl_context_update_precision(context, precision) bind(C)
+       use, intrinsic :: iso_c_binding
+       integer (c_int64_t), intent(in), value :: context
+       integer (c_int32_t), intent(in), value :: precision
+     end function qmckl_context_update_precision
+  end interface
       #+END_SRC
 
+***** TODO Tests                                                   :noexport:
+**** =qmckl_context_update_range=
+     Modifies the parameter for the numerical range in a given context.
+     #+BEGIN_SRC C :comments org :tangle qmckl.h
+qmckl_exit_code qmckl_context_update_range(const qmckl_context context, const int range);
+      #+END_SRC
+
+***** Source
       #+BEGIN_SRC C :tangle qmckl_context.c
-qmckl_context qmckl_context_set_range(const qmckl_context context, int range) {
+qmckl_exit_code qmckl_context_update_range(const qmckl_context context, const int range) {
-  qmckl_context_struct* ctx;
 
-  if (range <  2) return (qmckl_context) 0;
+  if (range <  2) return QMCKL_FAILURE;
-  if (range > 11) return (qmckl_context) 0;
+  if (range > 11) return QMCKL_FAILURE;
+
+  qmckl_context_struct* ctx = (qmckl_context_struct*) context;
+  if (ctx == NULL) return QMCKL_FAILURE;
 
-  ctx = (qmckl_context_struct*) qmckl_context_copy(context);
   ctx->range = range;
-  return (qmckl_context) ctx;
+  return QMCKL_SUCCESS;
 }
       #+END_SRC
 
-
+***** Fortran interface
-
+      #+BEGIN_SRC f90 :tangle qmckl_f.f90
-** =qmckl_context_get_precision=
+  interface
-
+     integer (c_int32_t) function qmckl_context_update_range(context, range) bind(C)
-   #+BEGIN_SRC C :tangle qmckl_context.h
+       use, intrinsic :: iso_c_binding
-int qmckl_context_get_precision(const qmckl_context context);
+       integer (c_int64_t), intent(in), value :: context
+       integer (c_int32_t), intent(in), value :: range
+     end function qmckl_context_update_range
+  end interface
       #+END_SRC
 
+***** TODO Tests                                                   :noexport:
+**** =qmckl_context_set_precision=
+     Returns a copy of the context with a different precision parameter.
+     #+BEGIN_SRC C :comments org :tangle qmckl.h
+qmckl_context qmckl_context_set_precision(const qmckl_context context, const int precision);
+     #+END_SRC
+
+***** Source
+      #+BEGIN_SRC C :tangle qmckl_context.c
+qmckl_context qmckl_context_set_precision(const qmckl_context context, const int precision) {
+  qmckl_context new_context = qmckl_context_copy(context);
+  if (new_context == 0) return 0;
+
+  if (qmckl_context_update_precision(context, precision) == QMCKL_FAILURE) return 0;
+
+  return new_context;
+}
+      #+END_SRC
+
+***** Fortran interface
+      #+BEGIN_SRC f90 :tangle qmckl_f.f90
+  interface
+     integer (c_int64_t) function qmckl_context_set_precision(context, precision) bind(C)
+       use, intrinsic :: iso_c_binding
+       integer (c_int64_t), intent(in), value :: context
+       integer (c_int32_t), intent(in), value :: precision
+     end function qmckl_context_set_precision
+  end interface
+      #+END_SRC
+
+***** TODO Tests                                                   :noexport:
+**** =qmckl_context_set_range=
+     Returns a copy of the context with a different precision parameter.
+     #+BEGIN_SRC C :comments org :tangle qmckl.h
+qmckl_context qmckl_context_set_range(const qmckl_context context, const int range);
+     #+END_SRC
+
+***** Source
+      #+BEGIN_SRC C :tangle qmckl_context.c
+qmckl_context qmckl_context_set_range(const qmckl_context context, const int range) {
+  qmckl_context new_context = qmckl_context_copy(context);
+  if (new_context == 0) return 0;
+
+  if (qmckl_context_update_range(context, range) == QMCKL_FAILURE) return 0;
+
+  return new_context;
+}
+      #+END_SRC
+
+***** Fortran interface
+      #+BEGIN_SRC f90 :tangle qmckl_f.f90
+  interface
+     integer (c_int64_t) function qmckl_context_set_range(context, range) bind(C)
+       use, intrinsic :: iso_c_binding
+       integer (c_int64_t), intent(in), value :: context
+       integer (c_int32_t), intent(in), value :: range
+     end function qmckl_context_set_range
+  end interface
+      #+END_SRC
+
+***** TODO Tests                                                   :noexport:
+
+**** =qmckl_context_get_precision=
+     Returns the value of the numerical precision in the context
+     #+BEGIN_SRC C :comments org :tangle qmckl.h
+int32_t qmckl_context_get_precision(const qmckl_context context);
+     #+END_SRC
+
+***** Source
       #+BEGIN_SRC C :tangle qmckl_context.c
 int qmckl_context_get_precision(const qmckl_context context) {
-  qmckl_context_struct* ctx;
+  const qmckl_context_struct* ctx = (qmckl_context_struct*) context;
-  ctx = (qmckl_context_struct*) context;
   return ctx->precision;
 }
       #+END_SRC
 
-** =qmckl_context_get_range=
+***** Fortran interface
-
+      #+BEGIN_SRC f90 :tangle qmckl_f.f90
-   #+BEGIN_SRC C :tangle qmckl_context.h
+  interface
-int qmckl_context_get_range(const qmckl_context context);
+     integer (c_int32_t) function qmckl_context_get_precision(context) bind(C)
+       use, intrinsic :: iso_c_binding
+       integer (c_int64_t), intent(in), value :: context
+     end function qmckl_context_get_precision
+  end interface
       #+END_SRC
 
+***** TODO Tests                                                   :noexport:
+**** =qmckl_context_get_range=
+     Returns the value of the numerical range in the context
+     #+BEGIN_SRC C :comments org :tangle qmckl.h
+int32_t qmckl_context_get_range(const qmckl_context context);
+     #+END_SRC
+
+***** Source
       #+BEGIN_SRC C :tangle qmckl_context.c
 int qmckl_context_get_range(const qmckl_context context) {
-  qmckl_context_struct* ctx;
+  const qmckl_context_struct* ctx = (qmckl_context_struct*) context;
-  ctx = (qmckl_context_struct*) context;
   return ctx->range;
 }
       #+END_SRC
 
+***** Fortran interface
+      #+BEGIN_SRC f90 :tangle qmckl_f.f90
+  interface
+     integer (c_int32_t) function qmckl_context_get_range(context) bind(C)
+       use, intrinsic :: iso_c_binding
+       integer (c_int64_t), intent(in), value :: context
+     end function qmckl_context_get_range
+  end interface
+      #+END_SRC
+
+***** TODO Tests                                                   :noexport:
+
+**** =qmckl_context_get_epsilon=
+     Returns $\epsilon = 2^{1-n}$ where =n= is the precision
+     #+BEGIN_SRC C :comments org :tangle qmckl.h
+double qmckl_context_get_epsilon(const qmckl_context context);
+     #+END_SRC
+
+***** Source
+      #+BEGIN_SRC C :tangle qmckl_context.c
+double qmckl_context_get_epsilon(const qmckl_context context) {
+  const qmckl_context_struct* ctx = (qmckl_context_struct*) context;
+  return pow(2.0,(double) 1-ctx->precision);
+}
+      #+END_SRC
+
+***** Fortran interface
+      #+BEGIN_SRC f90 :tangle qmckl_f.f90
+  interface
+     real (c_double) function qmckl_context_get_epsilon(context) bind(C)
+       use, intrinsic :: iso_c_binding
+       integer (c_int64_t), intent(in), value :: context
+     end function qmckl_context_get_epsilon
+  end interface
+      #+END_SRC
+
+***** TODO Tests                                                   :noexport:
+
+      
+
+*** End of files                                                   :noexport:
+
+***** Test
+      #+BEGIN_SRC C :comments link :tangle test_qmckl_context.c
+return MUNIT_OK;
+}
+      #+END_SRC
+
+  
+
+ # -*- mode: org -*-
+ # vim: syntax=c

src/qmckl_distance.org

@@ -0,0 +1,356 @@
+** Computation of distances
+
+ Function for the computation of distances between particles.
+
+ 3 files are produced:
+ - a source file : =qmckl_distance.f90=
+ - a C test file : =test_qmckl_distance.c=
+ - a Fortran test file : =test_qmckl_distance_f.f90=
+
+**** Headers                                                       :noexport:
+     #+BEGIN_SRC C :comments link :tangle test_qmckl_distance.c
+#include "qmckl.h"
+#include "munit.h"
+MunitResult test_qmckl_distance() {
+  qmckl_context context;
+  context = qmckl_context_create();
+
+     #+END_SRC
+
+
+*** Squared distance
+
+**** =qmckl_distance_sq=
+
+     Computes the matrix of the squared distances between all pairs of
+     points in two sets, one point within each set:
+     \[
+     C_{ij} = \sum_{k=1}^3 (A_{k,i}-B_{k,j})^2
+     \]
+
+***** Arguments
+
+      | =context=  | input  | Global state                                 |
+      | =transa=   | input  | Array =A= is =N=: Normal, =T=: Transposed   |
+      | =transb=   | input  | Array =B= is =N=: Normal, =T=: Transposed   |
+      | =m=        | input  | Number of points in the first set            |
+      | =n=        | input  | Number of points in the second set           |
+      | =A(lda,3)= | input  | Array containing the $m \times 3$ matrix $A$ |
+      | =lda=      | input  | Leading dimension of array =A=               |
+      | =B(ldb,3)= | input  | Array containing the $n \times 3$ matrix $B$ |
+      | =ldb=      | input  | Leading dimension of array =B=               |
+      | =C(ldc,n)= | output | Array containing the $m \times n$ matrix $C$ |
+      | =ldc=      | input  | Leading dimension of array =C=               |
+
+***** Requirements
+
+      - =context= is not 0
+      - =m= > 0
+      - =n= > 0
+      - =lda= >= 3 if =transa= is =N=
+      - =lda= >= m if =transa= is =T=
+      - =ldb= >= 3 if =transb= is =N=
+      - =ldb= >= n if =transb= is =T=
+      - =ldc= >= m if =transa= is =
+      - =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 $m \times n \times 8$ bytes
+
+***** Performance 
+
+      This function might be more efficient when =A= and =B= are
+      transposed.
+
+    #+BEGIN_SRC C :comments org :tangle qmckl.h
+qmckl_exit_code qmckl_distance_sq(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,
+                                  const double *C, const int64_t ldc);
+      #+END_SRC
+
+***** Source
+      #+BEGIN_SRC f90 :tangle qmckl_distance.f90
+integer function qmckl_distance_sq_f(context, transa, transb, m, n, A, LDA, B, LDB, C, LDC) result(info)
+  implicit none
+  integer*8  , 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*8 :: i,j
+  real*8    :: x, y, z
+  integer   :: transab
+
+  info = 0
+
+  if (context == 0_8) then
+     info = -1
+     return
+  endif
+
+  if (m <= 0_8) then
+     info = -2
+     return
+  endif
+
+  if (n <= 0_8) then
+     info = -3
+     return
+  endif
+
+  if (transa == 'N' .or. transa == 'n') then
+     transab = 0
+  else if (transa == 'T' .or. transa == 't') then
+     transab = 1
+  else
+     transab = -100
+  endif
+
+  if (transb == 'N' .or. transb == 'n') then
+     continue
+  else if (transa == 'T' .or. transa == 't') then
+     transab = transab + 2
+  else
+     transab = -100
+  endif
+
+  if (transab < 0) then
+     info = -4
+     return 
+  endif
+
+  if (iand(transab,1) == 0 .and. LDA < 3) then
+     info = -5
+     return
+  endif
+
+  if (iand(transab,1) == 1 .and. LDA < m) then
+     info = -6
+     return
+  endif
+
+  if (iand(transab,2) == 0 .and. LDA < 3) then
+     info = -6
+     return
+  endif
+
+  if (iand(transab,2) == 2 .and. LDA < m) then
+     info = -7
+     return
+  endif
+
+
+  select case (transab)
+     
+  case(0)
+
+     do j=1,n
+        do i=1,m
+           x = A(1,i) - B(1,j)
+           y = A(2,i) - B(2,j)
+           z = A(3,i) - B(3,j)
+           C(i,j) = x*x + y*y + z*z
+        end do
+     end do
+
+  case(1)
+
+     do j=1,n
+        do i=1,m
+           x = A(i,1) - B(1,j)
+           y = A(i,2) - B(2,j)
+           z = A(i,3) - B(3,j)
+           C(i,j) = x*x + y*y + z*z
+        end do
+     end do
+
+  case(2)
+
+     do j=1,n
+        do i=1,m
+           x = A(1,i) - B(j,1)
+           y = A(2,i) - B(j,2)
+           z = A(3,i) - B(j,3)
+           C(i,j) = x*x + y*y + z*z
+        end do
+     end do
+
+  case(3)
+
+     do j=1,n
+        do i=1,m
+           x = A(i,1) - B(j,1)
+           y = A(i,2) - B(j,2)
+           z = A(i,3) - B(j,3)
+           C(i,j) = x*x + y*y + z*z
+        end do
+     end do
+
+  end select
+  
+end function qmckl_distance_sq_f
+      #+END_SRC
+
+***** C interface                                                  :noexport:
+      #+BEGIN_SRC f90 :tangle qmckl_distance.f90
+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 (c_char)  , intent(in) , value :: transa, transb
+  integer (c_int64_t) , intent(in) , value :: m, n
+  integer (c_int64_t) , intent(in) , value :: lda
+  real    (c_double)  , intent(in)         :: A(lda,3)
+  integer (c_int64_t) , intent(in) , value :: ldb
+  real    (c_double)  , intent(in)         :: B(ldb,3)
+  integer (c_int64_t) , intent(in) , value :: ldc
+  real    (c_double)  , intent(out)        :: C(ldc,n)
+
+  integer, 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
+
+      #+BEGIN_SRC f90 :tangle qmckl_f.f90
+  interface
+     integer(c_int32_t) function qmckl_distance_sq(context, transa, transb, m, n, A, LDA, B, LDB, C, LDC) &
+          bind(C) 
+       use, intrinsic :: iso_c_binding
+       implicit none
+       integer (c_int64_t) , intent(in) , value :: context
+       character (c_char)  , intent(in) , value :: transa, transb
+       integer (c_int64_t) , intent(in) , value :: m, 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,3)
+       real    (c_double)  , intent(in)         :: B(ldb,3)
+       real    (c_double)  , intent(out)        :: C(ldc,n)
+     end function qmckl_distance_sq
+  end interface
+      #+END_SRC
+
+***** Test                                                         :noexport:
+      #+BEGIN_SRC f90 :tangle test_qmckl_distance_f.f90
+integer(c_int32_t) function test_qmckl_distance_sq(context) bind(C)
+  use qmckl
+  implicit none
+  integer(c_int64_t), intent(in), value :: context
+
+  double precision, allocatable :: A(:,:), B(:,:), C(:,:)
+  integer*8                     :: m, n, LDA, LDB, LDC
+  double precision              :: x
+  integer*8                     :: i,j 
+
+  m = 5
+  n = 6
+  LDA = m
+  LDB = n
+  LDC = 5
+
+  allocate( A(LDA,m), B(LDB,n), C(LDC,n) )
+
+  do j=1,m
+     do i=1,m
+        A(i,j) = -10.d0 + dble(i+j)
+     end do
+  end do
+  do j=1,n
+     do i=1,n
+        B(i,j) = -1.d0 + dble(i*j)
+     end do
+  end do
+
+  test_qmckl_distance_sq = qmckl_distance_sq(context, 'X', 't', m, n, A, LDA, B, LDB, C, LDC)
+  if (test_qmckl_distance_sq == 0) return 
+
+  test_qmckl_distance_sq = qmckl_distance_sq(context, 't', 'X', m, n, A, LDA, B, LDB, C, LDC)
+  if (test_qmckl_distance_sq == 0) return 
+
+  test_qmckl_distance_sq = qmckl_distance_sq(context, 'T', 't', m, n, A, LDA, B, LDB, C, LDC)
+  if (test_qmckl_distance_sq /= 0) return 
+
+  test_qmckl_distance_sq = -1
+
+  do j=1,n
+     do i=1,m
+        x =  (A(i,1)-B(j,1))**2 + &
+             (A(i,2)-B(j,2))**2 + &
+             (A(i,3)-B(j,3))**2
+        if ( dabs(1.d0 - C(i,j)/x) > 1.d-14 ) return
+     end do
+  end do
+  
+  test_qmckl_distance_sq = qmckl_distance_sq(context, 'n', 'T', m, n, A, LDA, B, LDB, C, LDC)
+  if (test_qmckl_distance_sq /= 0) return 
+
+  test_qmckl_distance_sq = -1
+
+  do j=1,n
+     do i=1,m
+        x =  (A(1,i)-B(j,1))**2 + &
+             (A(2,i)-B(j,2))**2 + &
+             (A(3,i)-B(j,3))**2
+        if ( dabs(1.d0 - C(i,j)/x) > 1.d-14 ) return
+     end do
+  end do
+
+  test_qmckl_distance_sq = qmckl_distance_sq(context, 'T', 'n', m, n, A, LDA, B, LDB, C, LDC)
+  if (test_qmckl_distance_sq /= 0) return 
+
+  test_qmckl_distance_sq = -1
+
+  do j=1,n
+     do i=1,m
+        x =  (A(i,1)-B(1,j))**2 + &
+             (A(i,2)-B(2,j))**2 + &
+             (A(i,3)-B(3,j))**2
+        if ( dabs(1.d0 - C(i,j)/x) > 1.d-14 ) return
+     end do
+  end do
+
+  test_qmckl_distance_sq = qmckl_distance_sq(context, 'n', 'N', m, n, A, LDA, B, LDB, C, LDC)
+  if (test_qmckl_distance_sq /= 0) return 
+
+  test_qmckl_distance_sq = -1
+
+  do j=1,n
+     do i=1,m
+        x =  (A(1,i)-B(1,j))**2 + &
+             (A(2,i)-B(2,j))**2 + &
+             (A(3,i)-B(3,j))**2
+        if ( dabs(1.d0 - C(i,j)/x) > 1.d-14 ) return
+     end do
+  end do
+  
+  test_qmckl_distance_sq = 0
+  
+  deallocate(A,B,C)
+end function test_qmckl_distance_sq
+      #+END_SRC
+    
+      #+BEGIN_SRC C :comments link :tangle test_qmckl_distance.c
+int test_qmckl_distance_sq(qmckl_context context);
+munit_assert_int(0, ==, test_qmckl_distance_sq(context));
+      #+END_SRC
+*** End of files                                                   :noexport:
+
+    #+BEGIN_SRC C :comments link :tangle test_qmckl_distance.c
+  if (qmckl_context_destroy(context) != QMCKL_SUCCESS)
+    return QMCKL_FAILURE;
+  return MUNIT_OK;
+}
+
+    #+END_SRC
+  
+
+ # -*- mode: org -*-
+ # vim: syntax=c

src/qmckl_footer.org

@@ -0,0 +1,18 @@
+* Acknowledgments
+
+  [[https://trex-coe.eu/sites/default/files/inline-images/euflag.jpg]]
+  [[https://trex-coe.eu][TREX: Targeting Real Chemical Accuracy at the Exascale]] project has received funding from the European Union’s Horizon 2020 - Research and Innovation program - under grant agreement no. 952165. The content of this document does not represent the opinion of the European Union, and the European Union is not responsible for any use that might be made of such content.
+
+* End of header files                                              :noexport:
+
+#+BEGIN_SRC C :tangle qmckl.h
+#endif
+#+END_SRC
+
+#+BEGIN_SRC f90 :tangle qmckl_f.f90
+end module qmckl
+#+END_SRC
+
+
+# -*- mode: org -*-
+  

src/qmckl_memory.org

@@ -0,0 +1,101 @@
+** Memory management
+
+ We override the allocation functions to enable the possibility of
+ optimized libraries to fine-tune the memory allocation.
+
+ 2 files are produced:
+ - a source file : =qmckl_memory.c=
+ - a test   file : =test_qmckl_memory.c=
+
+*** Headers                                                        :noexport:
+    #+BEGIN_SRC C :tangle qmckl_memory.c
+#include "qmckl.h"
+    #+END_SRC
+
+    #+BEGIN_SRC C :tangle test_qmckl_memory.c
+#include "qmckl.h"
+#include "munit.h"
+MunitResult test_qmckl_memory() {
+    #+END_SRC
+
+*** =qmckl_malloc=
+
+    Memory allocation function, letting the library choose how the
+    memory will be allocated, and a pointer is returned to the user.
+
+    #+BEGIN_SRC C :tangle qmckl.h
+void* qmckl_malloc(const qmckl_context ctx, const size_t size);
+    #+END_SRC
+
+    #+BEGIN_SRC f90 :tangle qmckl_f.f90
+  interface
+     type (c_ptr) function qmckl_malloc (context, size) bind(C)
+       use, intrinsic :: iso_c_binding
+       integer (c_int64_t), intent(in), value :: context
+       integer (c_int64_t), intent(in), value :: size
+     end function qmckl_malloc
+  end interface
+    #+END_SRC
+
+**** Source
+     #+BEGIN_SRC C :tangle qmckl_memory.c
+void* qmckl_malloc(const qmckl_context ctx, const size_t size) {
+  if (ctx == (qmckl_context) 0) {
+ /* Avoids unused parameter error */
+    return malloc( (size_t) size );
+  }
+  return malloc( (size_t) size );
+}
+
+     #+END_SRC
+
+**** Test                                                          :noexport:
+     #+BEGIN_SRC C :tangle test_qmckl_memory.c
+int *a;
+a = (int*) qmckl_malloc( (qmckl_context) 1, 3*sizeof(int));
+a[0] = 1;
+a[1] = 2;
+a[2] = 3;
+munit_assert_int(a[0], ==, 1);
+munit_assert_int(a[1], ==, 2);
+munit_assert_int(a[2], ==, 3);
+     #+END_SRC
+
+*** =qmckl_free=
+
+    #+BEGIN_SRC C :tangle qmckl.h
+void qmckl_free(void *ptr);
+    #+END_SRC
+
+    #+BEGIN_SRC f90 :tangle qmckl_f.f90
+  interface
+     subroutine qmckl_free (ptr) bind(C)
+       use, intrinsic :: iso_c_binding
+       type (c_ptr), intent(in), value :: ptr
+     end subroutine qmckl_free
+  end interface
+    #+END_SRC
+**** Source
+     #+BEGIN_SRC C :tangle qmckl_memory.c
+void qmckl_free(void *ptr) {
+  free(ptr);
+}
+     #+END_SRC
+
+**** Test                                                          :noexport:
+     #+BEGIN_SRC C :tangle test_qmckl_memory.c
+qmckl_free(a);
+     #+END_SRC
+
+*** End of files                                                   :noexport:
+
+**** Test
+   #+BEGIN_SRC C :comments org :tangle test_qmckl_memory.c
+  return MUNIT_OK;
+}
+
+   #+END_SRC
+ 
+
+ # -*- mode: org -*-
+ # vim: syntax=c

src/test_qmckl.org

@@ -0,0 +1,85 @@
+* QMCkl test                                                       :noexport:
+
+  This file is the main program of the unit tests. The tests rely on the
+  $\mu$unit framework, which is provided as a git submodule.
+
+  First, we use a script to find the list of all the produced test files:
+  #+NAME: test-files
+  #+BEGIN_SRC sh :exports none :results value
+grep BEGIN_SRC *.org | \
+    grep test_qmckl_ | \
+    rev | \
+    cut -d  ' ' -f 1 | \
+    rev | \
+    sort | \
+    uniq
+  #+END_SRC
+
+  #+RESULTS: test-files
+  | test_qmckl_ao.c       |
+  | test_qmckl_context.c  |
+  | test_qmckl_distance.c |
+  | test_qmckl_memory.c   |
+
+  We generate the function headers
+  #+BEGIN_SRC sh :var files=test-files :exports output :results raw
+echo "#+NAME: headers"
+echo "#+BEGIN_SRC C :tangle no"
+for file in $files
+do
+    routine=${file%.c}
+    echo "MunitResult ${routine}();"
+done
+echo "#+END_SRC"
+  #+END_SRC
+
+  #+RESULTS:
+  #+NAME: headers
+  #+BEGIN_SRC C :tangle no
+MunitResult test_qmckl_ao();
+MunitResult test_qmckl_context();
+MunitResult test_qmckl_distance();
+MunitResult test_qmckl_memory();
+  #+END_SRC
+
+  and the required function calls:
+  #+BEGIN_SRC sh :var files=test-files :exports output :results raw
+echo "#+NAME: calls"
+echo "#+BEGIN_SRC C :tangle no"
+for file in $files
+do
+    routine=${file%.c}
+    echo "  { (char*) \"${routine}\", ${routine}, NULL,NULL,MUNIT_TEST_OPTION_NONE,NULL},"
+done
+echo "#+END_SRC"
+  #+END_SRC
+  
+  #+RESULTS:
+  #+NAME: calls
+  #+BEGIN_SRC C :tangle no
+  { (char*) "test_qmckl_ao", test_qmckl_ao, NULL,NULL,MUNIT_TEST_OPTION_NONE,NULL},
+  { (char*) "test_qmckl_context", test_qmckl_context, NULL,NULL,MUNIT_TEST_OPTION_NONE,NULL},
+  { (char*) "test_qmckl_distance", test_qmckl_distance, NULL,NULL,MUNIT_TEST_OPTION_NONE,NULL},
+  { (char*) "test_qmckl_memory", test_qmckl_memory, NULL,NULL,MUNIT_TEST_OPTION_NONE,NULL},
+  #+END_SRC
+
+  #+BEGIN_SRC C :comments link :noweb yes :tangle test_qmckl.c 
+#include "qmckl.h"
+#include "munit.h"
+<<headers>>
+
+int main(int argc, char* argv[MUNIT_ARRAY_PARAM(argc + 1)]) {
+  static MunitTest test_suite_tests[] =
+    {
+     <<calls>>
+     { NULL, NULL, NULL, NULL, MUNIT_TEST_OPTION_NONE, NULL }
+    };
+
+  static const MunitSuite test_suite =
+    {
+     (char*) "", test_suite_tests, NULL, 1, MUNIT_SUITE_OPTION_NONE
+    };
+
+  return munit_suite_main(&test_suite, (void*) "µnit", argc, argv);
+  }
+  #+END_SRC