B
#+TITLE: Atomic Orbitals
#+SETUPFILE: ../tools/theme.setup
#+INCLUDE: ../tools/lib.org

* Introduction

The atomic basis set is defined as a list of shells. Each shell $s$ is
centered on a nucleus $A$, possesses a given angular momentum $l$ and a
radial function $R_s$. The radial function is a linear combination of
/primitive/ functions that can be of type Slater ($p=1$)  or
Gaussian ($p=2$):

\[
  R_s(\mathbf{r}) = \mathcal{N}_s |\mathbf{r}-\mathbf{R}_A|^{n_s}
  \sum_{k=1}^{N_{\text{prim}}} a_{ks}\, f_{ks}
 \exp \left( - \gamma_{ks} | \mathbf{r}-\mathbf{R}_A | ^p \right).
\]

In the case of Gaussian functions, $n_s$ is always zero. The
normalization factor $\mathcal{N}_s$ ensures that all the functions of
the shell are normalized (integrate) to unity. Usually, basis sets are
given a combination of normalized primitives, so the normalization
coefficients of the primitives, $f_{ks}$, need also to be provided.

Atomic orbitals (AOs) are defined as

\[
\chi_i (\mathbf{r}) = \mathcal{M}_i\, P_{\eta(i)}(\mathbf{r})\, R_{\theta(i)} (\mathbf{r})
\]

where $\theta(i)$ returns the shell on which the AO is expanded, and
$\eta(i)$ denotes which angular function is chosen and $P$ are the
generating functions of the given angular momentum $\eta(i)$.  Here,
the parameter $\mathcal{M}_i$ is an extra parameter which allows the
normalization of the different functions of the same shell to be
different, as in GAMESS for example.

In this section we describe first how the basis set is stored in the
context, and then we present the kernels used to compute the values,
gradients and Laplacian of the atomic basis functions.

* Headers                                                          :noexport:
  #+begin_src elisp :noexport :results none
(org-babel-lob-ingest "../tools/lib.org")
  #+end_src


  #+begin_src c :tangle (eval h_private_func)
#ifndef QMCKL_AO_HPF
#define QMCKL_AO_HPF

  #+end_src

  #+begin_src c :tangle (eval h_private_type)
#ifndef QMCKL_AO_HPT
#define QMCKL_AO_HPT

#include <stdbool.h>
  #+end_src

  #+begin_src f90 :tangle (eval f) :noweb yes
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
  #+end_src

  #+begin_src c :tangle (eval c_test) :noweb yes
#include "qmckl.h"
#include "assert.h"
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include <stdio.h>
#include <math.h>
#include <string.h>

#include "chbrclf.h"
#include "qmckl_ao_private_func.h"

int main() {
    qmckl_context context;
    context = qmckl_context_create();
  #+end_src

  #+begin_src c :tangle (eval c)
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#ifdef HAVE_STDINT_H
#include <stdint.h>
#elif HAVE_INTTYPES_H
#include <inttypes.h>
#endif

#include <stdlib.h>
#include <string.h>
#include <stdbool.h>
#include <assert.h>
#include <stdio.h>
#include <math.h>

#include "qmckl.h"
#include "qmckl_context_private_type.h"
#include "qmckl_memory_private_type.h"
#include "qmckl_memory_private_func.h"
#include "qmckl_ao_private_type.h"
#include "qmckl_ao_private_func.h"
  #+end_src

* Context

** Constant data

   The following arrays are stored in the context, and need to be set
   when initializing the library:

   #+NAME: constant_data
   |---------------------+----------------------+----------------------------------------------------------------------|
   | Variable            | Type                 | Description                                                          |
   |---------------------+----------------------+----------------------------------------------------------------------|
   | ~type~              | ~char~               | Gaussian (~'G'~) or Slater (~'S'~)                                   |
   | ~shell_num~         | ~int64_t~            | Number of shells                                                     |
   | ~prim_num~          | ~int64_t~            | Total number of primitives                                           |
   | ~nucleus_index~     | ~int64_t[nucl_num]~  | Index of the first shell of each nucleus                             |
   | ~nucleus_shell_num~ | ~int64_t[nucl_num]~  | Number of shells per nucleus                                         |
   | ~shell_ang_mom~     | ~int32_t[shell_num]~ | Angular momentum of each shell                                       |
   | ~shell_prim_num~    | ~int64_t[shell_num]~ | Number of primitives in each shell                                   |
   | ~shell_prim_index~  | ~int64_t[shell_num]~ | Address of the first primitive of each shell in the ~EXPONENT~ array |
   | ~shell_factor~      | ~double[shell_num]~  | Normalization factor for each shell                                  |
   | ~exponent~          | ~double[prim_num]~   | Array of exponents                                                   |
   | ~coefficient~       | ~double[prim_num]~   | Array of coefficients                                                |
   | ~prim_factor~       | ~double[prim_num]~   | Normalization factors of the primtives                               |
   | ~ao_num~            | ~int64_t~            | Number of AOs                                                        |
   | ~ao_cartesian~      | ~bool~               | If true, use polynomials. Otherwise, use spherical harmonics         |
   | ~ao_factor~         | ~double[ao_num]~     | Normalization factor of the AO                                       |

    For H_2 with the following basis set,

    #+NAME: basis
    #+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:

    #+NAME: params
    #+BEGIN_EXAMPLE
type = 'G'
shell_num = 12
prim_num = 20
ao_num = 38

nucleus_index = [0 , 6]

shell_ang_mom = [0, 0, 0, 1, 1, 2, 0, 0, 0, 1, 1, 2]
shell_factor  = [ 1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.]
shell_prim_num = [5, 1, 1, 1, 1, 1, 5, 1, 1, 1, 1, 1]
shell_prim_index = [0 , 5 , 6 , 7 , 8 , 9 , 10, 15, 16, 17, 18, 19]

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]

prim_factor = [ 1.0006253235944540e+01, 2.4169531573445120e+00, 7.9610924849766440e-01
  3.0734305383061117e-01, 1.2929684417481876e-01, 3.0734305383061117e-01,
  1.2929684417481876e-01, 2.1842769845268308e+00, 4.3649547399719840e-01,
  1.8135965626177861e+00, 1.0006253235944540e+01, 2.4169531573445120e+00,
  7.9610924849766440e-01, 3.0734305383061117e-01, 1.2929684417481876e-01,
  3.0734305383061117e-01, 1.2929684417481876e-01, 2.1842769845268308e+00,
  4.3649547399719840e-01, 1.8135965626177861e+00 ]
    #+END_EXAMPLE

   A scalar variable ~$V$~ present in this table can be set or get by
   calling the functions:

   #+NAME: template_scalar_c
   #+begin_src C :tangle no
qmckl_exit_code qmckl_set_ao_basis_$V$ ( qmckl_context context,
                                         const $type_of_V$ $V$);

qmckl_exit_code qmckl_get_ao_basis_$V$ ( const qmckl_context context,
                                         $type_of_V$ const $V$);
   #+end_src

   #+begin_src f90 :tangle no
interface
   integer(c_int32_t) function qmckl_set_ao_basis_$V$ (context, $V$) &
        bind(C)
     use, intrinsic :: iso_c_binding
     import
     implicit none
     integer (c_int64_t) , intent(in)  , value :: context
     $f_type_of_V$       , intent(in)  , value :: $V$
   end function qmckl_set_ao_basis_$V$
end interface

interface
   integer(c_int32_t) function qmckl_get_ao_basis_$V$ (context, $V$) &
        bind(C)
     use, intrinsic :: iso_c_binding
     import
     implicit none
     integer (c_int64_t) , intent(in)  , value :: context
     $f_type_of_V$       , intent(out)         :: $V$
   end function qmckl_get_ao_basis_$V$
end interface
   #+end_src
   For array variables, use the rule:

   #+NAME: template_array_c
   #+begin_src C :tangle no
qmckl_exit_code qmckl_set_ao_basis_$V$ ( qmckl_context context,
                                         const $type_of_V$ $V$,
                                         const int64_t size_max);

qmckl_exit_code qmckl_get_ao_basis_$V$ ( const qmckl_context context,
                                         $type_of_V$ const $V$,
                                         const int64_t size_max);
   #+end_src

   #+begin_src f90 :tangle no
interface
   integer(c_int32_t) function qmckl_set_ao_basis_$V$ (context, &
        $V$, size_max) bind(C)
     use, intrinsic :: iso_c_binding
     import
     implicit none
     integer (c_int64_t) , intent(in)  , value :: context
     $f_type_of_V$       , intent(in)  , value :: $V$
     integer (c_int64_t) , intent(in)  , value :: size_max
   end function qmckl_set_ao_basis_$V$
end interface

interface
   integer(c_int32_t) function qmckl_get_ao_basis_$V$ (context, &
        $V$, size_max) bind(C)
     use, intrinsic :: iso_c_binding
     import
     implicit none
     integer (c_int64_t) , intent(in)  , value :: context
     $f_type_of_V$       , intent(out)         :: $V$
     integer (c_int64_t) , intent(in)  , value :: size_max
   end function qmckl_get_ao_basis_$V$
end interface
   #+end_src

*** Data structure                                                 :noexport:

    #+begin_src c :comments org :tangle (eval h_private_type) :noweb yes
typedef struct qmckl_ao_basis_struct {
  int64_t   shell_num;
  int64_t   prim_num;
  int64_t   ao_num;
  int64_t * nucleus_index;
  int64_t * nucleus_shell_num;
  int32_t * shell_ang_mom;
  int64_t * shell_prim_num;
  int64_t * shell_prim_index;
  double  * shell_factor;
  double  * exponent;
  double  * coefficient;
  double  * prim_factor;
  double  * ao_factor;

  int64_t * nucleus_prim_index;
  double  * coefficient_normalized;
  int32_t * nucleus_max_ang_mom;
  double  * nucleus_range;
  double  * primitive_vgl;
  uint64_t  primitive_vgl_date;
  double  * shell_vgl;
  uint64_t  shell_vgl_date;
  double  * ao_vgl;
  uint64_t  ao_vgl_date;

  int32_t   uninitialized;
  bool      provided;
  bool      ao_cartesian;
  char      type;
#ifdef HAVE_HPC
  <<HPC_struct>>
#endif
} qmckl_ao_basis_struct;
    #+end_src

    The ~uninitialized~ integer contains one bit set to one for each
    initialization function which has not been called. It becomes equal
    to zero after all initialization functions have been called. The
    struct is then initialized and ~provided == true~.
    Some values are initialized by default, and are not concerned by
    this mechanism.

    #+begin_src c :comments org :tangle (eval h_private_func)
qmckl_exit_code qmckl_init_ao_basis(qmckl_context context);
    #+end_src

    #+begin_src c :comments org :tangle (eval c)
qmckl_exit_code qmckl_init_ao_basis(qmckl_context context) {

  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_CONTEXT,
                           "qmckl_init_ao_basis",
                           NULL);
  }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

  ctx->ao_basis.uninitialized = (1 << 14) - 1;

  /* Default values */
  ctx->ao_basis.ao_cartesian = true;

  return QMCKL_SUCCESS;
}
    #+end_src

*** Initialization functions

    ~size_max~ is the dimension of the input array, which should be
    equal of larger than the value given in the table of section [[Context]].

**** C interface

     #+NAME:pre2
     #+begin_src c  :exports none
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_CONTEXT,
                           "qmckl_get_ao_basis_*",
                           NULL);
 }

qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
     #+end_src

     #+NAME:post2
     #+begin_src c  :exports none
ctx->ao_basis.uninitialized &= ~mask;
ctx->ao_basis.provided = (ctx->ao_basis.uninitialized == 0);
if (ctx->ao_basis.provided) {
  qmckl_exit_code rc_ = qmckl_finalize_basis(context);
  if (rc_ != QMCKL_SUCCESS) return rc_;
 }

return QMCKL_SUCCESS;
     #+end_src

     To set the basis set, all the following functions need to be
     called.

     #+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_set_ao_basis_type (qmckl_context context,
                         const char basis_type);
     #+end_src

     #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_set_ao_basis_type(qmckl_context context,
                        const char basis_type)
{
  <<pre2>>

  if (basis_type != 'G' && basis_type != 'S') {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_set_ao_basis_type",
                           NULL);
  }

  int32_t mask = 1;
  ctx->ao_basis.type = basis_type;

  <<post2>>
}
     #+end_src


     #+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_set_ao_basis_shell_num (qmckl_context context,
                              const int64_t shell_num);
     #+end_src

     #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_set_ao_basis_shell_num (qmckl_context context,
                              const int64_t shell_num)
{
  <<pre2>>

  if (shell_num <= 0) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_set_ao_basis_shell_num",
                           "shell_num <= 0");
  }

  const int64_t prim_num = ctx->ao_basis.prim_num;

  if (0L < prim_num && prim_num < shell_num) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_set_ao_basis_shell_num",
                           "shell_num > prim_num");
  }

  int32_t mask = 1 << 1;
  ctx->ao_basis.shell_num = shell_num;

  <<post2>>
}
     #+end_src


     #+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_set_ao_basis_prim_num (qmckl_context context,
                             const int64_t prim_num);
     #+end_src

     #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_set_ao_basis_prim_num (qmckl_context context,
                             const int64_t prim_num)
{
  <<pre2>>

  if (prim_num <= 0) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_set_ao_basis_shell_num",
                           "prim_num must be positive");
  }

  const int64_t shell_num = ctx->ao_basis.shell_num;

  if (shell_num <= 0L) {
    return qmckl_failwith( context,
                           QMCKL_FAILURE,
                           "qmckl_set_ao_basis_shell_num",
                           "shell_num is not set");
  }

  if (prim_num < shell_num) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_set_ao_basis_shell_num",
                           "prim_num < shell_num");
  }

  int32_t mask = 1 << 2;
  ctx->ao_basis.prim_num = prim_num;

  <<post2>>
}
     #+end_src


     #+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_set_ao_basis_ao_num (qmckl_context context,
                           const int64_t ao_num);
     #+end_src

     #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_set_ao_basis_ao_num (qmckl_context context,
                           const int64_t ao_num)
{
  <<pre2>>

  if (ao_num <= 0) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_set_ao_basis_shell_num",
                           "ao_num must be positive");
  }

  const int64_t shell_num = ctx->ao_basis.shell_num;
  if (shell_num <= 0L) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_set_ao_basis_shell_num",
                           "shell_num is not set");
  }

  if (ao_num < shell_num) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_set_ao_basis_shell_num",
                           "ao_num < shell_num");
  }

  int32_t mask = 1 << 12;
  ctx->ao_basis.ao_num = ao_num;

  <<post2>>
}
     #+end_src


     #+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_set_ao_basis_nucleus_index (qmckl_context context,
                                  const int64_t* nucleus_index,
                                  const int64_t size_max);
     #+end_src

     #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_set_ao_basis_nucleus_index (qmckl_context context,
                                  const int64_t* nucleus_index,
                                  const int64_t size_max)
{
  <<pre2>>

  int32_t mask = 1 << 4;

  const int64_t nucl_num = ctx->nucleus.num;

  if (nucl_num <= 0L) {
    return qmckl_failwith( context,
                           QMCKL_FAILURE,
                           "qmckl_set_ao_basis_nucleus_index",
                           "nucl_num is not set");
  }

  if (size_max < nucl_num) {
    return qmckl_failwith( context,
                           QMCKL_FAILURE,
                           "qmckl_set_ao_basis_nucleus_index",
                           "input array too small");
  }

  if (ctx->ao_basis.nucleus_index != NULL) {
    qmckl_exit_code rc = qmckl_free(context, ctx->ao_basis.nucleus_index);
    if (rc != QMCKL_SUCCESS) {
      return qmckl_failwith( context, rc,
                             "qmckl_set_ao_basis_nucleus_index",
                             NULL);
    }
  }

  qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
  mem_info.size = nucl_num * sizeof(int64_t);
  int64_t* new_array = (int64_t*) qmckl_malloc(context, mem_info);

  if (new_array == NULL) {
    return qmckl_failwith( context,
                           QMCKL_ALLOCATION_FAILED,
                           "qmckl_set_ao_basis_nucleus_index",
                           NULL);
  }

  memcpy(new_array, nucleus_index, mem_info.size);

  ctx->ao_basis.nucleus_index = new_array;

  <<post2>>
}
     #+end_src


     #+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_set_ao_basis_nucleus_shell_num (qmckl_context context,
                                      const int64_t* nucleus_shell_num,
                                      const int64_t size_max);
     #+end_src

     #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_set_ao_basis_nucleus_shell_num (qmckl_context context,
                                      const int64_t* nucleus_shell_num,
                                      const int64_t size_max)
{
  <<pre2>>

  int32_t mask = 1 << 3;

  const int64_t nucl_num = ctx->nucleus.num;

  if (nucl_num <= 0L) {
    return qmckl_failwith( context,
                           QMCKL_FAILURE,
                           "qmckl_set_ao_basis_nucleus_shell_num",
                           "shell_num is not set");
  }

  if (size_max < nucl_num) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_3,
                           "qmckl_set_ao_basis_nucleus_shell_num",
                           "input array too small");
  }

  if (ctx->ao_basis.nucleus_shell_num != NULL) {
    qmckl_exit_code rc = qmckl_free(context, ctx->ao_basis.nucleus_shell_num);
    if (rc != QMCKL_SUCCESS) {
      return qmckl_failwith( context, rc,
                             "qmckl_set_ao_basis_nucleus_shell_num",
                             NULL);
    }
  }

  qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
  mem_info.size = nucl_num * sizeof(int64_t);
  int64_t* new_array = (int64_t*) qmckl_malloc(context, mem_info);

  if (new_array == NULL) {
    return qmckl_failwith( context,
                           QMCKL_ALLOCATION_FAILED,
                           "qmckl_set_ao_basis_nucleus_shell_num",
                           NULL);
  }

  memcpy(new_array, nucleus_shell_num, mem_info.size);

  ctx->ao_basis.nucleus_shell_num = new_array;

  <<post2>>
}
     #+end_src


     #+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_set_ao_basis_shell_ang_mom (qmckl_context context,
                                  const int32_t* shell_ang_mom,
                                  const int64_t size_max);
     #+end_src

     #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_set_ao_basis_shell_ang_mom (qmckl_context context,
                                  const int32_t* shell_ang_mom,
                                  const int64_t size_max)
{
  <<pre2>>

  int32_t mask = 1 << 5;

  const int64_t shell_num = ctx->ao_basis.shell_num;

  if (shell_num == 0L) {
    return qmckl_failwith( context,
                           QMCKL_FAILURE,
                           "qmckl_set_ao_basis_shell_ang_mom",
                           "shell_num is not set");
  }

  if (size_max < shell_num) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_3,
                           "qmckl_set_ao_basis_shell_ang_mom",
                           "input array too small");
  }

  if (ctx->ao_basis.shell_ang_mom != NULL) {
    qmckl_exit_code rc = qmckl_free(context, ctx->ao_basis.shell_ang_mom);
    if (rc != QMCKL_SUCCESS) {
      return qmckl_failwith( context, rc,
                             "qmckl_set_ao_basis_shell_ang_mom",
                             NULL);
    }
  }


  qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
  mem_info.size = shell_num * sizeof(int32_t);
  int32_t * new_array = (int32_t*) qmckl_malloc(context, mem_info);

  if (new_array == NULL) {
    return qmckl_failwith( context,
                           QMCKL_ALLOCATION_FAILED,
                           "qmckl_set_ao_basis_shell_ang_mom",
                           NULL);
  }

  memcpy(new_array, shell_ang_mom, mem_info.size);

  ctx->ao_basis.shell_ang_mom = new_array;

  <<post2>>
}
     #+end_src


     #+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_set_ao_basis_shell_prim_num (qmckl_context context,
                                   const int64_t* shell_prim_num,
                                   const int64_t size_max);
     #+end_src

     #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_set_ao_basis_shell_prim_num (qmckl_context context,
                                   const int64_t* shell_prim_num,
                                   const int64_t size_max)
{
  <<pre2>>

  int32_t mask = 1 << 6;

  const int64_t shell_num = ctx->ao_basis.shell_num;

  if (shell_num <= 0L) {
    return qmckl_failwith( context,
                           QMCKL_FAILURE,
                           "qmckl_set_ao_basis_shell_prim_num",
                           "shell_num is not set");
  }

  if (size_max < shell_num) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_3,
                           "qmckl_set_ao_basis_shell_prim_num",
                           "input array too small");
  }

  if (ctx->ao_basis.shell_prim_num != NULL) {
    qmckl_exit_code rc = qmckl_free(context, ctx->ao_basis.shell_prim_num);
    if (rc != QMCKL_SUCCESS) {
      return qmckl_failwith( context, rc,
                             "qmckl_set_ao_basis_shell_prim_num",
                             NULL);
    }
  }


  qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
  mem_info.size = shell_num * sizeof(int64_t);
  int64_t* new_array = (int64_t*) qmckl_malloc(context, mem_info);

  if (new_array == NULL) {
    return qmckl_failwith( context,
                           QMCKL_ALLOCATION_FAILED,
                           "qmckl_set_ao_basis_shell_prim_num",
                           NULL);
  }

  memcpy(new_array, shell_prim_num, mem_info.size);

  ctx->ao_basis.shell_prim_num = new_array;

  <<post2>>
}
     #+end_src

     
     #+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_set_ao_basis_shell_prim_index (qmckl_context context,
                                     const int64_t* shell_prim_index,
                                     const int64_t size_max);
     #+end_src
     
     #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_set_ao_basis_shell_prim_index (qmckl_context context,
                                     const int64_t* shell_prim_index,
                                     const int64_t size_max)
{
  <<pre2>>

  int32_t mask = 1 << 7;

  const int64_t shell_num = ctx->ao_basis.shell_num;

  if (shell_num <= 0L) {
    return qmckl_failwith( context,
                           QMCKL_FAILURE,
                           "qmckl_set_ao_basis_shell_prim_index",
                           "shell_num is not set");
  }

  if (size_max < shell_num) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_3,
                           "qmckl_set_ao_basis_shell_prim_index",
                           "input array too small");
  }

  if (ctx->ao_basis.shell_prim_index != NULL) {
    qmckl_exit_code rc = qmckl_free(context, ctx->ao_basis.shell_prim_index);
    if (rc != QMCKL_SUCCESS) {
      return qmckl_failwith( context, rc,
                             "qmckl_set_ao_basis_shell_prim_index",
                             NULL);
    }
  }

  qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
  mem_info.size = shell_num * sizeof(int64_t);
  int64_t* new_array = (int64_t*) qmckl_malloc(context, mem_info);

  if (new_array == NULL) {
    return qmckl_failwith( context,
                           QMCKL_ALLOCATION_FAILED,
                           "qmckl_set_ao_basis_shell_prim_index",
                           NULL);
  }

  memcpy(new_array, shell_prim_index, mem_info.size);

  ctx->ao_basis.shell_prim_index = new_array;

  <<post2>>
}
     #+end_src
     
     
     #+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_set_ao_basis_shell_factor (qmckl_context context,
                                 const double* shell_factor,
                                 const int64_t size_max);
     #+end_src
     
     #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_set_ao_basis_shell_factor (qmckl_context context,
                                 const double* shell_factor,
                                 const int64_t size_max)
{
  <<pre2>>

  int32_t mask = 1 << 8;

  const int64_t shell_num = ctx->ao_basis.shell_num;

  if (shell_num <= 0L) {
    return qmckl_failwith( context,
                           QMCKL_FAILURE,
                           "qmckl_set_ao_basis_shell_factor",
                           "shell_num is not set");
  }

  if (size_max < shell_num) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_3,
                           "qmckl_set_ao_basis_shell_factor",
                           "input array too small");
  }

  if (ctx->ao_basis.shell_factor != NULL) {
    qmckl_exit_code rc = qmckl_free(context, ctx->ao_basis.shell_factor);
    if (rc != QMCKL_SUCCESS) {
      return qmckl_failwith( context, rc,
                             "qmckl_set_ao_basis_shell_factor",
                             NULL);
    }
  }

  qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
  mem_info.size = shell_num * sizeof(double);
  double* new_array = (double*) qmckl_malloc(context, mem_info);

  if (new_array == NULL) {
    return qmckl_failwith( context,
                           QMCKL_ALLOCATION_FAILED,
                           "qmckl_set_ao_basis_shell_factor",
                           NULL);
  }

  memcpy(new_array, shell_factor, mem_info.size);

  ctx->ao_basis.shell_factor = new_array;

  <<post2>>
}
     #+end_src
     
     
     #+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_set_ao_basis_exponent (qmckl_context context,
                             const double* exponent,
                             const int64_t size_max);
     #+end_src
     
     #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_set_ao_basis_exponent (qmckl_context context,
                             const double* exponent,
                             const int64_t size_max)
{
  <<pre2>>

  int32_t mask = 1 << 9;

  const int64_t prim_num = ctx->ao_basis.prim_num;

  if (prim_num <= 0L) {
    return qmckl_failwith( context,
                           QMCKL_FAILURE,
                           "qmckl_set_ao_basis_exponent",
                           "prim_num is not set");
  }

  if (size_max < prim_num) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_3,
                           "qmckl_set_ao_basis_exponent",
                           "input array too small");
  }

  if (ctx->ao_basis.exponent != NULL) {
    qmckl_exit_code rc = qmckl_free(context, ctx->ao_basis.exponent);
    if (rc != QMCKL_SUCCESS) {
      return qmckl_failwith( context, rc,
                             "qmckl_set_ao_basis_exponent",
                             NULL);
    }
  }

  qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
  mem_info.size = prim_num * sizeof(double);
  double* new_array = (double*) qmckl_malloc(context, mem_info);

  if (new_array == NULL) {
    return qmckl_failwith( context,
                           QMCKL_ALLOCATION_FAILED,
                           "qmckl_set_ao_basis_exponent",
                           NULL);
  }

  memcpy(new_array, exponent, mem_info.size);

  ctx->ao_basis.exponent = new_array;

  <<post2>>
}
     #+end_src


     #+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_set_ao_basis_coefficient (qmckl_context context,
                                const double* coefficient,
                                const int64_t size_max);
     #+end_src

     #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_set_ao_basis_coefficient (qmckl_context context,
                                const double* coefficient,
                                const int64_t size_max)
{
  <<pre2>>

  int32_t mask = 1 << 10;

  const int64_t prim_num = ctx->ao_basis.prim_num;

  if (prim_num <= 0L) {
    return qmckl_failwith( context,
                           QMCKL_FAILURE,
                           "qmckl_set_ao_basis_coefficient",
                           "prim_num is not set");
  }

  if (size_max < prim_num) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_3,
                           "qmckl_set_ao_basis_coefficient",
                           "input array too small");
  }

  if (ctx->ao_basis.coefficient != NULL) {
    qmckl_exit_code rc = qmckl_free(context, ctx->ao_basis.coefficient);
    if (rc != QMCKL_SUCCESS) {
      return qmckl_failwith( context, rc,
                             "qmckl_set_ao_basis_coefficient",
                             NULL);
    }
  }

  qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
  mem_info.size = prim_num * sizeof(double);
  double* new_array = (double*) qmckl_malloc(context, mem_info);

  if (new_array == NULL) {
    return qmckl_failwith( context,
                           QMCKL_ALLOCATION_FAILED,
                           "qmckl_set_ao_basis_coefficient",
                           NULL);
  }

  memcpy(new_array, coefficient, mem_info.size);

  ctx->ao_basis.coefficient = new_array;

  <<post2>>
}
     #+end_src


     #+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_set_ao_basis_prim_factor (qmckl_context context,
                                const double* prim_factor,
                                const int64_t size_max);
     #+end_src

     #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_set_ao_basis_prim_factor (qmckl_context context,
                                const double* prim_factor,
                                const int64_t size_max)
{
  <<pre2>>

  int32_t mask = 1 << 11;

  const int64_t prim_num = ctx->ao_basis.prim_num;

  if (prim_num <= 0L) {
    return qmckl_failwith( context,
                           QMCKL_FAILURE,
                           "qmckl_set_ao_basis_prim_factor",
                           "prim_num is not set");
  }

  if (size_max < prim_num) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_3,
                           "qmckl_set_ao_basis_prim_factor",
                           "input array too small");
  }

  if (ctx->ao_basis.prim_factor != NULL) {
    qmckl_exit_code rc = qmckl_free(context, ctx->ao_basis.prim_factor);
    if (rc != QMCKL_SUCCESS) {
      return qmckl_failwith( context, rc,
                             "qmckl_set_ao_basis_prim_factor",
                             NULL);
    }
  }

  qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
  mem_info.size = prim_num * sizeof(double);
  double* new_array = (double*) qmckl_malloc(context, mem_info);

  if (new_array == NULL) {
    return qmckl_failwith( context,
                           QMCKL_ALLOCATION_FAILED,
                           "qmckl_set_ao_basis_prim_factor",
                           NULL);
  }

  memcpy(new_array, prim_factor, mem_info.size);

  ctx->ao_basis.prim_factor = new_array;

  <<post2>>
}
     #+end_src


     #+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_set_ao_basis_ao_factor (qmckl_context context,
                              const double* ao_factor,
                              const int64_t size_max);
     #+end_src

     #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_set_ao_basis_ao_factor (qmckl_context context,
                              const double* ao_factor,
                              const int64_t size_max)
{
  <<pre2>>

  int32_t mask = 1 << 13;

  const int64_t ao_num = ctx->ao_basis.ao_num;

  if (ao_num <= 0L) {
    return qmckl_failwith( context,
                           QMCKL_FAILURE,
                           "qmckl_set_ao_basis_ao_factor",
                           "ao_num is not set");
  }

  if (size_max < ao_num) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_3,
                           "qmckl_set_ao_basis_ao_factor",
                           "input array too small");
  }

  if (ctx->ao_basis.ao_factor != NULL) {
    qmckl_exit_code rc = qmckl_free(context, ctx->ao_basis.ao_factor);
    if (rc != QMCKL_SUCCESS) {
      return qmckl_failwith( context, rc,
                             "qmckl_set_ao_basis_ao_factor",
                             NULL);
    }
  }

  qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
  mem_info.size = ao_num * sizeof(double);
  double* new_array = (double*) qmckl_malloc(context, mem_info);

  if (new_array == NULL) {
    return qmckl_failwith( context,
                           QMCKL_ALLOCATION_FAILED,
                           "qmckl_set_ao_basis_ao_factor",
                           NULL);
  }

  memcpy(new_array, ao_factor, mem_info.size);

  ctx->ao_basis.ao_factor = new_array;

  <<post2>>
}
     #+end_src


     #+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_set_ao_basis_cartesian (qmckl_context context,
                              const bool cartesian);
     #+end_src

     #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_set_ao_basis_cartesian (qmckl_context context,
                              const bool cartesian)
{
  <<pre2>>

  int32_t mask = 1;
  ctx->ao_basis.ao_cartesian = cartesian;

  <<post2>>
}
     #+end_src

**** Fortran interface

      #+begin_src f90 :tangle (eval fh_func) :comments org
interface
  integer(c_int32_t) function qmckl_set_ao_basis_type (context, &
       basis_type) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    character(c_char)   , intent(in)  , value :: basis_type
  end function qmckl_set_ao_basis_type
end interface

interface
  integer(c_int32_t) function qmckl_set_ao_basis_shell_num(context, &
       num) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    integer (c_int64_t) , intent(in)  , value :: num
  end function qmckl_set_ao_basis_shell_num
end interface

interface
  integer(c_int32_t) function qmckl_set_ao_basis_prim_num(context, &
       num) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    integer (c_int64_t) , intent(in)  , value :: num
  end function qmckl_set_ao_basis_prim_num
end interface

interface
  integer(c_int32_t) function qmckl_set_ao_basis_nucleus_index(context, &
       idx, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    integer (c_int64_t) , intent(in)          :: idx(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_set_ao_basis_nucleus_index
end interface

interface
  integer(c_int32_t) function qmckl_set_ao_basis_nucleus_shell_num(context, &
       shell_num, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    integer (c_int64_t) , intent(in)          :: shell_num(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_set_ao_basis_nucleus_shell_num
end interface

interface
  integer(c_int32_t) function qmckl_set_ao_basis_shell_ang_mom(context, &
       shell_ang_mom, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    integer (c_int32_t) , intent(in)          :: shell_ang_mom(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_set_ao_basis_shell_ang_mom
end interface

interface
  integer(c_int32_t) function qmckl_set_ao_basis_shell_prim_num(context, &
       shell_prim_num, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    integer (c_int64_t) , intent(in)          :: shell_prim_num(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_set_ao_basis_shell_prim_num
end interface

interface
  integer(c_int32_t) function qmckl_set_ao_basis_shell_prim_index(context, &
       shell_prim_index, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    integer (c_int64_t) , intent(in)          :: shell_prim_index(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_set_ao_basis_shell_prim_index
end interface

interface
  integer(c_int32_t) function qmckl_set_ao_basis_shell_factor(context, &
       shell_factor, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    real    (c_double)  , intent(in)          :: shell_factor(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_set_ao_basis_shell_factor
end interface

interface
  integer(c_int32_t) function qmckl_set_ao_basis_exponent(context, &
       exponent, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    real    (c_double)  , intent(in)          :: exponent(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_set_ao_basis_exponent
end interface

interface
  integer(c_int32_t) function qmckl_set_ao_basis_coefficient(context, &
       coefficient, size_max)  bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    real    (c_double)  , intent(in)          :: coefficient(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_set_ao_basis_coefficient
end interface

interface
  integer(c_int32_t) function qmckl_set_ao_basis_prim_factor(context, &
       prim_factor, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    real    (c_double)  , intent(in)          :: prim_factor(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_set_ao_basis_prim_factor
end interface

interface
  integer(c_int32_t) function qmckl_set_ao_basis_ao_num(context, &
       num) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    integer (c_int64_t) , intent(in)  , value :: num
  end function qmckl_set_ao_basis_ao_num
end interface

interface
  integer(c_int32_t) function qmckl_set_ao_basis_cartesian(context, &
       cartesian) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    logical (c_bool)    , intent(in)  , value :: cartesian
  end function qmckl_set_ao_basis_cartesian
end interface

interface
  integer(c_int32_t) function qmckl_set_ao_basis_ao_factor(context, &
       ao_factor, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    real    (c_double)  , intent(in)          :: ao_factor(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_set_ao_basis_ao_factor
end interface
     #+end_src

*** Access functions

    ~size_max~ is the dimension of the input array, which should be
    equal of larger than the value given in the table of section [[Context]].

**** C interface

    #+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_get_ao_basis_type (const qmckl_context context,
                         char* const basis_type);
    #+end_src

    #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_ao_basis_type (const qmckl_context context,
                         char* const basis_type)
{

  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_CONTEXT,
                           "qmckl_get_ao_basis_type",
                           NULL);
  }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

  int32_t mask = 1;

  if ( (ctx->ao_basis.uninitialized & mask) != 0) {
    return qmckl_failwith( context,
                           QMCKL_NOT_PROVIDED,
                           "qmckl_get_ao_basis_type",
                           NULL);

  }

  assert (ctx->ao_basis.type != (char) 0);

  basis_type[0] = ctx->ao_basis.type;
  return QMCKL_SUCCESS;
}

    #+end_src


    #+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_get_ao_basis_shell_num (const qmckl_context context,
                              int64_t* const shell_num);
    #+end_src

    #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_ao_basis_shell_num (const qmckl_context context,
                              int64_t* const shell_num)
{
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_CONTEXT,
                           "qmckl_get_ao_basis_shell_factor",
                           NULL);
  }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

  int32_t mask = 1 << 1;

  if ( (ctx->ao_basis.uninitialized & mask) != 0) {
    return qmckl_failwith( context,
                           QMCKL_NOT_PROVIDED,
                           "qmckl_get_ao_basis_shell_num",
                           NULL);
  }

  assert (ctx->ao_basis.shell_num > (int64_t) 0);
  ,*shell_num = ctx->ao_basis.shell_num;
  return QMCKL_SUCCESS;
}

    #+end_src


    #+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_get_ao_basis_prim_num (const qmckl_context context,
                             int64_t* const prim_num);
    #+end_src

    #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_ao_basis_prim_num (const qmckl_context context,
                             int64_t* const prim_num)
{
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_CONTEXT,
                           "qmckl_get_ao_basis_prim_num",
                           NULL);
  }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

  int32_t mask = 1 << 2;

  if ( (ctx->ao_basis.uninitialized & mask) != 0) {
    return qmckl_failwith( context,
                           QMCKL_NOT_PROVIDED,
                           "qmckl_get_ao_basis_prim_num",
                           NULL);
  }

  assert (ctx->ao_basis.prim_num > (int64_t) 0);

  ,*prim_num = ctx->ao_basis.prim_num;
  return QMCKL_SUCCESS;
}

    #+end_src


    #+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_get_ao_basis_nucleus_shell_num (const qmckl_context context,
                                      int64_t* const nucleus_shell_num,
                                      const int64_t size_max);
    #+end_src

    #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_ao_basis_nucleus_shell_num (const qmckl_context context,
                                      int64_t* const nucleus_shell_num,
                                      const int64_t size_max)
{
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_CONTEXT,
                           "qmckl_get_ao_basis_nucleus_shell_num",
                           NULL);
  }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

  int32_t mask = 1 << 3;

  if ( (ctx->ao_basis.uninitialized & mask) != 0) {
    return qmckl_failwith( context,
                           QMCKL_NOT_PROVIDED,
                           "qmckl_get_ao_basis_nucleus_shell_num",
                           NULL);
  }

  if (nucleus_shell_num == NULL) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_get_ao_basis_nucleus_shell_num",
                           "NULL pointer");
  }

  if (size_max < ctx->nucleus.num) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_3,
                           "qmckl_get_ao_basis_nucleus_shell_num",
                           "Array too small. Expected nucl_num");
  }

  assert (ctx->ao_basis.nucleus_shell_num != NULL);
  memcpy(nucleus_shell_num, ctx->ao_basis.nucleus_shell_num,
         (size_t) ctx->nucleus.num * sizeof(int64_t));
  return QMCKL_SUCCESS;
}

    #+end_src


    #+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_get_ao_basis_nucleus_index (const qmckl_context context,
                                  int64_t* const nucleus_index,
                                  const int64_t size_max);
    #+end_src

    #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_ao_basis_nucleus_index (const qmckl_context context,
                                  int64_t* const nucleus_index,
                                  const int64_t size_max)
{
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_CONTEXT,
                           "qmckl_get_ao_basis_nucleus_index",
                           NULL);
  }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

  int32_t mask = 1 << 4;

  if ( (ctx->ao_basis.uninitialized & mask) != 0) {
    return qmckl_failwith( context,
                           QMCKL_NOT_PROVIDED,
                           "qmckl_get_ao_basis_nucleus_index",
                           NULL);
  }

  if (nucleus_index == NULL) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_get_ao_basis_nucleus_index",
                           "NULL pointer");
  }

  if (size_max < ctx->nucleus.num) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_3,
                           "qmckl_get_ao_basis_nucleus_index",
                           "Array too small. Expected shell_num");
  }

  assert (ctx->ao_basis.nucleus_index != NULL);
  memcpy(nucleus_index, ctx->ao_basis.nucleus_index,
         (size_t) ctx->nucleus.num * sizeof(int64_t));
  return QMCKL_SUCCESS;
}


    #+end_src


    #+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_get_ao_basis_shell_ang_mom (const qmckl_context context,
                                  int32_t* const shell_ang_mom,
                                  const int64_t size_max);
    #+end_src

    #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_ao_basis_shell_ang_mom (const qmckl_context context,
                                  int32_t* const shell_ang_mom,
                                  const int64_t size_max)
{
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_CONTEXT,
                           "qmckl_get_ao_basis_shell_ang_mom",
                           NULL);
  }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

  int32_t mask = 1 << 5;

  if ( (ctx->ao_basis.uninitialized & mask) != 0) {
    return qmckl_failwith( context,
                           QMCKL_NOT_PROVIDED,
                           "qmckl_get_ao_basis_shell_ang_mom",
                           NULL);
  }

  if (shell_ang_mom == NULL) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_get_ao_basis_shell_ang_mom",
                           "NULL pointer");
  }

  if (size_max < ctx->ao_basis.shell_num) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_3,
                           "qmckl_get_ao_basis_shell_ang_mom",
                           "Array too small. Expected shell_num");
  }

  assert (ctx->ao_basis.shell_ang_mom != NULL);
  memcpy(shell_ang_mom, ctx->ao_basis.shell_ang_mom,
         (size_t) ctx->ao_basis.shell_num * sizeof(int32_t));
  return QMCKL_SUCCESS;
}


    #+end_src


    #+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_get_ao_basis_shell_prim_num (const qmckl_context context,
                                   int64_t* const shell_prim_num,
                                   const int64_t size_max);
    #+end_src

    #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_ao_basis_shell_prim_num (const qmckl_context context,
                                   int64_t* const shell_prim_num,
                                   const int64_t size_max)
{
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_CONTEXT,
                           "qmckl_get_ao_basis_shell_prim_num",
                           NULL);
  }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

  int32_t mask = 1 << 6;

  if ( (ctx->ao_basis.uninitialized & mask) != 0) {
    return qmckl_failwith( context,
                           QMCKL_NOT_PROVIDED,
                           "qmckl_get_ao_basis_shell_prim_num",
                           NULL);
  }

  if (shell_prim_num == NULL) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_get_ao_basis_shell_prim_num",
                           "NULL pointer");
  }

  if (size_max < ctx->ao_basis.shell_num) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_3,
                           "qmckl_get_ao_basis_shell_prim_num",
                           "Array too small. Expected shell_num");
  }

  assert (ctx->ao_basis.shell_prim_num != NULL);
  memcpy(shell_prim_num, ctx->ao_basis.shell_prim_num,
         (size_t) ctx->ao_basis.shell_num * sizeof(int64_t));
  return QMCKL_SUCCESS;
}


    #+end_src


    #+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_get_ao_basis_shell_prim_index (const qmckl_context context,
                                     int64_t* const shell_prim_index,
                                     const int64_t size_max);
    #+end_src

    #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_ao_basis_shell_prim_index (const qmckl_context context,
                                     int64_t* const shell_prim_index,
                                     const int64_t size_max)
{
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_CONTEXT,
                           "qmckl_get_ao_basis_shell_prim_index",
                           NULL);
  }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

  int32_t mask = 1 << 7;

  if ( (ctx->ao_basis.uninitialized & mask) != 0) {
    return qmckl_failwith( context,
                           QMCKL_NOT_PROVIDED,
                           "qmckl_get_ao_basis_shell_prim_index",
                           NULL);
  }

  if (shell_prim_index == NULL) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_get_ao_basis_shell_prim_index",
                           "NULL pointer");
  }

  if (size_max < ctx->ao_basis.shell_num) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_3,
                           "qmckl_get_ao_basis_shell_prim_index",
                           "Array too small. Expected shell_num");
  }

  assert (ctx->ao_basis.shell_prim_index != NULL);
  memcpy(shell_prim_index, ctx->ao_basis.shell_prim_index,
         (size_t) ctx->ao_basis.shell_num * sizeof(int64_t));
  return QMCKL_SUCCESS;
}
    #+end_src


    #+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_get_ao_basis_shell_factor (const qmckl_context context,
                                 double*  const shell_factor,
                                 const int64_t size_max);
    #+end_src

    #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_ao_basis_shell_factor (const qmckl_context context,
                                 double* const shell_factor,
                                 const int64_t size_max)
{
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_CONTEXT,
                           "qmckl_get_ao_basis_shell_factor",
                           NULL);
  }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

  int32_t mask = 1 << 8;

  if ( (ctx->ao_basis.uninitialized & mask) != 0) {
    return qmckl_failwith( context,
                           QMCKL_NOT_PROVIDED,
                           "qmckl_get_ao_basis_shell_factor",
                           NULL);
  }

  if (shell_factor == NULL) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_get_ao_basis_shell_factor",
                           "NULL pointer");
  }

  if (size_max < ctx->ao_basis.shell_num) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_3,
                           "qmckl_get_ao_basis_shell_factor",
                           "Array too small. Expected shell_num");
  }

  assert (ctx->ao_basis.shell_factor != NULL);
  memcpy(shell_factor, ctx->ao_basis.shell_factor,
         (size_t) ctx->ao_basis.shell_num * sizeof(double));
  return QMCKL_SUCCESS;
}


    #+end_src


    #+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_get_ao_basis_exponent (const qmckl_context context,
                             double*  const exponent,
                             const int64_t size_max);
    #+end_src

    #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_ao_basis_exponent (const qmckl_context context,
                             double* const exponent,
                             const int64_t size_max)
{
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_CONTEXT,
                           "qmckl_get_ao_basis_exponent",
                           NULL);
  }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

  int32_t mask = 1 << 9;

  if ( (ctx->ao_basis.uninitialized & mask) != 0) {
    return qmckl_failwith( context,
                           QMCKL_NOT_PROVIDED,
                           "qmckl_get_ao_basis_exponent",
                           NULL);
  }

  if (exponent == NULL) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_get_ao_basis_exponent",
                           "NULL pointer");
  }

  if (size_max < ctx->ao_basis.prim_num) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_3,
                           "qmckl_get_ao_basis_exponent",
                           "Array too small. Expected prim_num");
  }

  assert (ctx->ao_basis.exponent != NULL);
  memcpy(exponent, ctx->ao_basis.exponent,
         (size_t) ctx->ao_basis.prim_num * sizeof(double));
  return QMCKL_SUCCESS;
}
    #+end_src


    #+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_get_ao_basis_coefficient (const qmckl_context context,
                                double*  const coefficient,
                                const int64_t size_max);
    #+end_src

    #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_ao_basis_coefficient (const qmckl_context context,
                                double* const coefficient,
                                const int64_t size_max)
{
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_CONTEXT,
                           "qmckl_get_ao_basis_coefficient",
                           NULL);
  }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

  int32_t mask = 1 << 10;

  if ( (ctx->ao_basis.uninitialized & mask) != 0) {
    return qmckl_failwith( context,
                           QMCKL_NOT_PROVIDED,
                           "qmckl_get_ao_basis_coefficient",
                           NULL);
  }

  if (coefficient == NULL) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_get_ao_basis_coefficient",
                           "NULL pointer");
  }

  if (size_max < ctx->ao_basis.prim_num) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_3,
                           "qmckl_get_ao_basis_coefficient",
                           "Array too small. Expected prim_num");
  }
  assert (ctx->ao_basis.coefficient != NULL);
  memcpy(coefficient, ctx->ao_basis.coefficient,
         (size_t) ctx->ao_basis.prim_num * sizeof(double));
  return QMCKL_SUCCESS;
}


    #+end_src


    #+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_get_ao_basis_prim_factor (const qmckl_context context,
                                double*  const prim_factor,
                                const int64_t size_max);
    #+end_src

    #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_ao_basis_prim_factor (const qmckl_context context,
                                double* const prim_factor,
                                const int64_t size_max)
{
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_CONTEXT,
                           "qmckl_get_ao_basis_prim_factor",
                           NULL);
  }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

  int32_t mask = 1 << 11;

  if ( (ctx->ao_basis.uninitialized & mask) != 0) {
    return qmckl_failwith( context,
                           QMCKL_NOT_PROVIDED,
                           "qmckl_get_ao_basis_prim_factor",
                           NULL);
  }

  if (prim_factor == NULL) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_get_ao_basis_prim_factor",
                           "NULL pointer");
  }

  if (size_max < ctx->ao_basis.prim_num) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_3,
                           "qmckl_get_ao_basis_prim_factor",
                           "Array too small. Expected prim_num");
  }

  assert (ctx->ao_basis.prim_factor != NULL);
  memcpy(prim_factor, ctx->ao_basis.prim_factor,
         (size_t) ctx->ao_basis.prim_num * sizeof(double));
  return QMCKL_SUCCESS;
}

    #+end_src


    #+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_get_ao_basis_ao_num (const qmckl_context context,
                           int64_t* const ao_num);
    #+end_src

    #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_ao_basis_ao_num (const qmckl_context context,
                           int64_t* const ao_num)
{
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_CONTEXT,
                           "qmckl_get_ao_basis_ao_num",
                           NULL);
  }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

  int32_t mask = 1 << 12;

  if ( (ctx->ao_basis.uninitialized & mask) != 0) {
    return qmckl_failwith( context,
                           QMCKL_NOT_PROVIDED,
                           "qmckl_get_ao_basis_ao_num",
                           NULL);
  }

  assert (ctx->ao_basis.ao_num > (int64_t) 0);

  *ao_num = ctx->ao_basis.ao_num;
  return QMCKL_SUCCESS;
}


    #+end_src


    #+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_get_ao_basis_ao_factor (const qmckl_context context,
                              double* const ao_factor,
                              const int64_t size_max);
    #+end_src

    #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_ao_basis_ao_factor (const qmckl_context context,
                              double* const ao_factor,
                              const int64_t size_max)
{
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_CONTEXT,
                           "qmckl_get_ao_basis_ao_factor",
                           NULL);
  }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

  int32_t mask = 1 << 13;

  if ( (ctx->ao_basis.uninitialized & mask) != 0) {
    return qmckl_failwith( context,
                           QMCKL_NOT_PROVIDED,
                           "qmckl_get_ao_basis_ao_factor",
                           NULL);
  }

  if (ao_factor == NULL) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_get_ao_basis_ao_factor",
                           "NULL pointer");
  }

  if (size_max < ctx->ao_basis.ao_num) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_3,
                           "qmckl_get_ao_basis_ao_factor",
                           "Array too small. Expected ao_num");
  }

  assert (ctx->ao_basis.ao_factor != NULL);
  memcpy(ao_factor, ctx->ao_basis.ao_factor, ctx->ao_basis.ao_num * sizeof(double));
  return QMCKL_SUCCESS;
}
    #+end_src


    When all the data for the AOs have been provided, the following
    function returns ~true~.

    #+begin_src c :comments org :tangle (eval h_func)
bool qmckl_ao_basis_provided (const qmckl_context context);
    #+end_src

    #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
bool qmckl_ao_basis_provided(const qmckl_context context) {

  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return false;
  }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

  return ctx->ao_basis.provided;
}
    #+end_src

**** Fortran interface

      #+begin_src f90 :tangle (eval fh_func) :comments org
interface
  integer(c_int32_t) function qmckl_get_ao_basis_type (context, &
       basis_type) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    character(c_char)   , intent(out)         :: basis_type
  end function qmckl_get_ao_basis_type
end interface

interface
  integer(c_int32_t) function qmckl_get_ao_basis_shell_num(context, &
       num) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    integer (c_int64_t) , intent(out)         :: num
  end function qmckl_get_ao_basis_shell_num
end interface

interface
  integer(c_int32_t) function qmckl_get_ao_basis_prim_num(context, &
       num) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    integer (c_int64_t) , intent(out)         :: num
  end function qmckl_get_ao_basis_prim_num
end interface

interface
  integer(c_int32_t) function qmckl_get_ao_basis_nucleus_shell_num(context, &
       shell_num, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    integer (c_int64_t) , intent(out)         :: shell_num(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_get_ao_basis_nucleus_shell_num
end interface

interface
  integer(c_int32_t) function qmckl_get_ao_basis_nucleus_index(context, &
       idx, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    integer (c_int64_t) , intent(out)         :: idx(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_get_ao_basis_nucleus_index
end interface

interface
  integer(c_int32_t) function qmckl_get_ao_basis_shell_ang_mom(context, &
       shell_ang_mom, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    integer (c_int32_t) , intent(out)         :: shell_ang_mom(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_get_ao_basis_shell_ang_mom
end interface

interface
  integer(c_int32_t) function qmckl_get_ao_basis_shell_prim_num(context, &
       shell_prim_num, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    integer (c_int64_t) , intent(out)         :: shell_prim_num(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_get_ao_basis_shell_prim_num
end interface

interface
  integer(c_int32_t) function qmckl_get_ao_basis_shell_prim_index(context, &
       shell_prim_index, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    integer (c_int64_t) , intent(out)         :: shell_prim_index(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_get_ao_basis_shell_prim_index
end interface

interface
  integer(c_int32_t) function qmckl_get_ao_basis_shell_factor(context, &
       shell_factor, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    real    (c_double)  , intent(out)         :: shell_factor(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_get_ao_basis_shell_factor
end interface

interface
  integer(c_int32_t) function qmckl_get_ao_basis_exponent(context, &
       exponent, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    real    (c_double)  , intent(out)         :: exponent(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_get_ao_basis_exponent
end interface

interface
  integer(c_int32_t) function qmckl_get_ao_basis_coefficient(context, &
       coefficient, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    real    (c_double)  , intent(out)         :: coefficient(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_get_ao_basis_coefficient
end interface

interface
  integer(c_int32_t) function qmckl_get_ao_basis_prim_factor(context, &
       prim_factor, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    real    (c_double)  , intent(out)         :: prim_factor(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_get_ao_basis_prim_factor
end interface

interface
  integer(c_int32_t) function qmckl_get_ao_basis_ao_num(context, &
       num) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    integer (c_int64_t) , intent(out)         :: num
  end function qmckl_get_ao_basis_ao_num
end interface

interface
  integer(c_int32_t) function qmckl_get_ao_basis_cartesian(context, &
       cartesian) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    logical (c_bool)    , intent(out)         :: cartesian
  end function qmckl_get_ao_basis_cartesian
end interface

interface
  integer(c_int32_t) function qmckl_get_ao_basis_ao_factor(context, &
       ao_factor, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    real    (c_double)  , intent(out)         :: ao_factor(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_get_ao_basis_ao_factor
end interface
     #+end_src

*** Test                                                           :noexport:

    #+begin_src c :tangle (eval c_test) :exports none :exports none
const int64_t   nucl_num      = chbrclf_nucl_num;
const double*   nucl_charge   = chbrclf_charge;
const double*   nucl_coord    = &(chbrclf_nucl_coord[0][0]);

qmckl_exit_code rc;
rc = qmckl_set_nucleus_num (context, nucl_num);
assert(rc == QMCKL_SUCCESS);

rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), 3*nucl_num);
assert(rc == QMCKL_SUCCESS);

rc = qmckl_set_nucleus_charge(context, nucl_charge, nucl_num);
assert(rc == QMCKL_SUCCESS);

assert(qmckl_nucleus_provided(context));


const int64_t    shell_num         =  chbrclf_shell_num;
const int64_t    prim_num          =  chbrclf_prim_num;
const int64_t    ao_num            =  chbrclf_ao_num;
const int64_t *  nucleus_index     =  &(chbrclf_basis_nucleus_index[0]);
const int64_t *  nucleus_shell_num =  &(chbrclf_basis_nucleus_shell_num[0]);
const int32_t *  shell_ang_mom     =  &(chbrclf_basis_shell_ang_mom[0]);
const int64_t *  shell_prim_num    =  &(chbrclf_basis_shell_prim_num[0]);
const int64_t *  shell_prim_index  =  &(chbrclf_basis_shell_prim_index[0]);
const double  *  shell_factor      =  &(chbrclf_basis_shell_factor[0]);
const double  *  exponent          =  &(chbrclf_basis_exponent[0]);
const double  *  coefficient       =  &(chbrclf_basis_coefficient[0]);
const double  *  prim_factor       =  &(chbrclf_basis_prim_factor[0]);
const double  *  ao_factor         =  &(chbrclf_basis_ao_factor[0]);

const char typ = 'G';

assert(!qmckl_ao_basis_provided(context));

rc = qmckl_set_ao_basis_type (context, typ);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));

rc = qmckl_set_ao_basis_shell_num (context, shell_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));

rc = qmckl_set_ao_basis_prim_num (context, prim_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));

rc = qmckl_set_ao_basis_nucleus_index (context, nucleus_index, nucl_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));

rc = qmckl_set_ao_basis_nucleus_shell_num (context, nucleus_shell_num, nucl_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));

rc = qmckl_set_ao_basis_shell_ang_mom (context, shell_ang_mom, shell_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));

rc = qmckl_set_ao_basis_shell_factor  (context, shell_factor, shell_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));

rc = qmckl_set_ao_basis_shell_prim_num (context, shell_prim_num, shell_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));

rc = qmckl_set_ao_basis_shell_prim_index (context, shell_prim_index, shell_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));

rc = qmckl_set_ao_basis_exponent      (context, exponent, prim_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));

rc = qmckl_set_ao_basis_coefficient   (context, coefficient, prim_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));

rc = qmckl_set_ao_basis_prim_factor (context, prim_factor, prim_num);
assert(rc == QMCKL_SUCCESS);

rc = qmckl_set_ao_basis_ao_num(context, ao_num);
assert(rc == QMCKL_SUCCESS);

rc = qmckl_set_ao_basis_ao_factor (context, ao_factor, ao_num);
assert(rc == QMCKL_SUCCESS);

assert(qmckl_ao_basis_provided(context));

int64_t    shell_num_test        ;
int64_t    prim_num_test         ;
int64_t    ao_num_test           ;
int64_t *  nucleus_index_test    ;
int64_t *  nucleus_shell_num_test;
int32_t *  shell_ang_mom_test    ;
int64_t *  shell_prim_num_test   ;
int64_t *  shell_prim_index_test ;
double  *  shell_factor_test     ;
double  *  exponent_test         ;
double  *  coefficient_test      ;
double  *  prim_factor_test      ;
double  *  ao_factor_test        ;
char    typ_test ;


rc = qmckl_get_ao_basis_type (context, &typ_test);
assert (rc == QMCKL_SUCCESS);
assert(typ == typ_test);

rc = qmckl_get_ao_basis_shell_num (context, &shell_num_test);
assert (rc == QMCKL_SUCCESS);
assert(shell_num == shell_num_test);

rc = qmckl_get_ao_basis_prim_num (context, &prim_num_test);
assert (rc == QMCKL_SUCCESS);
assert(prim_num == prim_num_test);

nucleus_index_test = (int64_t*) malloc (nucl_num * sizeof(int64_t));
rc = qmckl_get_ao_basis_nucleus_index (context, nucleus_index_test, nucl_num);
assert (rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i < nucl_num ; ++i) {
  assert(nucleus_index_test[i] == nucleus_index[i]);
 }
free(nucleus_index_test);

nucleus_shell_num_test = (int64_t*) malloc ( nucl_num * sizeof(int64_t));
rc = qmckl_get_ao_basis_nucleus_shell_num (context, nucleus_shell_num_test, nucl_num);
assert (rc == QMCKL_SUCCESS);

for (int64_t i=0 ; i < nucl_num ; ++i) {
  assert(nucleus_shell_num_test[i] == nucleus_shell_num[i]);
 }

shell_ang_mom_test = (int32_t*) malloc ( shell_num * sizeof(int32_t));
rc = qmckl_get_ao_basis_shell_ang_mom (context, shell_ang_mom_test, shell_num);
assert (rc == QMCKL_SUCCESS);

for (int64_t i=0 ; i < shell_num ; ++i) {
  assert(shell_ang_mom_test[i] == shell_ang_mom[i]);
 }
free(shell_ang_mom_test);

shell_factor_test = (double*) malloc ( shell_num * sizeof(double));
rc = qmckl_get_ao_basis_shell_factor  (context, shell_factor_test, shell_num);
assert (rc == QMCKL_SUCCESS);

for (int64_t i=0 ; i < shell_num ; ++i) {
  assert(shell_factor_test[i] == shell_factor[i]);
}
free(shell_factor_test);

shell_prim_num_test = (int64_t*) malloc ( shell_num * sizeof(int64_t));
rc = qmckl_get_ao_basis_shell_prim_num (context, shell_prim_num_test, shell_num);
assert (rc == QMCKL_SUCCESS);

for (int64_t i=0 ; i < shell_num ; ++i) {
  assert(shell_prim_num_test[i] == shell_prim_num[i]);
}
free(shell_prim_num_test);

shell_prim_index_test = (int64_t*) malloc ( shell_num * sizeof(int64_t));
rc = qmckl_get_ao_basis_shell_prim_index (context, shell_prim_index_test, shell_num);
assert (rc == QMCKL_SUCCESS);

for (int64_t i=0 ; i < shell_num ; ++i) {
  assert(shell_prim_index_test[i] == shell_prim_index[i]);
}
free(shell_prim_index_test);

exponent_test = (double*) malloc ( prim_num * sizeof(double));
rc = qmckl_get_ao_basis_exponent(context, exponent_test, prim_num);
assert (rc == QMCKL_SUCCESS);

for (int64_t i=0 ; i < prim_num ; ++i) {
  assert(exponent_test[i] == exponent[i]);
}
free(exponent_test);

coefficient_test = (double*) malloc ( prim_num * sizeof(double));
rc = qmckl_get_ao_basis_coefficient(context, coefficient_test, prim_num);
assert (rc == QMCKL_SUCCESS);

for (int64_t i=0 ; i < prim_num ; ++i) {
  assert(coefficient_test[i] == coefficient[i]);
}
free(coefficient_test);

prim_factor_test = (double*) malloc ( prim_num * sizeof(double));
rc = qmckl_get_ao_basis_prim_factor (context, prim_factor_test, prim_num);
assert (rc == QMCKL_SUCCESS);

for (int64_t i=0 ; i < prim_num ; ++i) {
  assert(prim_factor_test[i] == prim_factor[i]);
}
free(prim_factor_test);

rc = qmckl_get_ao_basis_ao_num(context, &ao_num_test);
assert(ao_num == ao_num_test);

ao_factor_test = (double*) malloc ( ao_num * sizeof(double));
rc = qmckl_get_ao_basis_ao_factor (context, ao_factor_test, ao_num);
assert (rc == QMCKL_SUCCESS);

for (int64_t i=0 ; i < ao_num ; ++i) {
  assert(ao_factor_test[i] == ao_factor[i]);
}
free(ao_factor_test);

    #+end_src

** Computed data

   The following data is computed as described in the next sections:

   |----------------------+-----------------------------------+----------------------------------------------------------------------------------------------|
   | Variable             | Type                              | Description                                                                                  |
   |----------------------+-----------------------------------+----------------------------------------------------------------------------------------------|
   | ~primitive_vgl~      | ~double[point_num][5][prim_num]~  | Value, gradients, Laplacian of the primitives at current  positions                          |
   | ~primitive_vgl_date~ | ~uint64_t~                        | Last modification date of Value, gradients, Laplacian of the primitives at current positions |
   | ~shell_vgl~          | ~double[point_num][5][shell_num]~ | Value, gradients, Laplacian of the primitives at current  positions                          |
   | ~shell_vgl_date~     | ~uint64_t~                        | Last modification date of Value, gradients, Laplacian of the AOs at current positions        |
   | ~ao_vgl~             | ~double[point_num][5][ao_num]~    | Value, gradients, Laplacian of the primitives at current  positions                          |
   | ~ao_vgl_date~        | ~uint64_t~                        | Last modification date of Value, gradients, Laplacian of the AOs at current positions        |

*** After initialization

    When the basis set is completely entered, extra data structures may be
    computed to accelerate the calculations. The primitives within each
    contraction are sorted in ascending order of their exponents, such
    that as soon as a primitive is zero all the following functions
    vanish. Also, it is possible to compute a nuclear radius beyond which
    all the primitives are zero up to the numerical accuracy defined in
    the context.

      #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_finalize_basis (qmckl_context context);
      #+end_src

      #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_finalize_basis(qmckl_context context) {

  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_CONTEXT,
                           "qmckl_finalize_basis",
                           NULL);
  }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

  int64_t nucl_num = 0;

  qmckl_exit_code rc = qmckl_get_nucleus_num(context, &nucl_num);
  if (rc != QMCKL_SUCCESS) return rc;

  /* nucleus_prim_index */
  {
    qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
    mem_info.size = (ctx->nucleus.num + (int64_t) 1) * sizeof(int64_t);

    ctx->ao_basis.nucleus_prim_index = (int64_t*) qmckl_malloc(context, mem_info);

    if (ctx->ao_basis.nucleus_prim_index == NULL) {
      return qmckl_failwith( context,
                             QMCKL_ALLOCATION_FAILED,
                             "ao_basis.nucleus_prim_index",
                             NULL);
    }

    for (int64_t i=0 ; i<nucl_num ; ++i) {
      int64_t shell_idx = ctx->ao_basis.nucleus_index[i];
      ctx->ao_basis.nucleus_prim_index[i] = ctx->ao_basis.shell_prim_index[shell_idx];
    }
    ctx->ao_basis.nucleus_prim_index[nucl_num] = ctx->ao_basis.prim_num;
  }


  /* Normalize coefficients */
  {
    qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
    mem_info.size = ctx->ao_basis.prim_num * sizeof(double);

    ctx->ao_basis.coefficient_normalized = (double *) qmckl_malloc(context, mem_info);

    if (ctx->ao_basis.coefficient_normalized == NULL) {
      return qmckl_failwith( context,
                             QMCKL_ALLOCATION_FAILED,
                             "ao_basis.coefficient_normalized",
                             NULL);
    }

    for (int64_t ishell=0 ; ishell < ctx->ao_basis.shell_num ; ++ishell) {
      for (int64_t iprim=ctx->ao_basis.shell_prim_index[ishell] ;
           iprim < ctx->ao_basis.shell_prim_index[ishell]+ctx->ao_basis.shell_prim_num[ishell] ;
           ++iprim) {
        ctx->ao_basis.coefficient_normalized[iprim] =
          ctx->ao_basis.coefficient[iprim] * ctx->ao_basis.prim_factor[iprim] *
          ctx->ao_basis.shell_factor[ishell];
      }
    }
  }


  /* Find max angular momentum on each nucleus */
  {
    qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
    mem_info.size = ctx->nucleus.num * sizeof(int32_t);

    ctx->ao_basis.nucleus_max_ang_mom = (int32_t *) qmckl_malloc(context, mem_info);

    if (ctx->ao_basis.nucleus_max_ang_mom == NULL) {
      return qmckl_failwith( context,
                             QMCKL_ALLOCATION_FAILED,
                             "ao_basis.nucleus_max_ang_mom",
                             NULL);
    }

    for (int64_t inucl=0 ; inucl < nucl_num ; ++inucl) {
      ctx->ao_basis.nucleus_max_ang_mom[inucl] = 0;
      for (int64_t ishell=ctx->ao_basis.nucleus_index[inucl]  ;
           ishell < ctx->ao_basis.nucleus_index[inucl] + ctx->ao_basis.nucleus_shell_num[inucl] ;
           ++ishell) {
        ctx->ao_basis.nucleus_max_ang_mom[inucl] =
          ctx->ao_basis.nucleus_max_ang_mom[inucl] > ctx->ao_basis.shell_ang_mom[ishell] ?
          ctx->ao_basis.nucleus_max_ang_mom[inucl] : ctx->ao_basis.shell_ang_mom[ishell] ;
      }
    }
  }

  /* Find distance beyond which all AOs are zero.
     The distance is obtained by sqrt(log(cutoff)*range) */
  {
    if (ctx->ao_basis.type == 'G') {
      qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
      mem_info.size = ctx->nucleus.num * sizeof(double);

      ctx->ao_basis.nucleus_range = (double *) qmckl_malloc(context, mem_info);

      if (ctx->ao_basis.nucleus_range == NULL) {
        return qmckl_failwith( context,
                               QMCKL_ALLOCATION_FAILED,
                               "ao_basis.nucleus_range",
                               NULL);
      }

      for (int64_t inucl=0 ; inucl < ctx->nucleus.num ; ++inucl) {
        ctx->ao_basis.nucleus_range[inucl] = 0.;
        for (int64_t ishell=ctx->ao_basis.nucleus_index[inucl]  ;
             ishell < ctx->ao_basis.nucleus_index[inucl] + ctx->ao_basis.nucleus_shell_num[inucl] ;
             ++ishell) {
          for (int64_t iprim=ctx->ao_basis.shell_prim_index[ishell] ;
               iprim < ctx->ao_basis.shell_prim_index[ishell] + ctx->ao_basis.shell_prim_num[ishell] ;
               ++iprim) {
            double range = 1./ctx->ao_basis.exponent[iprim];
            ctx->ao_basis.nucleus_range[inucl] =
              ctx->ao_basis.nucleus_range[inucl] > range ?
              ctx->ao_basis.nucleus_range[inucl] : range;
          }
        }
      }
    }
  }

  rc = QMCKL_SUCCESS;
#ifdef HAVE_HPC
  rc = qmckl_finalize_basis_hpc(context);
#endif

  return rc;
}
      #+end_src

*** HPC-specific data structures

    For faster access, we provide extra arrays for the shell information as:
    - $C_{psa}$ = =coef_per_nucleus (iprim, ishell, inucl)=
    - $\gamma_{pa}$ =expo_per_nucleus (iprim, inucl)=

    such that the computation of the radial parts can be vectorized.
    
    Exponents are sorted in increasing order to exit loops quickly,
    and only unique exponents are kept. This also allows to pack the
    exponents to enable vectorization of exponentials.

    The computation of the radial part is made as
    \[
      R_{sa} = \sum_p C_{psa} \gamma_{pa}
    \]
    which is a matrix-vector product.
    
    #+NAME:HPC_struct 
    #+begin_src c :comments org :exports none
 /* HPC specific data structures */
 qmckl_tensor coef_per_nucleus;
 qmckl_matrix expo_per_nucleus;
    #+end_src

    #+begin_src c :comments org :tangle (eval h_private_func) :exports none
qmckl_exit_code qmckl_finalize_basis_hpc (qmckl_context context);
    #+end_src
    
    #+begin_src c :comments org :tangle (eval c) :exports none
qmckl_exit_code qmckl_finalize_basis_hpc (qmckl_context context)
{
  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;

  int64_t shell_max = 0;
  int64_t prim_max  = 0;
  int64_t nucl_num = ctx->nucleus.num;
  for (int inucl=0 ; inucl < nucl_num ; ++inucl) {
    shell_max = ctx->ao_basis.nucleus_shell_num[inucl] > shell_max ?
                ctx->ao_basis.nucleus_shell_num[inucl] : shell_max;

    int64_t prim_num = 0;
    const int64_t ishell_start = ctx->ao_basis.nucleus_index[inucl];
    const int64_t ishell_end   = ctx->ao_basis.nucleus_index[inucl] + ctx->ao_basis.nucleus_shell_num[inucl];
    for (int64_t ishell = ishell_start ; ishell < ishell_end ; ++ishell) {
      prim_num += ctx->ao_basis.shell_prim_num[ishell];
    }
    prim_max = prim_num > prim_max ?
               prim_num : prim_max;
  }

  int64_t size[3] = { prim_max, shell_max, nucl_num };
  ctx->ao_basis.coef_per_nucleus = qmckl_tensor_alloc( context, 3, size );
  ctx->ao_basis.coef_per_nucleus = qmckl_tensor_set(ctx->ao_basis.coef_per_nucleus, 0.);

  ctx->ao_basis.expo_per_nucleus = qmckl_matrix_alloc( context, prim_max, nucl_num );
  ctx->ao_basis.expo_per_nucleus = qmckl_matrix_set(ctx->ao_basis.expo_per_nucleus, 0.);

  double expo[prim_max];
  double coef[shell_max][prim_max];
  for (int inucl=0 ; inucl < nucl_num ; ++inucl) {
    memset(expo, 0, sizeof(expo));
    memset(coef, 0, sizeof(expo));

    int64_t idx = 0;
    const int64_t ishell_start = ctx->ao_basis.nucleus_index[inucl];
    const int64_t ishell_end   = ctx->ao_basis.nucleus_index[inucl] + ctx->ao_basis.nucleus_shell_num[inucl];
    for (int64_t ishell = ishell_start ; ishell < ishell_end ; ++ishell) {
      
      const int64_t iprim_start = ctx->ao_basis.shell_prim_index[ishell];
      const int64_t iprim_end   = ctx->ao_basis.shell_prim_index[ishell] + ctx->ao_basis.shell_prim_num[ishell];
      for (int64_t iprim = iprim_start ; iprim < iprim_end ; ++iprim) {
        expo[idx] = ctx->ao_basis.coefficient_normalized[iprim];
        coef[ishell - ishell_start][idx] = ctx->ao_basis.coefficient_normalized[iprim];
        idx += 1;
      }
    }

    for (int64_t i=0 ; i<prim_max ; ++i) {
      qmckl_mat( ctx->ao_basis.expo_per_nucleus, i, inucl ) = expo[i];
    }

    for (int64_t j=0 ; j<shell_max ; ++j) {
      for (int64_t i=0 ; i<prim_max ; ++i) {
        qmckl_ten3( ctx->ao_basis.coef_per_nucleus, i, j, inucl ) = coef[j][i];
      }
    }

  }
  return QMCKL_SUCCESS;
}
  
    #+end_src

*** Access functions

    #+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code
qmckl_get_ao_basis_primitive_vgl (qmckl_context context,
                                  double* const primitive_vgl,
                                  const int64_t size_max);
    #+end_src

    Returns the array of values, gradients an Laplacian of primitive
    basis functions evaluated at the current coordinates.
    See section [[Computation of primitives]].

    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
qmckl_exit_code
qmckl_get_ao_basis_primitive_vgl (qmckl_context context,
                                  double* const primitive_vgl,
                                  const int64_t size_max)
{
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_CONTEXT,
                           "qmckl_get_ao_basis_primitive_vgl",
                           NULL);
  }

  if (size_max <= 0) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_3,
                           "qmckl_get_ao_basis_primitive_vgl",
                           "size_max <= 0");
  }

  qmckl_exit_code rc;

  rc = qmckl_provide_ao_basis_primitive_vgl(context);
  if (rc != QMCKL_SUCCESS) return rc;

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

  int64_t sze = ctx->ao_basis.prim_num * 5 * ctx->point.num;
  if (size_max < sze) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_3,
                           "qmckl_get_ao_basis_primitive_vgl",
                           "input array too small");
  }
  memcpy(primitive_vgl, ctx->ao_basis.primitive_vgl, (size_t) sze * sizeof(double));

  return QMCKL_SUCCESS;
}
    #+end_src

    #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
    interface
      integer(c_int32_t) function qmckl_get_ao_basis_primitive_vgl &
           (context, primitive_vgl, size_max) &
           bind(C)
        use, intrinsic :: iso_c_binding
        import
        implicit none

        integer (c_int64_t) , intent(in)  , value :: context
        double precision,     intent(out)         :: primitive_vgl(*)
        integer (c_int64_t) , intent(in)  , value :: size_max
      end function
    end interface
    #+end_src


    #+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code
qmckl_get_ao_basis_shell_vgl (qmckl_context context,
                              double* const shell_vgl,
                              const int64_t size_max);
    #+end_src

    Returns the array of values, gradients an Laplacian of contracted shells
    evaluated at the current coordinates. See section [[Computation of shells]].

    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
qmckl_exit_code
qmckl_get_ao_basis_shell_vgl (qmckl_context context,
                              double* const shell_vgl,
                              const int64_t size_max)
{
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_CONTEXT,
                           "qmckl_get_ao_basis_shell_vgl",
                           NULL);
  }

  qmckl_exit_code rc;

  rc = qmckl_provide_ao_basis_shell_vgl(context);
  if (rc != QMCKL_SUCCESS) return rc;

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

  int64_t sze = ctx->ao_basis.shell_num * 5 * ctx->point.num;
  if (size_max < sze) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_3,
                           "qmckl_get_ao_basis_shell_vgl",
                           "input array too small");
  }
  memcpy(shell_vgl, ctx->ao_basis.shell_vgl, (size_t)sze * sizeof(double));

  return QMCKL_SUCCESS;
}
    #+end_src

    #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
    interface
      integer(c_int32_t) function qmckl_get_ao_basis_shell_vgl &
           (context, shell_vgl, size_max) &
           bind(C)
        use, intrinsic :: iso_c_binding
        import
        implicit none

        integer (c_int64_t) , intent(in)  , value :: context
        double precision,     intent(out)         :: shell_vgl(*)
        integer (c_int64_t) , intent(in)  , value :: size_max
      end function
    end interface
    #+end_src


    #+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code
qmckl_get_ao_basis_ao_vgl (qmckl_context context,
                           double* const ao_vgl,
                           const int64_t size_max);
    #+end_src

    Returns the array of values, gradients an Laplacian of the atomic orbitals
    evaluated at the current coordinates.
    See section [[Combining radial and polynomial parts]].

    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
qmckl_exit_code
qmckl_get_ao_basis_ao_vgl (qmckl_context context,
                           double* const ao_vgl,
                           const int64_t size_max)
{

  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_CONTEXT,
                           "qmckl_get_ao_basis_ao_vgl",
                           NULL);
  }

  qmckl_exit_code rc;

  rc = qmckl_provide_ao_vgl(context);
  if (rc != QMCKL_SUCCESS) return rc;

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

  int64_t sze = ctx->ao_basis.ao_num * 5 * ctx->point.num;
  if (size_max < sze) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_3,
                           "qmckl_get_ao_basis_ao_vgl",
                           "input array too small");
  }
  memcpy(ao_vgl, ctx->ao_basis.ao_vgl, (size_t) sze * sizeof(double));

  return QMCKL_SUCCESS;
}
    #+end_src

    #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
  interface
     integer(c_int32_t) function qmckl_get_ao_basis_ao_vgl (context, &
          ao_vgl, size_max) bind(C)
       use, intrinsic :: iso_c_binding
       import
       implicit none
       integer (c_int64_t) , intent(in)  , value :: context
       double precision,     intent(out)         :: ao_vgl(*)
       integer (c_int64_t) , intent(in)  , value :: size_max
     end function qmckl_get_ao_basis_ao_vgl
  end interface
    #+end_src

    Uses the give array to compute the VGL.

    #+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code
qmckl_get_ao_basis_ao_vgl_inplace (qmckl_context context,
                                   double* const ao_vgl,
                                   const int64_t size_max);
    #+end_src

    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
qmckl_exit_code
qmckl_get_ao_basis_ao_vgl_inplace (qmckl_context context,
                                   double* const ao_vgl,
                                   const int64_t size_max)
{

  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_CONTEXT,
                           "qmckl_get_ao_basis_ao_vgl",
                           NULL);
  }

  qmckl_exit_code rc;

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

  int64_t sze = ctx->ao_basis.ao_num * 5 * ctx->point.num;
  if (size_max < sze) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_3,
                           "qmckl_get_ao_basis_ao_vgl",
                           "input array too small");
  }

  rc = qmckl_context_touch(context);
  if (rc != QMCKL_SUCCESS) return rc;

  double* old_array = ctx->ao_basis.ao_vgl;

  ctx->ao_basis.ao_vgl = ao_vgl;

  rc = qmckl_provide_ao_vgl(context);
  if (rc != QMCKL_SUCCESS) return rc;

  ctx->ao_basis.ao_vgl = old_array;

  return QMCKL_SUCCESS;
}
    #+end_src

    #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
  interface
     integer(c_int32_t) function qmckl_get_ao_basis_ao_vgl_inplace (context, &
          ao_vgl, size_max) bind(C)
       use, intrinsic :: iso_c_binding
       import
       implicit none
       integer (c_int64_t) , intent(in)  , value :: context
       double precision,     intent(out)         :: ao_vgl(*)
       integer (c_int64_t) , intent(in)  , value :: size_max
     end function qmckl_get_ao_basis_ao_vgl_inplace
  end interface
    #+end_src

* Radial part

** General functions for 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 \]

   |--------------+------------------+------------------------------------------------------|
   | Variable     | Type             | Description                                          |
   |--------------+------------------+------------------------------------------------------|
   | ~context~    | ~qmckl_context~  | Global state                                         |
   | ~X(3)~       | ~double[3]~      | Array containing the coordinates of the points       |
   | ~R(3)~       | ~double[3]~      | Array containing the x,y,z coordinates of the center |
   | ~n~          | ~int64_t~        | Number of computed Gaussians                         |
   | ~A(n)~       | ~double[n]~      | Exponents of the Gaussians                           |
   | ~VGL(ldv,5)~ | ~double[5][ldv]~ | Value, gradients and Laplacian of the Gaussians      |
   | ~ldv~        | ~int64_t~        | Leading dimension of array ~VGL~                     |
   |--------------+------------------+------------------------------------------------------|

   Requirements:

        - ~context~ \ne 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

   #+begin_src c :tangle (eval h_func)
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

   #+begin_src f90 :tangle (eval f)
integer function qmckl_ao_gaussian_vgl_f(context, X, R, n, A, VGL, ldv) result(info)
  use qmckl
  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 = QMCKL_SUCCESS

  if (context == QMCKL_NULL_CONTEXT) then
     info = QMCKL_INVALID_CONTEXT
     return
  endif

  if (n <= 0) then
     info = QMCKL_INVALID_ARG_4
     return
  endif

  if (ldv < n) then
     info = QMCKL_INVALID_ARG_7
     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

   #+begin_src f90 :tangle (eval f) :exports none
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 (eval fh_func) :exports none
  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 (eval f_test)
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, z
  double precision, allocatable :: VGL(:,:), A(:)
  double precision              :: epsilon

  epsilon = 3.d0 * qmckl_get_numprec_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
     z = dabs(1.d0 - VGL(i,1) / (dexp(-A(i) * r2)) )
     if ( z > epsilon ) then
        print *, z, epsilon
        return
     end if

     test_qmckl_ao_gaussian_vgl = -12
     z = dabs(1.d0 - VGL(i,2) / (&
          -2.d0 * A(i) * Y(1) * dexp(-A(i) * r2) ))
     if ( z > epsilon ) then
        print *, z, epsilon
        return
     end if

     test_qmckl_ao_gaussian_vgl = -13
     z = dabs(1.d0 - VGL(i,3) / (&
          -2.d0 * A(i) * Y(2) * dexp(-A(i) * r2) ))
     if ( z > epsilon ) then
        print *, z, epsilon
        return
     end if

     test_qmckl_ao_gaussian_vgl = -14
     z = dabs(1.d0 - VGL(i,4) / (&
          -2.d0 * A(i) * Y(3) * dexp(-A(i) * r2) ))
     if ( z > epsilon ) then
        print *, z, epsilon
        return
     end if

     test_qmckl_ao_gaussian_vgl = -15
     z = dabs(1.d0 - VGL(i,5) / (&
          A(i) * (4.d0*r2*A(i) - 6.d0) * dexp(-A(i) * r2) ))
     if ( z > epsilon ) then
        print *, z, epsilon
        return
     end if
  end do

  test_qmckl_ao_gaussian_vgl = 0

  deallocate(VGL)
end function test_qmckl_ao_gaussian_vgl
    #+end_src

    #+begin_src c :tangle (eval c_test) :exports none
int test_qmckl_ao_gaussian_vgl(qmckl_context context);
assert(0 == test_qmckl_ao_gaussian_vgl(context));
    #+end_src

** TODO General functions for Slater basis functions               :noexport:
** TODO General functions for Radial functions on a grid           :noexport:
** TODO Helper functions to accelerate calculations                :noexport:

  |--------------------------+---------------------+--------------------------------------------------------|
  | Variable                 | Type                | Description                                            |
  |--------------------------+---------------------+--------------------------------------------------------|
  | ~coefficient_normalized~ | ~double[prim_num]~  | Normalized primitive coefficients                      |
  | ~nucleus_prim_index~     | ~int64_t[nucl_num]~ | Index of the first primitive for each nucleus          |
  | ~nucleus_max_ang_mom~    | ~int32_t[nucl_num]~ | Maximum angular momentum for each nucleus              |
  | ~nucleus_range~          | ~double[nucl_num]~  | Distance beyond which all the AOs are zero             |
  |--------------------------+---------------------+--------------------------------------------------------|
  | ~nucl_shell_index~       | ~int64_t[nucl_num]~ | Index of the first shell for each nucleus              |
  | ~exponent_sorted~        | ~double[prim_num]~  | Array of exponents for sorted primitives               |
  | ~coeff_norm_sorted~      | ~double[prim_num]~  | Array of normalized coefficients for sorted primitives |
  | ~prim_factor_sorted~     | ~double[prim_num]~  | Normalization factors of the sorted primtives          |

** Computation of primitives
   :PROPERTIES:
   :Name:     qmckl_compute_ao_basis_primitive_gaussian_vgl
   :CRetType: qmckl_exit_code
   :FRetType: qmckl_exit_code
   :END:

   #+NAME: qmckl_ao_basis_primitive_gaussian_vgl_args
   | Variable             | Type                             | In/Out | Description                                      |
   |----------------------+----------------------------------+--------+--------------------------------------------------|
   | ~context~            | ~qmckl_context~                  | in     | Global state                                     |
   | ~prim_num~           | ~int64_t~                        | in     | Number of primitives                             |
   | ~point_num~          | ~int64_t~                        | in     | Number of points considered                      |
   | ~nucl_num~           | ~int64_t~                        | in     | Number of nuclei                                 |
   | ~nucleus_prim_index~ | ~int64_t[nucl_num]~              | in     | Index of the 1st primitive of each nucleus       |
   | ~coord~              | ~double[3][point_num]~           | in     | Coordinates                                      |
   | ~nucl_coord~         | ~double[3][nucl_num]~            | in     | Nuclear coordinates                              |
   | ~expo~               | ~double[prim_num]~               | in     | Exponents of the primitives                      |
   | ~primitive_vgl~      | ~double[point_num][5][prim_num]~ | out    | Value, gradients and Laplacian of the primitives |

   #+CALL: generate_c_header(table=qmckl_ao_basis_primitive_gaussian_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_ao_basis_primitive_gaussian_vgl"))

   #+RESULTS:
   #+begin_src c :tangle (eval h_func) :comments org
   qmckl_exit_code qmckl_compute_ao_basis_primitive_gaussian_vgl (
         const qmckl_context context,
         const int64_t prim_num,
         const int64_t point_num,
         const int64_t nucl_num,
         const int64_t* nucleus_prim_index,
         const double* coord,
         const double* nucl_coord,
         const double* expo,
         double* const primitive_vgl );
   #+end_src


   #+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_ao_basis_primitive_gaussian_vgl_f( &
     context, prim_num, point_num, nucl_num,                       &
     nucleus_prim_index, coord, nucl_coord,                  &
     expo, primitive_vgl)                                         &
     result(info)

  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer*8             , intent(in)  :: prim_num
  integer*8             , intent(in)  :: nucl_num
  integer*8             , intent(in)  :: point_num
  integer*8             , intent(in)  :: nucleus_prim_index(nucl_num+1)
  double precision      , intent(in)  :: coord(point_num,3)
  double precision      , intent(in)  :: nucl_coord(nucl_num,3)
  double precision      , intent(in)  :: expo(prim_num)
  double precision      , intent(out) :: primitive_vgl(prim_num,5,point_num)

  integer*8 :: inucl, iprim, ipoint
  double precision :: x, y, z, two_a, ar2, r2, v, cutoff

  info = QMCKL_SUCCESS

  ! Don't compute exponentials when the result will be almost zero.
  cutoff = -dlog(1.d-15)

  do inucl=1,nucl_num
     ! C is zero-based, so shift bounds by one
     do iprim = nucleus_prim_index(inucl)+1, nucleus_prim_index(inucl+1)
        do ipoint = 1, point_num
           x = coord(ipoint,1) - nucl_coord(inucl,1)
           y = coord(ipoint,2) - nucl_coord(inucl,2)
           z = coord(ipoint,3) - nucl_coord(inucl,3)

           r2 = x*x + y*y + z*z
           ar2 = expo(iprim)*r2
           if (ar2 > cutoff) cycle

           v = dexp(-ar2)
           two_a = -2.d0 * expo(iprim) * v

           primitive_vgl(iprim, 1, ipoint) = v
           primitive_vgl(iprim, 2, ipoint) = two_a * x
           primitive_vgl(iprim, 3, ipoint) = two_a * y
           primitive_vgl(iprim, 4, ipoint) = two_a * z
           primitive_vgl(iprim, 5, ipoint) = two_a * (3.d0 - 2.d0*ar2)

        end do
     end do
  end do

end function qmckl_compute_ao_basis_primitive_gaussian_vgl_f
   #+end_src

   #+CALL: generate_c_interface(table=qmckl_ao_basis_primitive_gaussian_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_ao_basis_primitive_gaussian_vgl"))

   #+RESULTS:
   #+begin_src f90 :tangle (eval f) :comments org :exports none
   integer(c_int32_t) function qmckl_compute_ao_basis_primitive_gaussian_vgl &
       (context, &
        prim_num, &
        point_num, &
        nucl_num, &
        nucleus_prim_index, &
        coord, &
        nucl_coord, &
        expo, &
        primitive_vgl) &
       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 :: prim_num
     integer (c_int64_t) , intent(in)  , value :: point_num
     integer (c_int64_t) , intent(in)  , value :: nucl_num
     integer (c_int64_t) , intent(in)          :: nucleus_prim_index(nucl_num)
     real    (c_double ) , intent(in)          :: coord(point_num,3)
     real    (c_double ) , intent(in)          :: nucl_coord(nucl_num,3)
     real    (c_double ) , intent(in)          :: expo(prim_num)
     real    (c_double ) , intent(out)         :: primitive_vgl(prim_num,5,point_num)

     integer(c_int32_t), external :: qmckl_compute_ao_basis_primitive_gaussian_vgl_f
     info = qmckl_compute_ao_basis_primitive_gaussian_vgl_f &
            (context, &
        prim_num, &
        point_num, &
        nucl_num, &
        nucleus_prim_index, &
        coord, &
        nucl_coord, &
        expo, &
        primitive_vgl)

   end function qmckl_compute_ao_basis_primitive_gaussian_vgl
   #+end_src

*** Provide                                                        :noexport:

    #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_provide_ao_basis_primitive_vgl(qmckl_context context);
    #+end_src

    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
qmckl_exit_code qmckl_provide_ao_basis_primitive_vgl(qmckl_context context)
{

  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_CONTEXT,
                           "qmckl_provide_get_ao_basis_primitive_vgl",
                           NULL);
  }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

  if (!ctx->ao_basis.provided) {
    return qmckl_failwith( context,
                           QMCKL_NOT_PROVIDED,
                           "qmckl_ao_basis_primitive_vgl",
                           NULL);
  }

  /* Compute if necessary */
  if (ctx->point.date > ctx->ao_basis.primitive_vgl_date) {

    /* Allocate array */
    if (ctx->ao_basis.primitive_vgl == NULL) {

      qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
      mem_info.size = ctx->ao_basis.prim_num * 5 * ctx->point.num *
        sizeof(double);
      double* primitive_vgl = (double*) qmckl_malloc(context, mem_info);

      if (primitive_vgl == NULL) {
        return qmckl_failwith( context,
                               QMCKL_ALLOCATION_FAILED,
                               "qmckl_ao_basis_primitive_vgl",
                               NULL);
      }
      ctx->ao_basis.primitive_vgl = primitive_vgl;
    }

    qmckl_exit_code rc;
    if (ctx->ao_basis.type == 'G') {
      rc = qmckl_compute_ao_basis_primitive_gaussian_vgl(context,
                                                         ctx->ao_basis.prim_num,
                                                         ctx->point.num,
                                                         ctx->nucleus.num,
                                                         ctx->ao_basis.nucleus_prim_index,
                                                         ctx->point.coord.data,
                                                         ctx->nucleus.coord.data,
                                                         ctx->ao_basis.exponent,
                                                         ctx->ao_basis.primitive_vgl);
    } else {
      return qmckl_failwith( context,
                             QMCKL_FAILURE,
                             "compute_ao_basis_primitive_vgl",
                             "Not yet implemented");
    }
    if (rc != QMCKL_SUCCESS) {
      return rc;
    }

    ctx->ao_basis.primitive_vgl_date = ctx->date;
  }

  return QMCKL_SUCCESS;
}
    #+end_src

*** Test                                                           :noexport:

    #+begin_src python :results output :exports none :exports none
import numpy as np

def f(a,x,y):
    return np.exp( -a*(np.linalg.norm(x-y))**2 )

def df(a,x,y,n):
    h0 = 1.e-6
    if   n == 1: h = np.array([h0,0.,0.])
    elif n == 2: h = np.array([0.,h0,0.])
    elif n == 3: h = np.array([0.,0.,h0])
    return ( f(a,x+h,y) - f(a,x-h,y) ) / (2.*h0)

def d2f(a,x,y,n):
    h0 = 1.e-6
    if   n == 1: h = np.array([h0,0.,0.])
    elif n == 2: h = np.array([0.,h0,0.])
    elif n == 3: h = np.array([0.,0.,h0])
    return ( f(a,x+h,y) - 2.*f(a,x,y) + f(a,x-h,y) ) / h0**2

def lf(a,x,y):
    return d2f(a,x,y,1) + d2f(a,x,y,2) + d2f(a,x,y,3)

elec_26_w1 = np.array( [  1.49050402641, 2.90106987953, -1.05920815468  ] )
nucl_1    = np.array( [ 1.096243353458458e+00, 8.907054016973815e-01, 7.777092280258892e-01 ] )
nucl_2    = np.array( [ 1.168459237342663e+00, 1.125660720053393e+00, 2.833370314829343e+00 ] )

#double prim_vgl[prim_num][5][elec_num];
a =  0.9059; x = elec_26_w1 ; y = nucl_1
print ( "[7][0][26] : %e"% f(a,x,y))
print ( "[7][1][26] : %e"% df(a,x,y,1))
print ( "[7][2][26] : %e"% df(a,x,y,2))
print ( "[7][3][26] : %e"% df(a,x,y,3))
print ( "[7][4][26] : %e"% lf(a,x,y))

    #+end_src

    #+RESULTS:
    : [7][0][26] : 1.050157e-03
    : [7][1][26] : -7.501497e-04
    : [7][2][26] : -3.825069e-03
    : [7][3][26] : 3.495056e-03
    : [7][4][26] : 2.040013e-02

    #+begin_src c :tangle (eval c_test) :exports none
{
#define walk_num 1 // chbrclf_walk_num
#define elec_num chbrclf_elec_num
#define prim_num chbrclf_prim_num

  int64_t elec_up_num   = chbrclf_elec_up_num;
  int64_t elec_dn_num   = chbrclf_elec_dn_num;
  double* elec_coord    = &(chbrclf_elec_coord[0][0][0]);

  rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num);
  assert (rc == QMCKL_SUCCESS);

  rc = qmckl_set_electron_walk_num (context, walk_num);
  assert (rc == QMCKL_SUCCESS);

  assert(qmckl_electron_provided(context));

  rc = qmckl_set_electron_coord (context, 'N', elec_coord, walk_num*elec_num*3);
  assert(rc == QMCKL_SUCCESS);



  double prim_vgl[elec_num*walk_num][5][prim_num];

  rc = qmckl_get_ao_basis_primitive_vgl(context, &(prim_vgl[0][0][0]),
          (int64_t) 5*elec_num*walk_num*prim_num );
  assert (rc == QMCKL_SUCCESS);

  assert( fabs(prim_vgl[26][0][7] - ( 1.0501570432064878E-003)) < 1.e-14 );
  assert( fabs(prim_vgl[26][1][7] - (-7.5014974095310560E-004)) < 1.e-14 );
  assert( fabs(prim_vgl[26][2][7] - (-3.8250692897610380E-003)) < 1.e-14 );
  assert( fabs(prim_vgl[26][3][7] - ( 3.4950559194080275E-003)) < 1.e-14 );
  assert( fabs(prim_vgl[26][4][7] - ( 2.0392163767356572E-002)) < 1.e-14 );

}


    #+end_src

*** Ideas for improvement                                          :noexport:

    #+begin_src c
// j : electrons
// l : primitives

k=0;
for (j=0 ; j<point_num ; ++j) {
  for (i=0 ; i<nucl_num ; ++i) {

    r2 = nucl_point_dist[i][j];

    if (r2 < nucl_radius2[i]) {

      for (l=0 ; l<prim_num ; ++l) {
        tmp[k].i   = i;
        tmp[k].j   = j;
        tmp[k].m   = m;
        tmp[k].ar2 = -expo[l] *r2;
        ++k;
      }
    }
  }
 }

// sort(tmp) in increasing ar2;
// Identify first ar2 above numerical accuracy threshold
// Compute vectorized exponentials on significant values
    #+end_src

** Computation of shells
   :PROPERTIES:
   :Name:     qmckl_compute_ao_basis_shell_gaussian_vgl
   :CRetType: qmckl_exit_code
   :FRetType: qmckl_exit_code
   :END:

   #+NAME: qmckl_ao_basis_shell_gaussian_vgl_args
   | Variable            | Type                              | In/Out | Description                                  |
   |---------------------+-----------------------------------+--------+----------------------------------------------|
   | ~context~           | ~qmckl_context~                   | in     | Global state                                 |
   | ~prim_num~          | ~int64_t~                         | in     | Number of primitives                         |
   | ~shell_num~         | ~int64_t~                         | in     | Number of shells                             |
   | ~point_num~         | ~int64_t~                         | in     | Number of points                             |
   | ~nucl_num~          | ~int64_t~                         | in     | Number of nuclei                             |
   | ~nucleus_shell_num~ | ~int64_t[nucl_num]~               | in     | Number of shells for each nucleus            |
   | ~nucleus_index~     | ~int64_t[nucl_num]~               | in     | Index of the 1st shell of each nucleus       |
   | ~nucleus_range~     | ~double[nucl_num]~                | in     | Range of the nucleus                         |
   | ~shell_prim_index~  | ~int64_t[shell_num]~              | in     | Index of the 1st primitive of each shell     |
   | ~shell_prim_num~    | ~int64_t[shell_num]~              | in     | Number of primitives per shell               |
   | ~coord~             | ~double[3][point_num]~            | in     | Coordinates                                  |
   | ~nucl_coord~        | ~double[3][nucl_num]~             | in     | Nuclear coordinates                          |
   | ~expo~              | ~double[prim_num]~                | in     | Exponents of the primitives                  |
   | ~coef_normalized~   | ~double[prim_num]~                | in     | Coefficients of the primitives               |
   | ~shell_vgl~         | ~double[point_num][5][shell_num]~ | out    | Value, gradients and Laplacian of the shells |

   #+CALL: generate_c_header(table=qmckl_ao_basis_shell_gaussian_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_ao_basis_shell_gaussian_vgl"))

   #+RESULTS:
   #+begin_src c :tangle (eval h_func) :comments org
   qmckl_exit_code qmckl_compute_ao_basis_shell_gaussian_vgl (
         const qmckl_context context,
         const int64_t prim_num,
         const int64_t shell_num,
         const int64_t point_num,
         const int64_t nucl_num,
         const int64_t* nucleus_shell_num,
         const int64_t* nucleus_index,
         const double* nucleus_range,
         const int64_t* shell_prim_index,
         const int64_t* shell_prim_num,
         const double* coord,
         const double* nucl_coord,
         const double* expo,
         const double* coef_normalized,
         double* const shell_vgl );
   #+end_src

   #+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_ao_basis_shell_gaussian_vgl_f( &
     context, prim_num, shell_num, point_num, nucl_num,       &
     nucleus_shell_num, nucleus_index, nucleus_range,         &
     shell_prim_index, shell_prim_num, coord, nucl_coord,     &
     expo, coef_normalized, shell_vgl)                        &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer*8             , intent(in)  :: prim_num
  integer*8             , intent(in)  :: shell_num
  integer*8             , intent(in)  :: nucl_num
  integer*8             , intent(in)  :: point_num
  integer*8             , intent(in)  :: nucleus_shell_num(nucl_num)
  integer*8             , intent(in)  :: nucleus_index(nucl_num)
  double precision      , intent(in)  :: nucleus_range(nucl_num)
  integer*8             , intent(in)  :: shell_prim_index(shell_num)
  integer*8             , intent(in)  :: shell_prim_num(shell_num)
  double precision      , intent(in)  :: coord(point_num,3)
  double precision      , intent(in)  :: nucl_coord(nucl_num,3)
  double precision      , intent(in)  :: expo(prim_num)
  double precision      , intent(in)  :: coef_normalized(prim_num)
  double precision      , intent(out) :: shell_vgl(shell_num,5,point_num)

  integer*8 :: inucl, iprim, ipoint, ishell
  integer*8 :: ishell_start, ishell_end
  integer*8 :: iprim_start , iprim_end
  double precision :: x, y, z, two_a, ar2, r2, v, cutoff

  info = QMCKL_SUCCESS

  ! Don't compute exponentials when the result will be almost zero.
  ! TODO : Use numerical precision here
  cutoff = -dlog(1.d-12)

  do ipoint = 1, point_num

     do inucl=1,nucl_num

        x = coord(ipoint,1) - nucl_coord(inucl,1)
        y = coord(ipoint,2) - nucl_coord(inucl,2)
        z = coord(ipoint,3) - nucl_coord(inucl,3)

        r2 = x*x + y*y + z*z

        if (r2 > cutoff*nucleus_range(inucl)) then
           cycle
        end if

        ! C is zero-based, so shift bounds by one
        ishell_start = nucleus_index(inucl) + 1
        ishell_end   = nucleus_index(inucl) + nucleus_shell_num(inucl)

        do ishell=ishell_start, ishell_end

           shell_vgl(ishell, 1, ipoint) = 0.d0
           shell_vgl(ishell, 2, ipoint) = 0.d0
           shell_vgl(ishell, 3, ipoint) = 0.d0
           shell_vgl(ishell, 4, ipoint) = 0.d0
           shell_vgl(ishell, 5, ipoint) = 0.d0

           iprim_start = shell_prim_index(ishell) + 1
           iprim_end   = shell_prim_index(ishell) + shell_prim_num(ishell)

           do iprim = iprim_start, iprim_end

              ar2 = expo(iprim)*r2
              if (ar2 > cutoff) then
                 cycle
              end if

              v = coef_normalized(iprim) * dexp(-ar2)
              two_a = -2.d0 * expo(iprim) * v

              shell_vgl(ishell, 1, ipoint) = &
                   shell_vgl(ishell, 1, ipoint) + v

              shell_vgl(ishell, 2, ipoint) = &
                   shell_vgl(ishell, 2, ipoint) + two_a * x

              shell_vgl(ishell, 3, ipoint) = &
                   shell_vgl(ishell, 3, ipoint) + two_a * y

              shell_vgl(ishell, 4, ipoint) = &
                   shell_vgl(ishell, 4, ipoint) + two_a * z

              shell_vgl(ishell, 5, ipoint) = &
                   shell_vgl(ishell, 5, ipoint) + two_a * (3.d0 - 2.d0*ar2)

           end do

        end do
     end do

  end do

end function qmckl_compute_ao_basis_shell_gaussian_vgl_f
   #+end_src

   #+CALL: generate_c_interface(table=qmckl_ao_basis_shell_gaussian_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_ao_basis_shell_gaussian_vgl"))

   #+RESULTS:
   #+begin_src f90 :tangle (eval f) :comments org :exports none
   integer(c_int32_t) function qmckl_compute_ao_basis_shell_gaussian_vgl &
       (context, &
        prim_num, &
        shell_num, &
        point_num, &
        nucl_num, &
        nucleus_shell_num, &
        nucleus_index, &
        nucleus_range, &
        shell_prim_index, &
        shell_prim_num, &
        coord, &
        nucl_coord, &
        expo, &
        coef_normalized, &
        shell_vgl) &
       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 :: prim_num
     integer (c_int64_t) , intent(in)  , value :: shell_num
     integer (c_int64_t) , intent(in)  , value :: point_num
     integer (c_int64_t) , intent(in)  , value :: nucl_num
     integer (c_int64_t) , intent(in)          :: nucleus_shell_num(nucl_num)
     integer (c_int64_t) , intent(in)          :: nucleus_index(nucl_num)
     real    (c_double ) , intent(in)          :: nucleus_range(nucl_num)
     integer (c_int64_t) , intent(in)          :: shell_prim_index(shell_num)
     integer (c_int64_t) , intent(in)          :: shell_prim_num(shell_num)
     real    (c_double ) , intent(in)          :: coord(point_num,3)
     real    (c_double ) , intent(in)          :: nucl_coord(nucl_num,3)
     real    (c_double ) , intent(in)          :: expo(prim_num)
     real    (c_double ) , intent(in)          :: coef_normalized(prim_num)
     real    (c_double ) , intent(out)         :: shell_vgl(shell_num,5,point_num)

     integer(c_int32_t), external :: qmckl_compute_ao_basis_shell_gaussian_vgl_f
     info = qmckl_compute_ao_basis_shell_gaussian_vgl_f &
            (context, &
        prim_num, &
        shell_num, &
        point_num, &
        nucl_num, &
        nucleus_shell_num, &
        nucleus_index, &
        nucleus_range, &
        shell_prim_index, &
        shell_prim_num, &
        coord, &
        nucl_coord, &
        expo, &
        coef_normalized, &
        shell_vgl)

   end function qmckl_compute_ao_basis_shell_gaussian_vgl
   #+end_src

*** Provide                                                        :noexport:

    #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_provide_ao_basis_shell_vgl(qmckl_context context);
    #+end_src

    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
qmckl_exit_code qmckl_provide_ao_basis_shell_vgl(qmckl_context context)
{

  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_CONTEXT,
                           "qmckl_provide_ao_basis_shell_vgl",
                           NULL);
  }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

  if (!ctx->ao_basis.provided) {
    return qmckl_failwith( context,
                           QMCKL_NOT_PROVIDED,
                           "qmckl_provide_ao_basis_shell_vgl",
                           NULL);
  }

  /* Compute if necessary */
  if (ctx->point.date > ctx->ao_basis.shell_vgl_date) {

    /* Allocate array */
    if (ctx->ao_basis.shell_vgl == NULL) {

      qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
      mem_info.size = ctx->ao_basis.shell_num * 5 * ctx->point.num * sizeof(double);
      double* shell_vgl = (double*) qmckl_malloc(context, mem_info);

      if (shell_vgl == NULL) {
        return qmckl_failwith( context,
                               QMCKL_ALLOCATION_FAILED,
                               "qmckl_ao_basis_shell_vgl",
                               NULL);
      }
      ctx->ao_basis.shell_vgl = shell_vgl;
    }

    qmckl_exit_code rc;
    if (ctx->ao_basis.type == 'G') {
      rc = qmckl_compute_ao_basis_shell_gaussian_vgl(context,
                                                     ctx->ao_basis.prim_num,
                                                     ctx->ao_basis.shell_num,
                                                     ctx->point.num,
                                                     ctx->nucleus.num,
                                                     ctx->ao_basis.nucleus_shell_num,
                                                     ctx->ao_basis.nucleus_index,
                                                     ctx->ao_basis.nucleus_range,
                                                     ctx->ao_basis.shell_prim_index,
                                                     ctx->ao_basis.shell_prim_num,
                                                     ctx->point.coord.data,
                                                     ctx->nucleus.coord.data,
                                                     ctx->ao_basis.exponent,
                                                     ctx->ao_basis.coefficient_normalized,
                                                     ctx->ao_basis.shell_vgl);
    } else {
      return qmckl_failwith( context,
                             QMCKL_FAILURE,
                             "compute_ao_basis_shell_vgl",
                             "Not yet implemented");
    }
    if (rc != QMCKL_SUCCESS) {
      return rc;
    }

    ctx->ao_basis.shell_vgl_date = ctx->date;
  }

  return QMCKL_SUCCESS;
}
    #+end_src

*** Test                                                           :noexport:

    #+begin_src python :results output :exports none
import numpy as np

def f(a,x,y):
    return np.sum( [c * np.exp( -b*(np.linalg.norm(x-y))**2) for b,c in a] )

def df(a,x,y,n):
    h0 = 1.e-6
    if   n == 1: h = np.array([h0,0.,0.])
    elif n == 2: h = np.array([0.,h0,0.])
    elif n == 3: h = np.array([0.,0.,h0])
    return ( f(a,x+h,y) - f(a,x-h,y) ) / (2.*h0)

def d2f(a,x,y,n):
    h0 = 1.e-6
    if   n == 1: h = np.array([h0,0.,0.])
    elif n == 2: h = np.array([0.,h0,0.])
    elif n == 3: h = np.array([0.,0.,h0])
    return ( f(a,x+h,y) - 2.*f(a,x,y) + f(a,x-h,y) ) / h0**2

def lf(a,x,y):
    return d2f(a,x,y,1) + d2f(a,x,y,2) + d2f(a,x,y,3)

elec_26_w1 = np.array( [  1.49050402641, 2.90106987953, -1.05920815468  ] )
elec_15_w2 = np.array( [  -2.20180344582,-1.9113150239,  2.2193744778600002 ] )
nucl_1    = np.array( [ 1.096243353458458e+00, 8.907054016973815e-01, 7.777092280258892e-01 ] )
nucl_2    = np.array( [ 1.168459237342663e+00, 1.125660720053393e+00, 2.833370314829343e+00 ] )

#double prim_vgl[prim_num][5][point_num];
x = elec_26_w1 ; y = nucl_1
a = [(  8.236000E+03,  -1.130000E-04 * 6.1616545431994848e+02 ),
     (  1.235000E+03,  -8.780000E-04 * 1.4847738511079908e+02 ),
     (  2.808000E+02,  -4.540000E-03 * 4.8888635917437597e+01 ),
     (  7.927000E+01,  -1.813300E-02 * 1.8933972232608955e+01 ),
     (  2.559000E+01,  -5.576000E-02 * 8.1089160941724145e+00 ),
     (  8.997000E+00,  -1.268950E-01 * 3.7024003863155635e+00 ),
     (  3.319000E+00,  -1.703520E-01 * 1.7525302846177560e+00 ),
     (  9.059000E-01,   1.403820E-01 * 6.6179013183966806e-01 ),
     (  3.643000E-01,   5.986840E-01 * 3.3419848027174592e-01 ),
     (  1.285000E-01,   3.953890E-01 * 1.5296336817449557e-01 )]

print ( "[1][0][26]  : %25.15e"% f(a,x,y))
print ( "[1][1][26]  : %25.15e"% df(a,x,y,1))
print ( "[1][2][26]  : %25.15e"% df(a,x,y,2))
print ( "[1][3][26]  : %25.15e"% df(a,x,y,3))
print ( "[1][4][26]  : %25.15e"% lf(a,x,y))

    #+end_src

    #+RESULTS:
    : [1][0][26]  :     3.564393437193867e-02
    : [1][1][26]  :    -6.030177988891605e-03
    : [1][2][26]  :    -3.074832579871845e-02
    : [1][3][26]  :     2.809546963133958e-02
    : [1][4][26]  :     1.903338597841753e-02

    #+begin_src c :tangle (eval c_test) :exports none
{
#define walk_num 1 // chbrclf_walk_num
#define elec_num chbrclf_elec_num
#define shell_num chbrclf_shell_num

  int64_t elec_up_num   = chbrclf_elec_up_num;
  int64_t elec_dn_num   = chbrclf_elec_dn_num;
  double* elec_coord    = &(chbrclf_elec_coord[0][0][0]);

  rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num);
  assert (rc == QMCKL_SUCCESS);

  rc = qmckl_set_electron_walk_num (context, walk_num);
  assert (rc == QMCKL_SUCCESS);

  assert(qmckl_electron_provided(context));

  rc = qmckl_set_electron_coord (context, 'N', elec_coord, walk_num*elec_num*3);
  assert(rc == QMCKL_SUCCESS);


  double shell_vgl[elec_num][5][shell_num];

  rc = qmckl_get_ao_basis_shell_vgl(context, &(shell_vgl[0][0][0]),
        (int64_t) 5*elec_num*shell_num);
  assert (rc == QMCKL_SUCCESS);

  printf(" shell_vgl[26][0][1]  %25.15e\n", shell_vgl[26][0][1]);
  printf(" shell_vgl[26][1][1]  %25.15e\n", shell_vgl[26][1][1]);
  printf(" shell_vgl[26][2][1]  %25.15e\n", shell_vgl[26][2][1]);
  printf(" shell_vgl[26][3][1]  %25.15e\n", shell_vgl[26][3][1]);
  printf(" shell_vgl[26][4][1]  %25.15e\n", shell_vgl[26][4][1]);

  assert( fabs(shell_vgl[26][0][1] - ( 3.564393437193868e-02)) < 1.e-14 );
  assert( fabs(shell_vgl[26][1][1] - (-6.030177987072189e-03)) < 1.e-14 );
  assert( fabs(shell_vgl[26][2][1] - (-3.074832579537582e-02)) < 1.e-14 );
  assert( fabs(shell_vgl[26][3][1] - ( 2.809546963519935e-02)) < 1.e-14 );
  assert( fabs(shell_vgl[26][4][1] - ( 1.896046117183968e-02)) < 1.e-14 );

}

    #+end_src

* Polynomial part

  Going from the atomic basis set to AOs implies a systematic
  construction of all the angular functions of each shell.  We
  consider two cases for the angular functions: the real-valued
  spherical harmonics, and the polynomials in Cartesian coordinates.
  In the case of spherical harmonics, the AOs are ordered in
  increasing magnetic quantum number ($-l \le m \le l$), and in the
  case of polynomials we choose the canonical ordering, i.e

  \begin{eqnarray}
  p & : & p_x, p_y, p_z \nonumber \\
  d & : & d_{xx}, d_{xy}, d_{xz}, d_{yy}, d_{yz}, d_{zz} \nonumber \\
  f & : & f_{xxx}, f_{xxy}, f_{xxz}, f_{xyy}, f_{xyz}, f_{xzz}, f_{yyy}, f_{yyz}, f_{yzz}, f_{zzz} \nonumber \\
  {\rm etc.} \nonumber
  \end{eqnarray}

** General functions for Powers of $x-X_i$
   :PROPERTIES:
   :Name:     qmckl_ao_power
   :CRetType: qmckl_exit_code
   :FRetType: qmckl_exit_code
   :END:

   The ~qmckl_ao_power~ function 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_{ik} = X_i^k \]

   #+NAME: qmckl_ao_power_args
   | Variable  | Type            | In/Out | Description                                       |
   |-----------+-----------------+--------+---------------------------------------------------|
   | ~context~ | ~qmckl_context~ | in     | Global state                                      |
   | ~n~       | int64_t         | in     | Number of values                                  |
   | ~X~       | double[n]       | in     | Array containing the input values                 |
   | ~LMAX~    | int32_t[n]      | in     | Array containing the maximum power for each value |
   | ~P~       | double[n][ldp]  | out    | Array containing all the powers of ~X~            |
   | ~ldp~     | int64_t         | in     | Leading dimension of array ~P~                    |

   Requirements:

   - ~context~ is not ~QMCKL_NULL_CONTEXT~
   - ~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]~

   #+CALL: generate_c_header(table=qmckl_ao_power_args,rettyp=get_value("CRetType"),fname="qmckl_ao_power")

   #+RESULTS:
   #+begin_src c :tangle (eval h_func) :comments org
   qmckl_exit_code qmckl_ao_power (
         const qmckl_context context,
         const int64_t n,
         const double* X,
         const int32_t* LMAX,
         double* const P,
         const int64_t ldp );
   #+end_src

   #+begin_src f90 :tangle (eval f)
integer function qmckl_ao_power_f(context, n, X, LMAX, P, ldp) result(info)
  use qmckl
  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,k

  info = QMCKL_SUCCESS

  if (context == QMCKL_NULL_CONTEXT) then
     info = QMCKL_INVALID_CONTEXT
     return
  endif

  if (n <= ldp) then
     info = QMCKL_INVALID_ARG_2
     return
  endif

  k = MAXVAL(LMAX)
  if (LDP < k) then
     info = QMCKL_INVALID_ARG_6
     return
  endif

  if (k <= 0) then
     info = QMCKL_INVALID_ARG_4
     return
  endif

  do i=1,n
     P(1,i) = X(i)
     do k=2,LMAX(i)
        P(k,i) = P(k-1,i) * X(i)
     end do
  end do

end function qmckl_ao_power_f
   #+end_src

    #+CALL: generate_f_interface(table=qmckl_ao_power_args,rettyp=get_value("CRetType"),fname="qmckl_ao_power")

    #+RESULTS:
    #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
    interface
      integer(c_int32_t) function qmckl_ao_power &
          (context, n, X, LMAX, P, ldp) &
          bind(C)
        use, intrinsic :: iso_c_binding
        import
        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

      end function qmckl_ao_power
    end interface
    #+end_src

   #+CALL: generate_c_interface(table=qmckl_ao_power_args,rettyp=get_value("CRetType"),fname="qmckl_ao_power")

   #+RESULTS:
   #+begin_src f90 :tangle (eval f) :comments org :exports none
   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(c_int32_t), external :: qmckl_ao_power_f
     info = qmckl_ao_power_f &
            (context, n, X, LMAX, P, ldp)

   end function qmckl_ao_power
   #+end_src

*** Test                                                           :noexport:

    #+begin_src f90 :tangle (eval f_test)
integer(c_int32_t) function test_qmckl_ao_power(context) bind(C)
  use qmckl
  implicit none

  integer(qmckl_context), 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_get_numprec_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 /= QMCKL_SUCCESS) return

  test_qmckl_ao_power = QMCKL_FAILURE

  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 = QMCKL_SUCCESS
  deallocate(X,P,LMAX)
end function test_qmckl_ao_power
    #+end_src

    #+begin_src c :tangle (eval c_test) :exports none
int  test_qmckl_ao_power(qmckl_context context);
assert(0 == test_qmckl_ao_power(context));
    #+end_src

** General functions for Value, Gradient and Laplacian of a polynomial
   :PROPERTIES:
   :Name:     qmckl_ao_polynomial_vgl
   :CRetType: qmckl_exit_code
   :FRetType: qmckl_exit_code
   :END:

   A polynomial is centered on a nucleus $\mathbf{R}_i$

   \[
   P_l(\mathbf{r},\mathbf{R}_i)  =   (x-X_i)^a (y-Y_i)^b (z-Z_i)^c
   \]

   The gradients with respect to electron coordinates are

   \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*}

   and the Laplacian is

   \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_polynomial_vgl~ computes the values, gradients and
   Laplacians at a given point in space, of all polynomials with an
   angular momentum up to ~lmax~.

   #+NAME: qmckl_ao_polynomial_vgl_args
   | Variable  | Type              | In/Out | Description                                          |
   |-----------+-------------------+--------+------------------------------------------------------|
   | ~context~ | ~qmckl_context~   | in     | Global state                                         |
   | ~X~       | ~double[3]~       | in     | Array containing the coordinates of the points       |
   | ~R~       | ~double[3]~       | in     | Array containing the x,y,z coordinates of the center |
   | ~lmax~    | ~int32_t~         | in     | Maximum angular momentum                             |
   | ~n~       | ~int64_t~         | inout  | Number of computed polynomials                       |
   | ~L~       | ~int32_t[n][ldl]~ | out    | Contains a,b,c for all ~n~ results                   |
   | ~ldl~     | ~int64_t~         | in     | Leading dimension of ~L~                             |
   | ~VGL~     | ~double[n][ldv]~  | out    | Value, gradients and Laplacian of the polynomials    |
   | ~ldv~     | ~int64_t~         | in     | Leading dimension of array ~VGL~                     |

   Requirements:

    - ~context~ \ne ~QMCKL_NULL_CONTEXT~
    - ~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):
      - Increasing 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"

    #+CALL: generate_c_header(table=qmckl_ao_polynomial_vgl_args,rettyp=get_value("CRetType"),fname=get_value("Name"))

    #+RESULTS:
    #+begin_src c :tangle (eval h_func) :comments org
    qmckl_exit_code qmckl_ao_polynomial_vgl (
          const qmckl_context context,
          const double* X,
          const double* R,
          const int32_t lmax,
          int64_t* n,
          int32_t* const L,
          const int64_t ldl,
          double* const VGL,
          const int64_t ldv );
    #+end_src

    #+CALL: generate_c_header(table=qmckl_ao_polynomial_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_ao_polynomial_vgl_doc")

    #+RESULTS:
    #+begin_src c :tangle (eval h_func) :comments org
    qmckl_exit_code qmckl_ao_polynomial_vgl_doc (
          const qmckl_context context,
          const double* X,
          const double* R,
          const int32_t lmax,
          int64_t* n,
          int32_t* const L,
          const int64_t ldl,
          double* const VGL,
          const int64_t ldv );
    #+end_src

    #+begin_src c :tangle (eval c) :comments org
qmckl_exit_code
qmckl_ao_polynomial_vgl (const qmckl_context context,
                         const double* X,
                         const double* R,
                         const int32_t lmax,
                         int64_t* n,
                         int32_t* const L,
                         const int64_t ldl,
                         double* const VGL,
                         const int64_t ldv )
{
#ifdef HAVE_HPC
  //return qmckl_ao_polynomial_vgl_hpc (context, X, R, lmax, n, L, ldl, VGL, ldv);
  return qmckl_ao_polynomial_vgl_doc (context, X, R, lmax, n, L, ldl, VGL, ldv);
#else
  return qmckl_ao_polynomial_vgl_doc (context, X, R, lmax, n, L, ldl, VGL, ldv);
#endif
}
    #+end_src

    #+begin_src f90 :tangle (eval f)
integer function qmckl_ao_polynomial_vgl_doc_f (context,    &
     X, R, lmax, n, L, ldl, VGL, ldv) result(info)
  use qmckl
  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)
  real*8            :: pows(-2:lmax,3)
  double precision  :: xy, yz, xz
  double precision  :: da, db, dc, dd

  info = 0

  if (context == QMCKL_NULL_CONTEXT) then
     info = QMCKL_INVALID_CONTEXT
     return
  endif

  if (lmax < 0) then
     info = QMCKL_INVALID_ARG_4
     return
  endif

  if (ldl < 3) then
     info = QMCKL_INVALID_ARG_7
     return
  endif

  if (ldv < 5) then
     info = QMCKL_INVALID_ARG_9
     return
  endif


  do i=1,3
     Y(i) = X(i) - R(i)
  end do

  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 = QMCKL_SUCCESS

end function qmckl_ao_polynomial_vgl_doc_f
    #+end_src

    #+CALL: generate_c_interface(table=qmckl_ao_polynomial_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_ao_polynomial_vgl_doc" )

    #+RESULTS:
    #+begin_src f90 :tangle (eval f) :comments org :exports none
    integer(c_int32_t) function qmckl_ao_polynomial_vgl_doc &
        (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)
      real    (c_double ) , intent(in)          :: R(3)
      integer (c_int32_t) , intent(in)  , value :: lmax
      integer (c_int64_t) , intent(inout)        :: n
      integer (c_int32_t) , intent(out)         :: L(ldl,n)
      integer (c_int64_t) , intent(in)  , value :: ldl
      real    (c_double ) , intent(out)         :: VGL(ldv,n)
      integer (c_int64_t) , intent(in)  , value :: ldv

      integer(c_int32_t), external :: qmckl_ao_polynomial_vgl_doc_f
      info = qmckl_ao_polynomial_vgl_doc_f &
             (context, X, R, lmax, n, L, ldl, VGL, ldv)

    end function qmckl_ao_polynomial_vgl_doc
    #+end_src

    #+CALL: generate_f_interface(table=qmckl_ao_polynomial_vgl_args,rettyp=get_value("FRetType"),fname="qmckl_ao_polynomial_vgl_doc" )

    #+RESULTS:
    #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
    interface
      integer(c_int32_t) function qmckl_ao_polynomial_vgl_doc &
          (context, X, R, lmax, n, L, ldl, VGL, ldv) &
          bind(C)
        use, intrinsic :: iso_c_binding
        import
        implicit none

        integer (c_int64_t) , intent(in)  , value :: context
        real    (c_double ) , intent(in)          :: X(3)
        real    (c_double ) , intent(in)          :: R(3)
        integer (c_int32_t) , intent(in)  , value :: lmax
        integer (c_int64_t) , intent(inout)        :: n
        integer (c_int32_t) , intent(out)         :: L(ldl,n)
        integer (c_int64_t) , intent(in)  , value :: ldl
        real    (c_double ) , intent(out)         :: VGL(ldv,n)
        integer (c_int64_t) , intent(in)  , value :: ldv

      end function qmckl_ao_polynomial_vgl_doc
    end interface
    #+end_src

    #+CALL: generate_f_interface(table=qmckl_ao_polynomial_vgl_args,rettyp=get_value("FRetType"),fname="qmckl_ao_polynomial_vgl" )

    #+RESULTS:
    #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
    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
        import
        implicit none

        integer (c_int64_t) , intent(in)  , value :: context
        real    (c_double ) , intent(in)          :: X(3)
        real    (c_double ) , intent(in)          :: R(3)
        integer (c_int32_t) , intent(in)  , value :: lmax
        integer (c_int64_t) , intent(inout)        :: n
        integer (c_int32_t) , intent(out)         :: L(ldl,n)
        integer (c_int64_t) , intent(in)  , value :: ldl
        real    (c_double ) , intent(out)         :: VGL(ldv,n)
        integer (c_int64_t) , intent(in)  , value :: ldv

      end function qmckl_ao_polynomial_vgl
    end interface
    #+end_src

    #+CALL: generate_c_header(table=qmckl_ao_polynomial_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_ao_polynomial_transp_vgl")

    #+RESULTS:
    #+begin_src c :tangle (eval h_func) :comments org
    qmckl_exit_code qmckl_ao_polynomial_transp_vgl (
          const qmckl_context context,
          const double* X,
          const double* R,
          const int32_t lmax,
          int64_t* n,
          int32_t* const L,
          const int64_t ldl,
          double* const VGL,
          const int64_t ldv );
    #+end_src

    #+CALL: generate_c_header(table=qmckl_ao_polynomial_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_ao_polynomial_transp_vgl_doc")

    #+RESULTS:
    #+begin_src c :tangle (eval h_func) :comments org
    qmckl_exit_code qmckl_ao_polynomial_transp_vgl_doc (
          const qmckl_context context,
          const double* X,
          const double* R,
          const int32_t lmax,
          int64_t* n,
          int32_t* const L,
          const int64_t ldl,
          double* const VGL,
          const int64_t ldv );
    #+end_src

    #+CALL: generate_c_header(table=qmckl_ao_polynomial_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_ao_polynomial_transp_vgl_hpc")

    #+RESULTS:
    #+begin_src c :tangle (eval h_func) :comments org
    qmckl_exit_code qmckl_ao_polynomial_transp_vgl_hpc (
          const qmckl_context context,
          const double* X,
          const double* R,
          const int32_t lmax,
          int64_t* n,
          int32_t* const L,
          const int64_t ldl,
          double* const VGL,
          const int64_t ldv );
    #+end_src

    #+begin_src c :tangle (eval c) :comments org
qmckl_exit_code
qmckl_ao_polynomial_transp_vgl (const qmckl_context context,
                                const double* X,
                                const double* R,
                                const int32_t lmax,
                                int64_t* n,
                                int32_t* const L,
                                const int64_t ldl,
                                double* const VGL,
                                const int64_t ldv )
{
#ifdef HAVE_HPC
  return qmckl_ao_polynomial_transp_vgl_hpc (context, X, R, lmax, n, L, ldl, VGL, ldv);
#else
  return qmckl_ao_polynomial_transp_vgl_doc (context, X, R, lmax, n, L, ldl, VGL, ldv);
#endif
}
    #+end_src

    #+begin_src f90 :tangle (eval f)
integer function qmckl_ao_polynomial_transp_vgl_doc_f (context,    &
     X, R, lmax, n, L, ldl, VGL, ldv) result(info)
  use qmckl
  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,5)
  integer*8 , intent(in)  :: ldv

  integer*8         :: i,j
  integer           :: a,b,c,d
  real*8            :: Y(3)
  real*8            :: pows(-2:21,3) ! lmax < 22
  double precision  :: xy, yz, xz
  double precision  :: da, db, dc, dd

  info = 0

  if (context == QMCKL_NULL_CONTEXT) then
     info = QMCKL_INVALID_CONTEXT
     return
  endif

  if (lmax < 0) then
     info = QMCKL_INVALID_ARG_4
     return
  endif

  if (ldl < 3) then
     info = QMCKL_INVALID_ARG_7
     return
  endif

  if (ldv < (lmax+1)*(lmax+2)*(lmax+3)/6) then
     info = QMCKL_INVALID_ARG_9
     return
  endif


  if (lmax > 0) then

     do i=1,3
        Y(i) = X(i) - R(i)
     end do
     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

     l  (1:3,1:4) = 0
     VGL(1:4,1:5) = 0.d0

     VGL(1  ,1  ) = 1.d0

     l  (1,2) = 1
     VGL(2,1) = Y(1)
     VGL(2,2) = 1.d0

     l  (2,3) = 1
     VGL(3,1) = Y(2)
     VGL(3,3) = 1.d0

     l  (3,4) = 1
     VGL(4,1) = Y(3)
     VGL(4,4) = 1.d0

     n=4
  else
     VGL(1,1) = 1.d0
     VGL(1,2:5) = 0.d0
     l(1:3,1) = 0
     n=1
     return
  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

           xy = pows(a,1) * pows(b,2)
           yz = pows(b,2) * pows(c,3)
           xz = pows(a,1) * pows(c,3)

           l(1,n) = a
           l(2,n) = b
           l(3,n) = c

           VGL(n,1) = xy * pows(c,3)

           xy = dc * xy
           xz = db * xz
           yz = da * yz

           VGL(n,2) = pows(a-1,1) * yz
           VGL(n,3) = pows(b-1,2) * xz
           VGL(n,4) = pows(c-1,3) * xy

           VGL(n,5) = &
                (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 = QMCKL_SUCCESS

end function qmckl_ao_polynomial_transp_vgl_doc_f
    #+end_src

    #+begin_src c :tangle (eval c) :comments org
qmckl_exit_code
qmckl_ao_polynomial_transp_vgl_hpc (const qmckl_context context,
                                    const double* X,
                                    const double* R,
                                    const int32_t lmax,
                                    int64_t* n,
                                    int32_t* const L,
                                    const int64_t ldl,
                                    double* const VGL,
                                    const int64_t ldv )
{
  const qmckl_context_struct* ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL && X != NULL && R != NULL && n != NULL && L != NULL && VGL != NULL);
  if (lmax < 0) return QMCKL_INVALID_ARG_4;
  if (ldl < 3) return QMCKL_INVALID_ARG_7;

  int32_t m;

  switch (lmax) {
  case 0:
    {
      if (ldv < 1) return QMCKL_INVALID_ARG_9;
      L[0] = 0; L[1] = 0; L[2] = 0;

      VGL[0         ] = 1.0;
      VGL[ldv       ] = 0.0;
      VGL[ldv<<1    ] = 0.0;
      VGL[(ldv<<1)+ldv] = 0.0;
      VGL[ldv<<2    ] = 0.0;
      m=1;
      break;
    }
  case 1:
    {
      if (ldv < 4) return QMCKL_INVALID_ARG_9;

      if (ldl == 3) {

        const int32_t ltmp[12] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
        for (int i=0 ; i<12 ; ++i)
          L[i] = ltmp[i];

      } else {

        int32_t* const l[4] = {L, L+ldl, L+(ldl<<1), L+ldl+(ldl<<1)};
        l[0][0] = 0; l[0][1] = 0; l[0][2] = 0;
        l[1][0] = 1; l[1][1] = 0; l[1][2] = 0;
        l[2][0] = 0; l[2][1] = 1; l[2][2] = 0;
        l[3][0] = 0; l[3][1] = 0; l[3][2] = 1;
      }


      if (ldv == 4) {

        const double tmp[20] = {1.0, X[0]-R[0], X[1]-R[1], X[2]-R[2],
                                0.0, 1.0, 0.0, 0.0,
                                0.0, 0.0, 1.0, 0.0,
                                0.0, 0.0, 0.0, 1.0,
                                0.0, 0.0, 0.0, 0.0};

        for (int i=0 ; i<20 ; ++i)
          VGL[i] = tmp[i];

      } else {

        double* const vgl1 = VGL;
        double* const vgl2 = VGL + ldv;
        double* const vgl3 = VGL + (ldv << 1);
        double* const vgl4 = VGL + ldv + (ldv << 1);
        double* const vgl5 = VGL + (ldv << 2);

        for (int32_t k=0 ; k<4 ; ++k) {
          vgl2[k] = 0.0;
          vgl3[k] = 0.0;
          vgl4[k] = 0.0;
          vgl5[k] = 0.0;
        }
        vgl1[0] = 1.0;
        vgl1[1] = X[0]-R[0];
        vgl1[2] = X[1]-R[1];
        vgl1[3] = X[2]-R[2];
        vgl2[1] = 1.0;
        vgl3[2] = 1.0;
        vgl4[3] = 1.0;
      }
      m=4;
      break;
    }
  case 2:
    {
      if (ldv < 10) return QMCKL_INVALID_ARG_9;
      if (ldl == 3) {
        const int32_t ltmp[30] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1,
                                  2, 0, 0, 1, 1, 0, 1, 0, 1, 0, 2, 0,
                                  0, 1, 1, 0, 0, 2};
        for (int i=0 ; i<30 ; ++i)
          L[i] = ltmp[i];

      } else {
        int32_t* l[10];
        for (int32_t i=0 ; i<10 ; ++i) {
          l[i] = L + i*ldl;
        }

        l[0][0] = 0; l[0][1] = 0; l[0][2] = 0;
        l[1][0] = 1; l[1][1] = 0; l[1][2] = 0;
        l[2][0] = 0; l[2][1] = 1; l[2][2] = 0;
        l[3][0] = 0; l[3][1] = 0; l[3][2] = 1;
        l[4][0] = 2; l[4][1] = 0; l[4][2] = 0;
        l[5][0] = 1; l[5][1] = 1; l[5][2] = 0;
        l[6][0] = 1; l[6][1] = 0; l[6][2] = 1;
        l[7][0] = 0; l[7][1] = 2; l[7][2] = 0;
        l[8][0] = 0; l[8][1] = 1; l[8][2] = 1;
        l[9][0] = 0; l[9][1] = 0; l[9][2] = 2;
      }

      const double Y[3] = { X[0]-R[0],
                            X[1]-R[1],
                            X[2]-R[2] };

      if (ldv == 50) {
         const double tmp[50] = {
            1.0, Y[0], Y[1], Y[2], Y[0] * Y[0],
            Y[0] * Y[1], Y[0] * Y[2], Y[1] * Y[1], Y[1] * Y[2], Y[2] * Y[2],
            0.0, 1.0, 0.0, 0.0, Y[0] + Y[0],
            Y[1], Y[2], 0.0, 0.0, 0.0,
            0.0, 0.0, 1.0, 0.0, 0.0,
            Y[0], 0.0, Y[1] + Y[1], Y[2], 0.0,
            0.0, 0.0, 0.0, 1.0, 0.0,
            0.0, Y[0], 0.0, Y[1], Y[2] + Y[2],
            0.0, 0.0, 0.0, 0.0, 2.0,
            0.0, 0.0, 2.0, 0., 2.0 };
         for (int i=0 ; i<50 ; ++i)
           VGL[i] = tmp[i];
      } else {
        double* const vgl1 = VGL;
        double* const vgl2 = VGL + ldv;
        double* const vgl3 = VGL + (ldv << 1);
        double* const vgl4 = VGL + 3*ldv;
        double* const vgl5 = VGL + (ldv << 2);

        vgl1[0] = 1.0      ; vgl1[1] = Y[0]     ; vgl1[2] = Y[1];
        vgl1[3] = Y[2]     ; vgl1[4] = Y[0]*Y[0]; vgl1[5] = Y[0]*Y[1];
        vgl1[6] = Y[0]*Y[2]; vgl1[7] = Y[1]*Y[1]; vgl1[8] = Y[1]*Y[2];
        vgl1[9] = Y[2]*Y[2];

        vgl2[0] = 0.0 ; vgl2[1] = 1.0      ; vgl2[2] = 0.0 ;
        vgl2[3] = 0.0 ; vgl2[4] = Y[0]+Y[0]; vgl2[5] = Y[1];
        vgl2[6] = Y[2]; vgl2[7] = 0.0      ; vgl2[8] = 0.0 ;
        vgl2[9] = 0.0 ;

        vgl3[0] = 0.0; vgl3[1] = 0.0      ; vgl3[2] = 1.0 ;
        vgl3[3] = 0.0; vgl3[4] = 0.0      ; vgl3[5] = Y[0];
        vgl3[6] = 0.0; vgl3[7] = Y[1]+Y[1]; vgl3[8] = Y[2];
        vgl3[9] = 0.0;

        vgl4[0] = 0.0      ; vgl4[1] = 0.0; vgl4[2] = 0.0 ;
        vgl4[3] = 1.0      ; vgl4[4] = 0.0; vgl4[5] = 0.0 ;
        vgl4[6] = Y[0]     ; vgl4[7] = 0.0; vgl4[8] = Y[1];
        vgl4[9] = Y[2]+Y[2];

        vgl5[0] = 0.0; vgl5[1] = 0.0; vgl5[2] = 0.0;
        vgl5[3] = 0.0; vgl5[4] = 2.0; vgl5[5] = 0.0;
        vgl5[6] = 0.0; vgl5[7] = 2.0; vgl5[8] = 0.0;
        vgl5[9] = 2.0;
      }
      m=10;
      break;
    }
  default:
    {
      const int32_t size_max = (lmax+1)*(lmax+2)*(lmax+3)/6;
      if (ldv < size_max) return QMCKL_INVALID_ARG_9;

      double* const vgl1 = VGL;
      double* const vgl2 = VGL + ldv;
      double* const vgl3 = VGL + (ldv<<1);
      double* const vgl4 = VGL + ldv + (ldv<<1);
      double* const vgl5 = VGL + (ldv<<2);

      const double Y[3] = { X[0]-R[0],
                            X[1]-R[1],
                            X[2]-R[2] };

      assert(size_max > lmax+3);
      double pows[3][size_max];

      for (int32_t i=0 ; i<3 ; ++i) {
        pows[0][i] = 1.0;
        pows[1][i] = 1.0;
        pows[2][i] = 1.0;
      }

      for (int32_t i=3 ; i<=lmax+2 ; ++i) {
        pows[0][i] = pows[0][i-1] * Y[0];
        pows[1][i] = pows[1][i-1] * Y[1];
        pows[2][i] = pows[2][i-1] * Y[2];
      }

      int32_t* l[size_max];
      for (int32_t i=0 ; i<size_max ; ++i) {
        l[i] = &(L[i*ldl]);
      }

      for (int32_t i=0 ; i<4 ; ++i) {
        l[i][0] = 0;
        l[i][1] = 0;
        l[i][2] = 0;
      }
      l[1][0] = 1;
      l[2][1] = 1;
      l[3][2] = 1;

      for (int32_t k=0 ; k<4 ; ++k) {
        vgl2[k] = 0.0;
        vgl3[k] = 0.0;
        vgl4[k] = 0.0;
        vgl5[k] = 0.0;
      }
      vgl1[0] = 1.0;
      vgl1[1] = Y[0];
      vgl1[2] = Y[1];
      vgl1[3] = Y[2];

      vgl2[1] = 1.0;
      vgl3[2] = 1.0;
      vgl4[3] = 1.0;
      m=4;

      double dd = 2.0;

      for (int32_t d=2 ; d<= lmax ; ++d) {
        double da = dd;

        for (int32_t a=d ; a>=0 ; --a) {
          double db = dd-da;

          for (int32_t b=d-a ; b>=0 ; --b) {
            const int32_t c  = d - a - b;
            const double dc = dd - da - db;

            double xy = pows[0][a+2] * pows[1][b+2];
            double yz = pows[1][b+2] * pows[2][c+2];
            double xz = pows[0][a+2] * pows[2][c+2];

            l[m][0] = a;
            l[m][1] = b;
            l[m][2] = c;

            vgl1[m] = xy * pows[2][c+2];

            xy *= dc;
            xz *= db;
            yz *= da;

            vgl2[m] = pows[0][a+1] * yz;
            vgl3[m] = pows[1][b+1] * xz;
            vgl4[m] = pows[2][c+1] * xy;

            vgl5[m] = (da-1.) * pows[0][a] * yz +
              (db-1.) * pows[1][b] * xz +
              (dc-1.) * pows[2][c] * xy;

            db -= 1.0;
            ++m;
          }
          da -= 1.0;
        }
        dd += 1.0;
      }
    }
  }
  ,*n = m;
  return QMCKL_SUCCESS;
}
    #+end_src

    #+CALL: generate_c_interface(table=qmckl_ao_polynomial_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_ao_polynomial_transp_vgl_doc")

    #+RESULTS:
    #+begin_src f90 :tangle (eval f) :comments org :exports none
    integer(c_int32_t) function qmckl_ao_polynomial_transp_vgl_doc &
        (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)
      real    (c_double ) , intent(in)          :: R(3)
      integer (c_int32_t) , intent(in)  , value :: lmax
      integer (c_int64_t) , intent(inout)        :: n
      integer (c_int32_t) , intent(out)         :: L(ldl,n)
      integer (c_int64_t) , intent(in)  , value :: ldl
      real    (c_double ) , intent(out)         :: VGL(ldv,n)
      integer (c_int64_t) , intent(in)  , value :: ldv

      integer(c_int32_t), external :: qmckl_ao_polynomial_transp_vgl_doc_f
      info = qmckl_ao_polynomial_transp_vgl_doc_f &
             (context, X, R, lmax, n, L, ldl, VGL, ldv)

    end function qmckl_ao_polynomial_transp_vgl_doc
    #+end_src

    #+CALL: generate_f_interface(table=qmckl_ao_polynomial_vgl_args,rettyp=get_value("FRetType"),fname="qmckl_ao_polynomial_transp_vgl_doc")

    #+RESULTS:
    #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
    interface
      integer(c_int32_t) function qmckl_ao_polynomial_transp_vgl_doc &
          (context, X, R, lmax, n, L, ldl, VGL, ldv) &
          bind(C)
        use, intrinsic :: iso_c_binding
        import
        implicit none

        integer (c_int64_t) , intent(in)  , value :: context
        real    (c_double ) , intent(in)          :: X(3)
        real    (c_double ) , intent(in)          :: R(3)
        integer (c_int32_t) , intent(in)  , value :: lmax
        integer (c_int64_t) , intent(inout)        :: n
        integer (c_int32_t) , intent(out)         :: L(ldl,n)
        integer (c_int64_t) , intent(in)  , value :: ldl
        real    (c_double ) , intent(out)         :: VGL(ldv,n)
        integer (c_int64_t) , intent(in)  , value :: ldv

      end function qmckl_ao_polynomial_transp_vgl_doc
    end interface
    #+end_src

    #+CALL: generate_f_interface(table=qmckl_ao_polynomial_vgl_args,rettyp=get_value("FRetType"),fname="qmckl_ao_polynomial_transp_vgl")

    #+RESULTS:
    #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
    interface
      integer(c_int32_t) function qmckl_ao_polynomial_transp_vgl &
          (context, X, R, lmax, n, L, ldl, VGL, ldv) &
          bind(C)
        use, intrinsic :: iso_c_binding
        import
        implicit none

        integer (c_int64_t) , intent(in)  , value :: context
        real    (c_double ) , intent(in)          :: X(3)
        real    (c_double ) , intent(in)          :: R(3)
        integer (c_int32_t) , intent(in)  , value :: lmax
        integer (c_int64_t) , intent(inout)        :: n
        integer (c_int32_t) , intent(out)         :: L(ldl,n)
        integer (c_int64_t) , intent(in)  , value :: ldl
        real    (c_double ) , intent(out)         :: VGL(ldv,n)
        integer (c_int64_t) , intent(in)  , value :: ldv

      end function qmckl_ao_polynomial_transp_vgl
    end interface
    #+end_src

*** Test                                                           :noexport:

    #+begin_src f90 :tangle (eval f_test)
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_get_numprec_epsilon(context)

  X = (/ 1.1 , 2.2 ,  3.3 /)
  R = (/ 0.1 , 1.2 , -2.3 /)
  Y(:) = X(:) - R(:)

  lmax = 4;
  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 /= QMCKL_SUCCESS) return
  if (n /= d) return

  do j=1,n
     test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE
     do i=1,3
        if (L(i,j) < 0) return
     end do
     test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE
     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 = QMCKL_FAILURE
     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 = QMCKL_FAILURE
     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 = QMCKL_FAILURE
     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 = QMCKL_FAILURE
     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 = QMCKL_SUCCESS

  deallocate(L,VGL)
end function test_qmckl_ao_polynomial_vgl
    #+end_src

    #+begin_src c :tangle (eval c_test)
int  test_qmckl_ao_polynomial_vgl(qmckl_context context);
assert(0 == test_qmckl_ao_polynomial_vgl(context));

double X[3] = { 1.1, 2.2, 3.3 };
double R[3] = { 0.2, 1.1, 3.0 };
int32_t ldv[8] = {1, 4, 10, 20, 35, 56, 84, 120};
for (int32_t ldl=3 ; ldl<=5 ; ++ldl) {
    int64_t n;
    int32_t L0[200][ldl];
    int32_t L1[200][ldl];
    printf("ldl=%d\n", ldl);
    for (int32_t lmax=0 ; lmax<=7 ; lmax++) {
      double VGL0[5][ldv[lmax]];
      double VGL1[5][ldv[lmax]];
      memset(&L0[0][0], 0, sizeof(L0));
      memset(&L1[0][0], 0, sizeof(L1));
      memset(&VGL0[0][0], 0, sizeof(VGL0));
      memset(&VGL1[0][0], 0, sizeof(VGL1));
      rc = qmckl_ao_polynomial_transp_vgl_doc (context, X, R, lmax, &n, &(L0[0][0]), ldl, &(VGL0[0][0]), ldv[lmax]);
      assert (rc == QMCKL_SUCCESS);
      rc = qmckl_ao_polynomial_transp_vgl_hpc (context, X, R, lmax, &n, &(L1[0][0]), ldl, &(VGL1[0][0]), ldv[lmax]);
      assert (rc == QMCKL_SUCCESS);
      printf("lmax=%d\n", lmax);
      for (int32_t l=0 ; l<n ; ++l) {
        for (int32_t k=0 ; k<3 ; ++k) {
          printf("L[%d][%d] = %d   %d\n", l, k, L0[l][k], L1[l][k]);
          assert( L0[l][k] == L1[l][k] );
        }
      }

      for (int32_t k=0 ; k<5 ; ++k) {
        for (int32_t l=0 ; l<n ; ++l) {
          printf("VGL[%d][%d] = %e   %e  %e\n", k, l, VGL0[k][l], VGL1[k][l], VGL0[k][l]-VGL1[k][l]);
          assert( fabs(1.-(VGL0[k][l]+1.e-100)/(VGL1[k][l]+1.e-100)) < 1.e-15 );
        }
      }
    }
}
    #+end_src

* Combining radial and polynomial parts
   :PROPERTIES:
   :Name:     qmckl_compute_ao_vgl
   :CRetType: qmckl_exit_code
   :FRetType: qmckl_exit_code
   :END:

** Unoptimized version
    #+NAME: qmckl_ao_vgl_args_doc
   | Variable              | Type                              | In/Out | Description                                  |
   |-----------------------+-----------------------------------+--------+----------------------------------------------|
   | ~context~             | ~qmckl_context~                   | in     | Global state                                 |
   | ~ao_num~              | ~int64_t~                         | in     | Number of AOs                                |
   | ~shell_num~           | ~int64_t~                         | in     | Number of shells                             |
   | ~point_num~           | ~int64_t~                         | in     | Number of points                             |
   | ~nucl_num~            | ~int64_t~                         | in     | Number of nuclei                             |
   | ~coord~               | ~double[3][point_num]~            | in     | Coordinates                                  |
   | ~nucl_coord~          | ~double[3][nucl_num]~             | in     | Nuclear  coordinates                         |
   | ~nucleus_index~       | ~int64_t[nucl_num]~               | in     | Index of the 1st shell of each nucleus       |
   | ~nucleus_shell_num~   | ~int64_t[nucl_num]~               | in     | Number of shells per nucleus                 |
   | ~nucleus_range~       | ~double[nucl_num]~                | in     | Range beyond which all is zero               |
   | ~nucleus_max_ang_mom~ | ~int32_t[nucl_num]~               | in     | Maximum angular momentum per nucleus         |
   | ~shell_ang_mom~       | ~int32_t[shell_num]~              | in     | Angular momentum of each shell               |
   | ~ao_factor~           | ~double[ao_num]~                  | in     | Normalization factor of the AOs              |
   | ~shell_vgl~           | ~double[point_num][5][shell_num]~ | in     | Value, gradients and Laplacian of the shells |
   | ~ao_vgl~              | ~double[point_num][5][ao_num]~    | out    | Value, gradients and Laplacian of the AOs    |

    #+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_ao_vgl_doc_f(context, &
     ao_num, shell_num, point_num, nucl_num, &
     coord, nucl_coord, nucleus_index, nucleus_shell_num, &
     nucleus_range, nucleus_max_ang_mom, shell_ang_mom, &
     ao_factor, shell_vgl, ao_vgl) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer*8             , intent(in)  :: ao_num
  integer*8             , intent(in)  :: shell_num
  integer*8             , intent(in)  :: point_num
  integer*8             , intent(in)  :: nucl_num
  double precision      , intent(in)  :: coord(point_num,3)
  double precision      , intent(in)  :: nucl_coord(nucl_num,3)
  integer*8             , intent(in)  :: nucleus_index(nucl_num)
  integer*8             , intent(in)  :: nucleus_shell_num(nucl_num)
  double precision      , intent(in)  :: nucleus_range(nucl_num)
  integer               , intent(in)  :: nucleus_max_ang_mom(nucl_num)
  integer               , intent(in)  :: shell_ang_mom(shell_num)
  double precision      , intent(in)  :: ao_factor(ao_num)
  double precision      , intent(in)  :: shell_vgl(shell_num,5,point_num)
  double precision      , intent(out) :: ao_vgl(ao_num,5,point_num)

  double precision  :: e_coord(3), n_coord(3)
  integer*8         :: n_poly
  integer           :: l, il, k
  integer*8         :: ipoint, inucl, ishell
  integer*8         :: ishell_start, ishell_end
  integer           :: lstart(0:20)
  double precision  :: x, y, z, r2
  double precision  :: cutoff
  integer, external :: qmckl_ao_polynomial_vgl_doc_f

  double precision, allocatable  :: poly_vgl(:,:)
  integer         , allocatable  :: powers(:,:), ao_index(:)

  allocate(poly_vgl(5,ao_num), powers(3,ao_num), ao_index(ao_num))

  ! Pre-computed data
  do l=0,20
     lstart(l) = l*(l+1)*(l+2)/6 +1
  end do

  k=1
  do inucl=1,nucl_num
     ishell_start = nucleus_index(inucl) + 1
     ishell_end   = nucleus_index(inucl) + nucleus_shell_num(inucl)
     do ishell = ishell_start, ishell_end
        l = shell_ang_mom(ishell)
        ao_index(ishell) = k
        k = k + lstart(l+1) - lstart(l)
     end do
  end do
  info = QMCKL_SUCCESS

  ! Don't compute polynomials when the radial part is zero.
  cutoff = -dlog(1.d-12)

  do ipoint = 1, point_num
     e_coord(1) = coord(ipoint,1)
     e_coord(2) = coord(ipoint,2)
     e_coord(3) = coord(ipoint,3)
     do inucl=1,nucl_num
        n_coord(1) = nucl_coord(inucl,1)
        n_coord(2) = nucl_coord(inucl,2)
        n_coord(3) = nucl_coord(inucl,3)

        ! Test if the point is in the range of the nucleus
        x = e_coord(1) - n_coord(1)
        y = e_coord(2) - n_coord(2)
        z = e_coord(3) - n_coord(3)

        r2 = x*x + y*y + z*z

        if (r2 > cutoff*nucleus_range(inucl)) then
           cycle
        end if

        ! Compute polynomials
        info = qmckl_ao_polynomial_vgl_doc_f(context, e_coord, n_coord, &
             nucleus_max_ang_mom(inucl), n_poly, powers, 3_8, &
             poly_vgl, 5_8)

        ! Loop over shells
        ishell_start = nucleus_index(inucl) + 1
        ishell_end   = nucleus_index(inucl) + nucleus_shell_num(inucl)
        do ishell = ishell_start, ishell_end
           k = ao_index(ishell)
           l = shell_ang_mom(ishell)
           do il = lstart(l), lstart(l+1)-1
              ! Value
              ao_vgl(k,1,ipoint) = &
                   poly_vgl(1,il) * shell_vgl(ishell,1,ipoint) * ao_factor(k)

              ! Grad_x
              ao_vgl(k,2,ipoint) = ( &
                   poly_vgl(2,il) * shell_vgl(ishell,1,ipoint) + &
                   poly_vgl(1,il) * shell_vgl(ishell,2,ipoint) &
                   ) * ao_factor(k)

              ! Grad_y
              ao_vgl(k,3,ipoint) = ( &
                   poly_vgl(3,il) * shell_vgl(ishell,1,ipoint) + &
                   poly_vgl(1,il) * shell_vgl(ishell,3,ipoint) &
                   ) * ao_factor(k)

              ! Grad_z
              ao_vgl(k,4,ipoint) = ( &
                   poly_vgl(4,il) * shell_vgl(ishell,1,ipoint) + &
                   poly_vgl(1,il) * shell_vgl(ishell,4,ipoint) &
                   ) * ao_factor(k)

              ! Lapl_z
              ao_vgl(k,5,ipoint) = ( &
                   poly_vgl(5,il) * shell_vgl(ishell,1,ipoint) + &
                   poly_vgl(1,il) * shell_vgl(ishell,5,ipoint) + &
                   2.d0 * ( &
                   poly_vgl(2,il) * shell_vgl(ishell,2,ipoint) + &
                   poly_vgl(3,il) * shell_vgl(ishell,3,ipoint) + &
                   poly_vgl(4,il) * shell_vgl(ishell,4,ipoint) ) &
                   ) * ao_factor(k)

              k = k+1
           end do
        end do
     end do
  end do

  deallocate(poly_vgl, powers)
end function qmckl_compute_ao_vgl_doc_f
    #+end_src

** HPC version
    #+NAME: qmckl_ao_vgl_args_hpc_gaussian
   | Variable              | Type                           | In/Out | Description                                  |
   |-----------------------+--------------------------------+--------+----------------------------------------------|
   | ~context~             | ~qmckl_context~                | in     | Global state                                 |
   | ~ao_num~              | ~int64_t~                      | in     | Number of AOs                                |
   | ~shell_num~           | ~int64_t~                      | in     | Number of shells                             |
   | ~prim_num~            | ~int64_t~                      | in     | Number of primitives                         |
   | ~point_num~           | ~int64_t~                      | in     | Number of points                             |
   | ~nucl_num~            | ~int64_t~                      | in     | Number of nuclei                             |
   | ~coord~               | ~double[3][point_num]~         | in     | Coordinates                                  |
   | ~nucl_coord~          | ~double[3][nucl_num]~          | in     | Nuclear  coordinates                         |
   | ~nucleus_index~       | ~int64_t[nucl_num]~            | in     | Index of the 1st shell of each nucleus       |
   | ~nucleus_shell_num~   | ~int64_t[nucl_num]~            | in     | Number of shells per nucleus                 |
   | ~nucleus_range~       | ~double[nucl_num]~             | in     | Range beyond which all is zero               |
   | ~nucleus_max_ang_mom~ | ~int32_t[nucl_num]~            | in     | Maximum angular momentum per nucleus         |
   | ~shell_ang_mom~       | ~int32_t[shell_num]~           | in     | Angular momentum of each shell               |
   | ~shell_prim_index~    | ~int64_t[shell_num]~           | in     | Index of the 1st primitive of each shell     |
   | ~shell_prim_num~      | ~int64_t[shell_num]~           | in     | Number of primitives per shell               |
   | ~ao_factor~           | ~double[ao_num]~               | in     | Normalization factor of the AOs              |
   | ~ao_expo~             | ~double[prim_num]~             | in     | Value, gradients and Laplacian of the shells |
   | ~coef_normalized~     | ~double[prim_num]~             | in     | Value, gradients and Laplacian of the shells |
   | ~ao_vgl~              | ~double[point_num][5][ao_num]~ | out    | Value, gradients and Laplacian of the AOs    |

    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
qmckl_exit_code
qmckl_compute_ao_vgl_hpc_gaussian (
          const qmckl_context context,
          const int64_t ao_num,
          const int64_t shell_num,
          const int64_t prim_num,
          const int64_t point_num,
          const int64_t nucl_num,
          const double* coord,
          const double* nucl_coord,
          const int64_t* nucleus_index,
          const int64_t* nucleus_shell_num,
          const double* nucleus_range,
          const int32_t* nucleus_max_ang_mom,
          const int32_t* shell_ang_mom,
          const int64_t* shell_prim_index,
          const int64_t* shell_prim_num,
          const double* ao_factor,
          const double* expo,
          const double* coef_normalized,
          double* const ao_vgl )
{
  int32_t lstart[32];
  for (int32_t l=0 ; l<32 ; ++l) {
    lstart[l] = l*(l+1)*(l+2)/6;
  }

  int64_t ao_index[shell_num+1];
  int64_t size_max = 0;
  {
     int64_t k=0;
     for (int inucl=0 ; inucl < nucl_num ; ++inucl) {
        const int64_t ishell_start = nucleus_index[inucl];
        const int64_t ishell_end   = nucleus_index[inucl] + nucleus_shell_num[inucl];
        for (int64_t ishell = ishell_start ; ishell < ishell_end ; ++ishell) {
           const int l = shell_ang_mom[ishell];
           ao_index[ishell] = k;
           k += lstart[l+1] - lstart[l];
           size_max = size_max < lstart[l+1] ? lstart[l+1] : size_max;
        }
     }
     ao_index[shell_num] = ao_num+1;
  }

  /* Don't compute polynomials when the radial part is zero. */
  double cutoff = -log(1.e-12);

#ifdef HAVE_OPENMP
  #pragma omp parallel
#endif
  {
    qmckl_exit_code rc;
    double c_[prim_num];
    double expo_[prim_num];
    double ar2[prim_num];
    int32_t powers[prim_num];
    double  poly_vgl_l1[4][4] = {{1.0, 0.0, 0.0, 0.0},
                                {0.0, 1.0, 0.0, 0.0},
                                {0.0, 0.0, 1.0, 0.0},
                                {0.0, 0.0, 0.0, 1.0}};
    double  poly_vgl_l2[5][10] = {{1., 0., 0., 0., 0., 0., 0., 0., 0., 0.},
                                  {0., 1., 0., 0., 0., 0., 0., 0., 0., 0.},
                                  {0., 0., 1., 0., 0., 0., 0., 0., 0., 0.},
                                  {0., 0., 0., 1., 0., 0., 0., 0., 0., 0.},
                                  {0., 0., 0., 0., 2., 0., 0., 2., 0., 2.}};
    double  poly_vgl[5][size_max];

#ifdef HAVE_OPENMP
    #pragma omp for
#endif
    for (int64_t ipoint=0 ; ipoint < point_num ; ++ipoint) {
       const double e_coord[3] = { coord[ipoint],
                                   coord[ipoint + point_num],
                                   coord[ipoint + 2*point_num] };

       for (int64_t inucl=0 ; inucl < nucl_num ; ++inucl) {
          const double n_coord[3] = { nucl_coord[inucl],
                                      nucl_coord[inucl + nucl_num],
                                      nucl_coord[inucl + 2*nucl_num] };

          /* Test if the point is in the range of the nucleus */
          const double x = e_coord[0] - n_coord[0];
          const double y = e_coord[1] - n_coord[1];
          const double z = e_coord[2] - n_coord[2];

          const double r2 = x*x + y*y + z*z;

          if (r2 > cutoff * nucleus_range[inucl]) {
             continue;
          }

          int64_t n_poly;
          switch (nucleus_max_ang_mom[inucl]) {
          case 0:
            break;

          case 1:
            poly_vgl_l1[0][1] = x;
            poly_vgl_l1[0][2] = y;
            poly_vgl_l1[0][3] = z;
            break;

          case 2:
            poly_vgl_l2[0][1] = x;
            poly_vgl_l2[0][2] = y;
            poly_vgl_l2[0][3] = z;
            poly_vgl_l2[0][4] = x*x;
            poly_vgl_l2[0][5] = x*y;
            poly_vgl_l2[0][6] = x*z;
            poly_vgl_l2[0][7] = y*y;
            poly_vgl_l2[0][8] = y*z;
            poly_vgl_l2[0][9] = z*z;
            poly_vgl_l2[1][4] = x+x;
            poly_vgl_l2[1][5] = y;
            poly_vgl_l2[1][6] = z;
            poly_vgl_l2[2][5] = x;
            poly_vgl_l2[2][7] = y+y;
            poly_vgl_l2[2][8] = z;
            poly_vgl_l2[3][6] = x;
            poly_vgl_l2[3][8] = y;
            poly_vgl_l2[3][9] = z+z;
            break;

          default:
              rc = qmckl_ao_polynomial_transp_vgl_hpc(context, e_coord, n_coord,
                                                   nucleus_max_ang_mom[inucl],
                                                   &n_poly, powers, (int64_t) 3,
                                                   &(poly_vgl[0][0]), size_max);

              assert (rc == QMCKL_SUCCESS);
              break;
          }

          const int64_t ishell_start = nucleus_index[inucl];
          const int64_t ishell_end   = nucleus_index[inucl] + nucleus_shell_num[inucl];
          for (int64_t ishell = ishell_start ; ishell < ishell_end ; ++ishell) {

             int64_t nidx = 0;
             const int64_t iprim_start = shell_prim_index[ishell];
             const int64_t iprim_end   = shell_prim_index[ishell] + shell_prim_num[ishell];
             for (int64_t iprim = iprim_start ; iprim < iprim_end ; ++iprim) {
                const double v = expo[iprim]*r2;
                if (v <= cutoff) {
                   ar2[nidx] = v;
                   c_[nidx] = coef_normalized[iprim];
                   expo_[nidx] = expo[iprim];
                   ++nidx;
                }
             }

             double s1_ = 0.;
             double s5_ = 0.;
             double s6_ = 0.;
             for (int idx=0 ; idx<nidx ; ++idx) {
                const double v = c_[idx] * exp(-ar2[idx]);
                const double t = v*expo_[idx];
                s1_ += v;
                s6_ -= t;
                s5_ += t*(4.0*ar2[idx]-6.0);
             }
             s6_ += s6_;
             const double s1 = s1_;
             const double s2 = s6_*x;
             const double s3 = s6_*y;
             const double s4 = s6_*z;
             const double s5 = s5_;

             const int64_t k = ao_index[ishell];
             double* __restrict__ const ao_vgl_1 = ao_vgl + ipoint*5*ao_num + k;

             const int32_t l = shell_ang_mom[ishell];
             const int32_t n = lstart[l+1]-lstart[l];

             double* __restrict__ const ao_vgl_2 = ao_vgl_1 + ao_num;
             double* __restrict__ const ao_vgl_3 = ao_vgl_1 + (ao_num<<1);
             double* __restrict__ const ao_vgl_4 = ao_vgl_1 + (ao_num<<1) + ao_num;
             double* __restrict__ const ao_vgl_5 = ao_vgl_1 + (ao_num<<2);

             double* __restrict__ poly_vgl_1;
             double* __restrict__ poly_vgl_2;
             double* __restrict__ poly_vgl_3;
             double* __restrict__ poly_vgl_4;
             double* __restrict__ poly_vgl_5;
             if (nidx > 0) {
                 const double* __restrict__ f = ao_factor + k;
                 const int64_t idx = lstart[l];

                 switch (nucleus_max_ang_mom[inucl]) {
                 case 0:
                   break;
                 case 1:
                   poly_vgl_1 = &(poly_vgl_l1[0][idx]);
                   poly_vgl_2 = &(poly_vgl_l1[1][idx]);
                   poly_vgl_3 = &(poly_vgl_l1[2][idx]);
                   poly_vgl_4 = &(poly_vgl_l1[3][idx]);
                   break;
                 case 2:
                   poly_vgl_1 = &(poly_vgl_l2[0][idx]);
                   poly_vgl_2 = &(poly_vgl_l2[1][idx]);
                   poly_vgl_3 = &(poly_vgl_l2[2][idx]);
                   poly_vgl_4 = &(poly_vgl_l2[3][idx]);
                   break;
                 default:
                   poly_vgl_1 = &(poly_vgl[0][idx]);
                   poly_vgl_2 = &(poly_vgl[1][idx]);
                   poly_vgl_3 = &(poly_vgl[2][idx]);
                   poly_vgl_4 = &(poly_vgl[3][idx]);
                   poly_vgl_5 = &(poly_vgl[4][idx]);
                 }
                 switch (n) {
                 case(1):
                     ao_vgl_1[0] = s1 * f[0];
                     ao_vgl_2[0] = s2 * f[0];
                     ao_vgl_3[0] = s3 * f[0];
                     ao_vgl_4[0] = s4 * f[0];
                     ao_vgl_5[0] = s5;
                     break;
                 case (3):
#ifdef HAVE_OPENMP
                     #pragma omp simd
#endif
                     for (int il=0 ; il<3 ; ++il) {
                       ao_vgl_1[il] =  poly_vgl_1[il] * s1 * f[il];
                       ao_vgl_2[il] = (poly_vgl_2[il] * s1 + poly_vgl_1[il] * s2) * f[il];
                       ao_vgl_3[il] = (poly_vgl_3[il] * s1 + poly_vgl_1[il] * s3) * f[il];
                       ao_vgl_4[il] = (poly_vgl_4[il] * s1 + poly_vgl_1[il] * s4) * f[il];
                       ao_vgl_5[il] = (poly_vgl_1[il] * s5 +
                                       2.0*(poly_vgl_2[il] * s2 +
                                            poly_vgl_3[il] * s3 +
                                            poly_vgl_4[il] * s4 )) * f[il];
                     }
                     break;
                 case(5):
#ifdef HAVE_OPENMP
                     #pragma omp simd
#endif
                     for (int il=0 ; il<5 ; ++il) {
                       ao_vgl_1[il] =  poly_vgl_1[il] * s1 * f[il];
                       ao_vgl_2[il] = (poly_vgl_2[il] * s1 + poly_vgl_1[il] * s2) * f[il];
                       ao_vgl_3[il] = (poly_vgl_3[il] * s1 + poly_vgl_1[il] * s3) * f[il];
                       ao_vgl_4[il] = (poly_vgl_4[il] * s1 + poly_vgl_1[il] * s4) * f[il];
                       ao_vgl_5[il] = (poly_vgl_1[il] * s5 +
                                       2.0*(poly_vgl_2[il] * s2 +
                                            poly_vgl_3[il] * s3 +
                                            poly_vgl_4[il] * s4 )) * f[il];
                     }
                     break;
                 default:
#ifdef HAVE_OPENMP
                     #pragma omp simd simdlen(8)
#endif
                     for (int il=0 ; il<n ; ++il) {
                       ao_vgl_1[il] =  poly_vgl_1[il] * s1 * f[il];
                       ao_vgl_2[il] = (poly_vgl_2[il] * s1 + poly_vgl_1[il] * s2) * f[il];
                       ao_vgl_3[il] = (poly_vgl_3[il] * s1 + poly_vgl_1[il] * s3) * f[il];
                       ao_vgl_4[il] = (poly_vgl_4[il] * s1 + poly_vgl_1[il] * s4) * f[il];
                       ao_vgl_5[il] = (poly_vgl_5[il] * s1 + poly_vgl_1[il] * s5 +
                                       2.0*(poly_vgl_2[il] * s2 +
                                             poly_vgl_3[il] * s3 +
                                             poly_vgl_4[il] * s4 )) * f[il];
                     }
                     break;
                 }
             } else {
                 for (int64_t il=0 ; il<n ; ++il) {
                     ao_vgl_1[il] = 0.0;
                     ao_vgl_2[il] = 0.0;
                     ao_vgl_3[il] = 0.0;
                     ao_vgl_4[il] = 0.0;
                     ao_vgl_5[il] = 0.0;
                 }
             }
          }

       }
    }
  }
  return QMCKL_SUCCESS;
}
    #+end_src

** Interfaces
#  #+CALL: generate_c_header(table=qmckl_ao_vgl_args_doc,rettyp=get_value("CRetType"),fname="qmckl_compute_ao_vgl"))
#  (Commented because the header needs to go into h_private_func)

    #+RESULTS:
    #+begin_src c :tangle (eval h_private_func) :comments org
    qmckl_exit_code qmckl_compute_ao_vgl_doc (
          const qmckl_context context,
          const int64_t ao_num,
          const int64_t shell_num,
          const int64_t point_num,
          const int64_t nucl_num,
          const double* coord,
          const double* nucl_coord,
          const int64_t* nucleus_index,
          const int64_t* nucleus_shell_num,
          const double* nucleus_range,
          const int32_t* nucleus_max_ang_mom,
          const int32_t* shell_ang_mom,
          const double* ao_factor,
          const double* shell_vgl,
          double* const ao_vgl );
    #+end_src
    #+begin_src c :tangle (eval h_private_func) :comments org
    qmckl_exit_code qmckl_compute_ao_vgl_hpc_gaussian (
          const qmckl_context context,
          const int64_t ao_num,
          const int64_t shell_num,
          const int64_t prim_num,
          const int64_t point_num,
          const int64_t nucl_num,
          const double* coord,
          const double* nucl_coord,
          const int64_t* nucleus_index,
          const int64_t* nucleus_shell_num,
          const double* nucleus_range,
          const int32_t* nucleus_max_ang_mom,
          const int32_t* shell_ang_mom,
          const int64_t* shell_prim_index,
          const int64_t* shell_prim_num,
          const double* ao_factor,
          const double* expo,
          const double* coef_normalized,
          double* const ao_vgl );
    #+end_src

    #+CALL: generate_c_interface(table=qmckl_ao_vgl_args_doc,rettyp=get_value("CRetType"),fname="qmckl_compute_ao_vgl_doc"))

    #+RESULTS:
    #+begin_src f90 :tangle (eval f) :comments org :exports none
    integer(c_int32_t) function qmckl_compute_ao_vgl_doc &
        (context, &
         ao_num, &
         shell_num, &
         point_num, &
         nucl_num, &
         coord, &
         nucl_coord, &
         nucleus_index, &
         nucleus_shell_num, &
         nucleus_range, &
         nucleus_max_ang_mom, &
         shell_ang_mom, &
         ao_factor, &
         shell_vgl, &
         ao_vgl) &
        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 :: ao_num
      integer (c_int64_t) , intent(in)  , value :: shell_num
      integer (c_int64_t) , intent(in)  , value :: point_num
      integer (c_int64_t) , intent(in)  , value :: nucl_num
      real    (c_double ) , intent(in)          :: coord(point_num,3)
      real    (c_double ) , intent(in)          :: nucl_coord(nucl_num,3)
      integer (c_int64_t) , intent(in)          :: nucleus_index(nucl_num)
      integer (c_int64_t) , intent(in)          :: nucleus_shell_num(nucl_num)
      real    (c_double ) , intent(in)          :: nucleus_range(nucl_num)
      integer (c_int32_t) , intent(in)          :: nucleus_max_ang_mom(nucl_num)
      integer (c_int32_t) , intent(in)          :: shell_ang_mom(shell_num)
      real    (c_double ) , intent(in)          :: ao_factor(ao_num)
      real    (c_double ) , intent(in)          :: shell_vgl(shell_num,5,point_num)
      real    (c_double ) , intent(out)         :: ao_vgl(ao_num,5,point_num)

      integer(c_int32_t), external :: qmckl_compute_ao_vgl_doc_f
      info = qmckl_compute_ao_vgl_doc_f &
             (context, &
         ao_num, &
         shell_num, &
         point_num, &
         nucl_num, &
         coord, &
         nucl_coord, &
         nucleus_index, &
         nucleus_shell_num, &
         nucleus_range, &
         nucleus_max_ang_mom, &
         shell_ang_mom, &
         ao_factor, &
         shell_vgl, &
         ao_vgl)

    end function qmckl_compute_ao_vgl_doc
    #+end_src


*** Provide                                                        :noexport:

    #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context);
    #+end_src

    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context)
{

  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_CONTEXT,
                           "qmckl_provide_ao_vgl",
                           NULL);
  }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

  if (!ctx->ao_basis.provided) {
    return qmckl_failwith( context,
                           QMCKL_NOT_PROVIDED,
                           "qmckl_ao_vgl",
                           NULL);
  }

  /* Compute if necessary */
  if (ctx->point.date > ctx->ao_basis.ao_vgl_date) {

    qmckl_exit_code rc;

    /* Provide required data */
#ifndef HAVE_HPC
    rc = qmckl_provide_ao_basis_shell_vgl(context);
    if (rc != QMCKL_SUCCESS) {
        return qmckl_failwith( context, rc, "qmckl_provide_ao_basis_shell_vgl", NULL);
    }
#endif

    /* Allocate array */
    if (ctx->ao_basis.ao_vgl == NULL) {

      qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
      mem_info.size = ctx->ao_basis.ao_num * 5 * ctx->point.num * sizeof(double);
      double* ao_vgl = (double*) qmckl_malloc(context, mem_info);

      if (ao_vgl == NULL) {
        return qmckl_failwith( context,
                               QMCKL_ALLOCATION_FAILED,
                               "qmckl_ao_basis_ao_vgl",
                               NULL);
      }
      ctx->ao_basis.ao_vgl = ao_vgl;
    }
#ifdef HAVE_HPC
    if (ctx->ao_basis.type == 'G') {
      rc = qmckl_compute_ao_vgl_hpc_gaussian(context,
                                  ctx->ao_basis.ao_num,
                                  ctx->ao_basis.shell_num,
                                  ctx->ao_basis.prim_num,
                                  ctx->point.num,
                                  ctx->nucleus.num,
                                  ctx->point.coord.data,
                                  ctx->nucleus.coord.data,
                                  ctx->ao_basis.nucleus_index,
                                  ctx->ao_basis.nucleus_shell_num,
                                  ctx->ao_basis.nucleus_range,
                                  ctx->ao_basis.nucleus_max_ang_mom,
                                  ctx->ao_basis.shell_ang_mom,
                                  ctx->ao_basis.shell_prim_index,
                                  ctx->ao_basis.shell_prim_num,
                                  ctx->ao_basis.ao_factor,
                                  ctx->ao_basis.exponent,
                                  ctx->ao_basis.coefficient_normalized,
                                  ctx->ao_basis.ao_vgl);
      /*
    } else if (ctx->ao_basis.type == 'S') {
      rc = qmck_compute_ao_vgl_hpc_slater(context,
                                  ctx->ao_basis.ao_num,
                                  ctx->ao_basis.shell_num,
                                  ctx->ao_basis.prim_num,
                                  ctx->point.num,
                                  ctx->nucleus.num,
                                  ctx->point.coord.data,
                                  ctx->nucleus.coord.data,
                                  ctx->ao_basis.nucleus_index,
                                  ctx->ao_basis.nucleus_shell_num,
                                  ctx->ao_basis.nucleus_range,
                                  ctx->ao_basis.nucleus_max_ang_mom,
                                  ctx->ao_basis.shell_ang_mom,
                                  ctx->ao_basis.shell_prim_index,
                                  ctx->ao_basis.shell_prim_num,
                                  ctx->ao_basis.ao_factor,
                                  ctx->ao_basis.exponent,
                                  ctx->ao_basis.coefficient_normalized,
                                  ctx->ao_basis.ao_vgl);
      ,*/
    } else {
      rc = qmckl_compute_ao_vgl_doc(context,
                                  ctx->ao_basis.ao_num,
                                  ctx->ao_basis.shell_num,
                                  ctx->point.num,
                                  ctx->nucleus.num,
                                  ctx->point.coord.data,
                                  ctx->nucleus.coord.data,
                                  ctx->ao_basis.nucleus_index,
                                  ctx->ao_basis.nucleus_shell_num,
                                  ctx->ao_basis.nucleus_range,
                                  ctx->ao_basis.nucleus_max_ang_mom,
                                  ctx->ao_basis.shell_ang_mom,
                                  ctx->ao_basis.ao_factor,
                                  ctx->ao_basis.shell_vgl,
                                  ctx->ao_basis.ao_vgl);
    }
#else
    rc = qmckl_compute_ao_vgl_doc(context,
                                  ctx->ao_basis.ao_num,
                                  ctx->ao_basis.shell_num,
                                  ctx->point.num,
                                  ctx->nucleus.num,
                                  ctx->point.coord.data,
                                  ctx->nucleus.coord.data,
                                  ctx->ao_basis.nucleus_index,
                                  ctx->ao_basis.nucleus_shell_num,
                                  ctx->ao_basis.nucleus_range,
                                  ctx->ao_basis.nucleus_max_ang_mom,
                                  ctx->ao_basis.shell_ang_mom,
                                  ctx->ao_basis.ao_factor,
                                  ctx->ao_basis.shell_vgl,
                                  ctx->ao_basis.ao_vgl);
#endif
    if (rc != QMCKL_SUCCESS) {
      return rc;
    }

    ctx->ao_basis.ao_vgl_date = ctx->date;
  }

  return QMCKL_SUCCESS;
}
    #+end_src

*** Test                                                           :noexport:

    #+begin_src python :results output :exports none
import numpy as np
from math import sqrt

def f(a,x,y):
    return np.sum( [c * np.exp( -b*(np.linalg.norm(x-y))**2) for b,c in a] )

def df(a,x,y,n):
    h0 = 1.e-6
    if   n == 1: h = np.array([h0,0.,0.])
    elif n == 2: h = np.array([0.,h0,0.])
    elif n == 3: h = np.array([0.,0.,h0])
    return ( f(a,x+h,y) - f(a,x-h,y) ) / (2.*h0)

def d2f(a,x,y,n):
    h0 = 1.e-6
    if   n == 1: h = np.array([h0,0.,0.])
    elif n == 2: h = np.array([0.,h0,0.])
    elif n == 3: h = np.array([0.,0.,h0])
    return ( f(a,x+h,y) - 2.*f(a,x,y) + f(a,x-h,y) ) / h0**2

def lf(a,x,y):
    return d2f(a,x,y,1) + d2f(a,x,y,2) + d2f(a,x,y,3)

elec_26_w1 = np.array( [  1.49050402641, 2.90106987953, -1.05920815468  ] )
elec_15_w2 = np.array( [  -2.20180344582,-1.9113150239,  2.2193744778600002 ] )
nucl_1    = np.array( [ -2.302574592081335e+00, -3.542027060505035e-01, -5.334129934317614e-02] )

#double ao_vgl[prim_num][5][elec_num];
x = elec_26_w1 ; y = nucl_1
a = [( 403.830000, 0.001473 * 5.9876577632594533e+04),
     ( 121.170000, 0.012672 * 7.2836806319891484e+03),
     (  46.345000, 0.058045 * 1.3549226646722386e+03),
     (  19.721000, 0.170510 * 3.0376315094739988e+02),
     (   8.862400, 0.318596 * 7.4924579607137730e+01),
     (   3.996200, 0.384502 * 1.8590543353806009e+01),
     (   1.763600, 0.273774 * 4.4423176930919421e+00),
     (   0.706190, 0.074397 * 8.9541051939952665e-01)]

norm = sqrt(3.)
print ( "[0][26][219] : %25.15e"%(f(a,x,y) * (x[0] - y[0])**2) )
print ( "[1][26][219] : %25.15e"%(df(a,x,y,1)* (x[0] - y[0]) * (x[1] - y[1]) + 2.*f(a,x,y) * (x[0] - y[0])) )

print ( "[0][26][220] : %25.15e"%(norm*f(a,x,y) * (x[0] - y[0]) * (x[1] - y[1]) ))
print ( "[1][26][220] : %25.15e"%(norm*df(a,x,y,1)* (x[0] - y[0]) * (x[1] - y[1]) + norm*f(a,x,y) * (x[1] - y[1])) )

print ( "[0][26][221] : %25.15e"%(norm*f(a,x,y) * (x[0] - y[0]) * (x[2] - y[2])) )
print ( "[1][26][221] : %25.15e"%(norm*df(a,x,y,1)* (x[0] - y[0]) * (x[2] - y[2]) + norm*f(a,x,y) * (x[2] - y[2])) )

print ( "[0][26][222] : %25.15e"%(f(a,x,y) * (x[1] - y[1]) * (x[1] - y[1])) )
print ( "[1][26][222] : %25.15e"%(df(a,x,y,1)* (x[1] - y[1]) * (x[1] - y[1])) )

print ( "[0][26][223] : %25.15e"%(norm*f(a,x,y) * (x[1] - y[1]) * (x[2] - y[2])) )
print ( "[1][26][223] : %25.15e"%(norm*df(a,x,y,1)* (x[1] - y[1]) * (x[2] - y[2])) )

print ( "[0][26][224] : %25.15e"%(f(a,x,y) * (x[2] - y[2]) * (x[2] - y[2])) )
print ( "[1][26][224] : %25.15e"%(df(a,x,y,1)* (x[2] - y[2]) * (x[2] - y[2])) )

    #+end_src

    #+RESULTS:
    #+begin_example
    [0][26][219] :     1.020302912653649e-08
    [1][26][219] :    -4.153046808203204e-08
    [0][26][220] :     1.516649653540510e-08
    [1][26][220] :    -7.725252615816528e-08
    [0][26][221] :    -4.686389780112468e-09
    [1][26][221] :     2.387073693851122e-08
    [0][26][222] :     7.514847283937212e-09
    [1][26][222] :    -4.025905373647693e-08
    [0][26][223] :    -4.021924592380977e-09
    [1][26][223] :     2.154652944642284e-08
    [0][26][224] :     7.175074806631758e-10
    [1][26][224] :    -3.843880138733679e-09
    #+end_example

     #+begin_src c :tangle (eval c_test) :exports none
{
#define walk_num 1 // chbrclf_walk_num
#define elec_num chbrclf_elec_num
#define shell_num chbrclf_shell_num
#define ao_num chbrclf_ao_num

int64_t elec_up_num   = chbrclf_elec_up_num;
int64_t elec_dn_num   = chbrclf_elec_dn_num;
double* elec_coord    = &(chbrclf_elec_coord[0][0][0]);

rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num);
assert (rc == QMCKL_SUCCESS);

rc = qmckl_set_electron_walk_num (context, walk_num);
assert (rc == QMCKL_SUCCESS);

assert(qmckl_electron_provided(context));

rc = qmckl_set_electron_coord (context, 'N', elec_coord, walk_num*elec_num*3);
assert(rc == QMCKL_SUCCESS);


double ao_vgl[elec_num][5][ao_num];

rc = qmckl_get_ao_basis_ao_vgl(context, &(ao_vgl[0][0][0]),
         (int64_t) 5*elec_num*ao_num);
assert (rc == QMCKL_SUCCESS);

printf("\n");
printf(" ao_vgl ao_vgl[0][26][219] %25.15e\n", ao_vgl[26][0][219]);
printf(" ao_vgl ao_vgl[1][26][219] %25.15e\n", ao_vgl[26][1][219]);
printf(" ao_vgl ao_vgl[0][26][220] %25.15e\n", ao_vgl[26][0][220]);
printf(" ao_vgl ao_vgl[1][26][220] %25.15e\n", ao_vgl[26][1][220]);
printf(" ao_vgl ao_vgl[0][26][221] %25.15e\n", ao_vgl[26][0][221]);
printf(" ao_vgl ao_vgl[1][26][221] %25.15e\n", ao_vgl[26][1][221]);
printf(" ao_vgl ao_vgl[0][26][222] %25.15e\n", ao_vgl[26][0][222]);
printf(" ao_vgl ao_vgl[1][26][222] %25.15e\n", ao_vgl[26][1][222]);
printf(" ao_vgl ao_vgl[0][26][223] %25.15e\n", ao_vgl[26][0][223]);
printf(" ao_vgl ao_vgl[1][26][223] %25.15e\n", ao_vgl[26][1][223]);
printf(" ao_vgl ao_vgl[0][26][224] %25.15e\n", ao_vgl[26][0][224]);
printf(" ao_vgl ao_vgl[1][26][224] %25.15e\n", ao_vgl[26][1][224]);
printf("\n");

assert( fabs(ao_vgl[26][0][219] - ( 1.020298798341620e-08)) < 1.e-14 );
assert( fabs(ao_vgl[26][1][219] - (-4.928035238010602e-08)) < 1.e-14 );
assert( fabs(ao_vgl[26][0][220] - ( 1.516643537739178e-08)) < 1.e-14 );
assert( fabs(ao_vgl[26][1][220] - (-7.725221462603871e-08)) < 1.e-14 );
assert( fabs(ao_vgl[26][0][221] - (-4.686370882518819e-09)) < 1.e-14 );
assert( fabs(ao_vgl[26][1][221] - ( 2.387064067626827e-08)) < 1.e-14 );
assert( fabs(ao_vgl[26][0][222] - ( 7.514816980753531e-09)) < 1.e-14 );
assert( fabs(ao_vgl[26][1][222] - (-4.025889138635182e-08)) < 1.e-14 );
assert( fabs(ao_vgl[26][0][223] - (-4.021908374204471e-09)) < 1.e-14 );
assert( fabs(ao_vgl[26][1][223] - ( 2.154644255710413e-08)) < 1.e-14 );
assert( fabs(ao_vgl[26][0][224] - ( 7.175045873560788e-10)) < 1.e-14 );
assert( fabs(ao_vgl[26][1][224] - (-3.843864637762753e-09)) < 1.e-14 );
}

     #+end_src

* End of files                                                     :noexport:

  #+begin_src c :tangle (eval h_private_type)
#endif
  #+end_src

  #+begin_src c :tangle (eval h_private_func)
#endif
  #+end_src

*** Test
  #+begin_src c :tangle (eval c_test)
    rc = qmckl_context_destroy(context);
    assert (rc == QMCKL_SUCCESS);

    return 0;
}
  #+end_src

*** Compute file names
    #+begin_src emacs-lisp
; The following is required to compute the file names

(setq pwd (file-name-directory buffer-file-name))
(setq name (file-name-nondirectory (substring buffer-file-name 0 -4)))
(setq f  (concat pwd name "_f.f90"))
(setq fh (concat pwd name "_fh.f90"))
(setq c  (concat pwd name ".c"))
(setq h  (concat name ".h"))
(setq h_private  (concat name "_private.h"))
(setq c_test  (concat pwd "test_" name ".c"))
(setq f_test  (concat pwd "test_" name "_f.f90"))

; Minted
(require 'ox-latex)
(setq org-latex-listings 'minted)
(add-to-list 'org-latex-packages-alist '("" "listings"))
(add-to-list 'org-latex-packages-alist '("" "color"))

    #+end_src

    #+RESULTS:
    |   | color    |
    |   | listings |


# -*- mode: org -*-
# vim: syntax=c


* TODO [0/1] Missing features                                      :noexport:
  - [ ] Error messages to tell what is missing when initializing