Possibility to compute only values

org/qmckl_ao.org

@@ -299,6 +299,8 @@ typedef struct qmckl_ao_basis_struct {
   uint64_t           shell_vgl_date;
   double  * restrict ao_vgl;
   uint64_t           ao_vgl_date;
+  double  * restrict ao_value;
+  uint64_t           ao_value_date;
 
   int32_t            uninitialized;
   bool               provided;
@@ -2490,8 +2492,10 @@ free(ao_factor_test);
    | ~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~             | ~double[point_num][5][ao_num]~    | Value, gradients, Laplacian of the AOs at current  positions                                 |
    | ~ao_vgl_date~        | ~uint64_t~                        | Last modification date of Value, gradients, Laplacian of the AOs at current positions        |
+   | ~ao_value~           | ~double[point_num][ao_num]~       | Values of the the AOs at current positions                                                   |
+   | ~ao_value_date~      | ~uint64_t~                        | Last modification date of the values of the AOs at current positions                         |
    |----------------------+-----------------------------------+----------------------------------------------------------------------------------------------|
 
 *** After initialization
@@ -3022,7 +3026,7 @@ qmckl_get_ao_basis_ao_vgl (qmckl_context context,
   end interface
     #+end_src
 
-    Uses the give array to compute the VGL.
+    Uses the given array to compute the VGL.
 
     #+begin_src c :comments org :tangle (eval h_func) :noweb yes
 qmckl_exit_code
@@ -3088,6 +3092,133 @@ qmckl_get_ao_basis_ao_vgl_inplace (qmckl_context context,
   end interface
     #+end_src
 
+
+    
+    #+begin_src c :comments org :tangle (eval h_func) :noweb yes
+qmckl_exit_code
+qmckl_get_ao_basis_ao_value (qmckl_context context,
+                             double* const ao_value,
+                             const int64_t size_max);
+    #+end_src
+
+    Returns the array of values 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_value (qmckl_context context,
+                             double* const ao_value,
+                             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_value",
+                           NULL);
+  }
+
+  qmckl_exit_code rc;
+
+  rc = qmckl_provide_ao_value(context);
+  if (rc != QMCKL_SUCCESS) return rc;
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
+  assert (ctx != NULL);
+
+  int64_t sze = ctx->ao_basis.ao_num * ctx->point.num;
+  if (size_max < sze) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_3,
+                           "qmckl_get_ao_basis_ao_value",
+                           "input array too small");
+  }
+  memcpy(ao_value, ctx->ao_basis.ao_value, (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_value (context, &
+          ao_value, 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_value(*)
+       integer (c_int64_t) , intent(in)  , value :: size_max
+     end function qmckl_get_ao_basis_ao_value
+  end interface
+    #+end_src
+
+    Uses the given array to compute the value.
+
+    #+begin_src c :comments org :tangle (eval h_func) :noweb yes
+qmckl_exit_code
+qmckl_get_ao_basis_ao_value_inplace (qmckl_context context,
+                                     double* const ao_value,
+                                     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_value_inplace (qmckl_context context,
+                                     double* const ao_value,
+                                     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_value",
+                           NULL);
+  }
+
+  qmckl_exit_code rc;
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
+  assert (ctx != NULL);
+
+  int64_t sze = ctx->ao_basis.ao_num * ctx->point.num;
+  if (size_max < sze) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_3,
+                           "qmckl_get_ao_basis_ao_value",
+                           "input array too small");
+  }
+
+  rc = qmckl_context_touch(context);
+  if (rc != QMCKL_SUCCESS) return rc;
+
+  double* old_array = ctx->ao_basis.ao_value;
+
+  ctx->ao_basis.ao_value = ao_value;
+
+  rc = qmckl_provide_ao_value(context);
+  if (rc != QMCKL_SUCCESS) return rc;
+
+  ctx->ao_basis.ao_value = 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_value_inplace (context, &
+          ao_value, 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_value(*)
+       integer (c_int64_t) , intent(in)  , value :: size_max
+     end function qmckl_get_ao_basis_ao_value_inplace
+  end interface
+    #+end_src
+
 * Radial part
 
 ** General functions for Gaussian basis functions
@@ -5285,13 +5416,760 @@ for (int32_t ldl=3 ; ldl<=5 ; ++ldl) {
     #+end_src
 
 * Combining radial and polynomial parts
+
+** Values only
+   :PROPERTIES:
+   :Name:     qmckl_compute_ao_value
+   :CRetType: qmckl_exit_code
+   :FRetType: qmckl_exit_code
+   :END:
+*** Unoptimized version
+     #+NAME: qmckl_ao_value_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_value~            | ~double[point_num][ao_num]~       | out    | Values of the AOs                            |
+
+     #+begin_src f90 :comments org :tangle (eval f) :noweb yes
+integer function qmckl_compute_ao_value_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_value) &
+     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_value(ao_num,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_value(k,ipoint) = &
+                   poly_vgl(1,il) * shell_vgl(ishell,1,ipoint) * ao_factor(k)
+              k = k+1
+           end do
+        end do
+     end do
+  end do
+
+  deallocate(poly_vgl, powers)
+end function qmckl_compute_ao_value_doc_f
+     #+end_src
+
+*** HPC version
+     #+NAME: qmckl_ao_value_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_value~            | ~double[point_num][ao_num]~ | out    | Values of the AOs                            |
+
+     #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
+#ifdef HAVE_HPC
+qmckl_exit_code
+qmckl_compute_ao_value_hpc_gaussian (const qmckl_context context,
+                                     const int64_t ao_num,
+                                     const int64_t shell_num,
+                                     const int32_t* restrict prim_num_per_nucleus,
+                                     const int64_t point_num,
+                                     const int64_t nucl_num,
+                                     const double* restrict coord,
+                                     const double* restrict nucl_coord,
+                                     const int64_t* restrict nucleus_index,
+                                     const int64_t* restrict nucleus_shell_num,
+                                     const double* nucleus_range,
+                                     const int32_t* restrict nucleus_max_ang_mom,
+                                     const int32_t* restrict shell_ang_mom,
+                                     const double* restrict ao_factor,
+                                     const qmckl_matrix expo_per_nucleus,
+                                     const qmckl_tensor coef_per_nucleus,
+                                     double* restrict const ao_value )
+{
+  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 prim_max = 0;
+  int64_t shell_max = 0;
+  {
+    int64_t k=0;
+    for (int inucl=0 ; inucl < nucl_num ; ++inucl) {
+      prim_max = prim_num_per_nucleus[inucl] > prim_max ?
+        prim_num_per_nucleus[inucl] : prim_max;
+      shell_max = nucleus_shell_num[inucl] > shell_max ?
+        nucleus_shell_num[inucl] : shell_max;
+      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 ar2[prim_max];
+    int32_t powers[3*size_max];
+    double  poly_vgl[5*size_max];
+
+    double exp_mat[prim_max];
+    double ce_mat[shell_max];
+
+    double coef_mat[nucl_num][shell_max][prim_max];
+    for (int i=0 ; i<nucl_num ; ++i) {
+      for (int j=0 ; j<shell_max; ++j) {
+        for (int k=0 ; k<prim_max; ++k) {
+          coef_mat[i][j][k] = qmckl_ten3(coef_per_nucleus,k, j, i);
+        }
+      }
+    }
+
+#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[0] = 0.;
+          poly_vgl[1] = x;
+          poly_vgl[2] = y;
+          poly_vgl[3] = z;
+          break;
+
+        case 2:
+          poly_vgl[0] = 0.;
+          poly_vgl[1] = x;
+          poly_vgl[2] = y;
+          poly_vgl[3] = z;
+          poly_vgl[4] = x*x;
+          poly_vgl[5] = x*y;
+          poly_vgl[6] = x*z;
+          poly_vgl[7] = y*y;
+          poly_vgl[8] = y*z;
+          poly_vgl[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, size_max);
+          assert (rc == QMCKL_SUCCESS);
+          break;
+        }
+
+ /* Compute all exponents */
+
+        int64_t nidx = 0;
+        for (int64_t iprim = 0 ; iprim < prim_num_per_nucleus[inucl] ; ++iprim) {
+          const double v = qmckl_mat(expo_per_nucleus, iprim, inucl) * r2;
+          if (v <= cutoff) {
+            ar2[iprim] = v;
+            ++nidx;
+          } else {
+            break;
+          }
+        }
+
+        for (int64_t iprim = 0 ; iprim < nidx ; ++iprim) {
+          exp_mat[iprim] = exp(-ar2[iprim]);
+        }
+
+        for (int i=0 ; i<nucleus_shell_num[inucl] ; ++i) {
+            ce_mat[i] = 0.;
+        }
+        for (int k=0 ; k<nidx; ++k) {
+          for (int i=0 ; i<nucleus_shell_num[inucl] ; ++i) {
+            ce_mat[i] += coef_mat[inucl][i][k] * exp_mat[k];
+          }
+        }
+
+        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 double s1 = ce_mat[ishell-ishell_start];
+          if (s1 == 0.0) continue;
+
+          const int64_t k = ao_index[ishell];
+          double* restrict const ao_value_1 = ao_value + ipoint*ao_num + k;
+
+          const int32_t l = shell_ang_mom[ishell];
+          const int32_t n = lstart[l+1]-lstart[l];
+
+          double* restrict poly_vgl_1 = NULL;
+          if (nidx > 0) {
+            const double* restrict f = ao_factor + k;
+            const int64_t idx = lstart[l];
+
+            poly_vgl_1 = &(poly_vgl[idx]);
+
+            switch (n) {
+            case(1):
+              ao_value_1[0] = s1 * f[0];
+              break;
+            case (3):
+#ifdef HAVE_OPENMP
+#pragma omp simd
+#endif
+              for (int il=0 ; il<3 ; ++il) {
+                ao_value_1[il] =  poly_vgl_1[il] * s1 * f[il];
+              }
+              break;
+            case(6):
+#ifdef HAVE_OPENMP
+#pragma omp simd
+#endif
+              for (int il=0 ; il<6 ; ++il) {
+                ao_value_1[il] =  poly_vgl_1[il] * s1 * f[il];
+              }
+              break;
+            default:
+#ifdef HAVE_OPENMP
+#pragma omp simd simdlen(8)
+#endif
+              for (int il=0 ; il<n ; ++il) {
+                ao_value_1[il] =  poly_vgl_1[il] * s1 * f[il];
+              }
+              break;
+            }
+          } else {
+            for (int64_t il=0 ; il<n ; ++il) {
+              ao_value_1[il] = 0.0;
+            }
+          }
+        }
+      }
+    }
+  }
+
+  return QMCKL_SUCCESS;
+}
+#endif
+     #+end_src
+
+*** Interfaces
+ #  #+CALL: generate_c_header(table=qmckl_ao_value_args_doc,rettyp=get_value("CRetType"),fname="qmckl_compute_ao_value"))
+ #  (Commented because the header needs to go into h_private_func)
+
+     #+begin_src c :tangle (eval h_private_func) :comments org
+    qmckl_exit_code qmckl_compute_ao_value_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_value );
+     #+end_src
+
+     #+begin_src c :tangle (eval h_private_func) :comments org
+#ifdef HAVE_HPC
+    qmckl_exit_code qmckl_compute_ao_value_hpc_gaussian (
+          const qmckl_context context,
+          const int64_t ao_num,
+          const int64_t shell_num,
+          const int32_t* prim_num_per_nucleus,
+          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 qmckl_matrix expo_per_nucleus,
+          const qmckl_tensor coef_per_nucleus,
+          double* const ao_value );
+#endif
+     #+end_src
+
+     #+CALL: generate_c_interface(table=qmckl_ao_value_args_doc,rettyp=get_value("CRetType"),fname="qmckl_compute_ao_value_doc"))
+
+     #+RESULTS:
+     #+begin_src f90 :tangle (eval f) :comments org :exports none
+    integer(c_int32_t) function qmckl_compute_ao_value_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_value) &
+        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_value(ao_num,point_num)
+
+      integer(c_int32_t), external :: qmckl_compute_ao_value_doc_f
+      info = qmckl_compute_ao_value_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_value)
+
+    end function qmckl_compute_ao_value_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_value(qmckl_context context);
+     #+end_src
+
+     #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
+qmckl_exit_code qmckl_provide_ao_value(qmckl_context context)
+{
+
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_CONTEXT,
+                           "qmckl_provide_ao_value",
+                           NULL);
+  }
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
+  assert (ctx != NULL);
+
+  if (!ctx->ao_basis.provided) {
+    return qmckl_failwith( context,
+                           QMCKL_NOT_PROVIDED,
+                           "qmckl_ao_value",
+                           NULL);
+  }
+
+ /* Compute if necessary */
+  if (ctx->point.date > ctx->ao_basis.ao_value_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_value == 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_value = (double*) qmckl_malloc(context, mem_info);
+
+      if (ao_value == NULL) {
+        return qmckl_failwith( context,
+                               QMCKL_ALLOCATION_FAILED,
+                               "qmckl_ao_basis_ao_value",
+                               NULL);
+      }
+      ctx->ao_basis.ao_value = ao_value;
+    }
+
+    if (ctx->ao_basis.ao_vgl_date == ctx->point.date) {
+
+      // ao_vgl has been computed at this step: Just copy the data.
+
+      double * v = &(ctx->ao_basis.ao_value[0]);
+      double * vgl = &(ctx->ao_basis.ao_vgl[0]);
+      for (int i=0 ; i<ctx->point.num ; ++i) {
+        for (int k=0 ; k<ctx->ao_basis.ao_num ; ++k) {
+          v[k] = vgl[k];
+        }
+        v   += ctx->ao_basis.ao_num;
+        vgl += ctx->ao_basis.ao_num * 5;
+      }
+
+    } else {
+
+#ifdef HAVE_HPC
+      if (ctx->ao_basis.type == 'G') {
+        rc = qmckl_compute_ao_value_hpc_gaussian(context,
+                                                 ctx->ao_basis.ao_num,
+                                                 ctx->ao_basis.shell_num,
+                                                 ctx->ao_basis.prim_num_per_nucleus,
+                                                 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.expo_per_nucleus,
+                                                 ctx->ao_basis.coef_per_nucleus,
+                                                 ctx->ao_basis.ao_value);
+        /*
+          } else if (ctx->ao_basis.type == 'S') {
+          rc = qmck_compute_ao_value_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_value);
+        ,*/
+      } else {
+        rc = qmckl_compute_ao_value_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_value);
+      }
+#else
+      rc = qmckl_compute_ao_value_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_value);
+#endif
+      if (rc != QMCKL_SUCCESS) {
+        return rc;
+      }
+
+    }
+
+    ctx->ao_basis.ao_value_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
+
+h0 = 1.e-4
+def f(a,x,y):
+    return np.sum( [c * np.exp( -b*(np.linalg.norm(x-y))**2) for b,c in a] )
+
+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_value[prim_num][5][elec_num];
+x = elec_26_w1 ; y = nucl_1
+a = [( 4.0382999999999998e+02, 1.4732000000000000e-03 * 5.9876577632594533e+04),
+     ( 1.2117000000000000e+02, 1.2672500000000000e-02 * 7.2836806319891484e+03),
+     ( 4.6344999999999999e+01, 5.8045100000000002e-02 * 1.3549226646722386e+03),
+     ( 1.9721000000000000e+01, 1.7051030000000000e-01 * 3.0376315094739988e+02),
+     ( 8.8623999999999992e+00, 3.1859579999999998e-01 * 7.4924579607137730e+01),
+     ( 3.9962000000000000e+00, 3.8450230000000002e-01 * 1.8590543353806009e+01),
+     ( 1.7636000000000001e+00, 2.7377370000000001e-01 * 4.4423176930919421e+00),
+     ( 7.0618999999999998e-01, 7.4396699999999996e-02 * 8.9541051939952665e-01)]
+
+norm = sqrt(3.)
+# x^2 * g(r)
+print ( "[26][0][219] : %25.15e"%(fx(a,x,y)) )
+print ( "[26][1][219] : %25.15e"%(df(a,x,y,1)) )
+print ( "[26][2][219] : %25.15e"%(df(a,x,y,2)) )
+print ( "[26][3][219] : %25.15e"%(df(a,x,y,3)) )
+print ( "[26][4][219] : %25.15e"%(lf(a,x,y)) )
+
+print ( "[26][0][220] : %25.15e"%(norm*f(a,x,y) * (x[0] - y[0]) * (x[1] - y[1]) ))
+print ( "[26][1][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 ( "[26][0][221] : %25.15e"%(norm*f(a,x,y) * (x[0] - y[0]) * (x[2] - y[2])) )
+print ( "[26][1][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 ( "[26][0][222] : %25.15e"%(f(a,x,y) * (x[1] - y[1]) * (x[1] - y[1])) )
+print ( "[26][1][222] : %25.15e"%(df(a,x,y,1)* (x[1] - y[1]) * (x[1] - y[1])) )
+
+print ( "[26][0][223] : %25.15e"%(norm*f(a,x,y) * (x[1] - y[1]) * (x[2] - y[2])) )
+print ( "[26][1][223] : %25.15e"%(norm*df(a,x,y,1)* (x[1] - y[1]) * (x[2] - y[2])) )
+
+print ( "[26][0][224] : %25.15e"%(f(a,x,y) * (x[2] - y[2]) * (x[2] - y[2])) )
+print ( "[26][1][224] : %25.15e"%(df(a,x,y,1)* (x[2] - y[2]) * (x[2] - y[2])) )
+
+     #+end_src
+
+     #+RESULTS:
+
+      #+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_value[elec_num][ao_num];
+
+rc = qmckl_get_ao_basis_ao_value(context, &(ao_value[0][0]),
+         (int64_t) elec_num*ao_num);
+assert (rc == QMCKL_SUCCESS);
+
+printf("\n");
+printf(" ao_value ao_value[26][219] %25.15e\n", ao_value[26][219]);
+printf(" ao_value ao_value[26][220] %25.15e\n", ao_value[26][220]);
+printf(" ao_value ao_value[26][221] %25.15e\n", ao_value[26][221]);
+printf(" ao_value ao_value[26][222] %25.15e\n", ao_value[26][222]);
+printf(" ao_value ao_value[26][223] %25.15e\n", ao_value[26][223]);
+printf(" ao_value ao_value[26][224] %25.15e\n", ao_value[26][224]);
+printf("\n");
+
+assert( fabs(ao_value[26][219] - (  1.020298798341620e-08)) < 1.e-14 );
+assert( fabs(ao_value[26][220] - (  1.516643537739178e-08)) < 1.e-14 );
+assert( fabs(ao_value[26][221] - ( -4.686370882518819e-09)) < 1.e-14 );
+assert( fabs(ao_value[26][222] - (  7.514816980753531e-09)) < 1.e-14 );
+assert( fabs(ao_value[26][223] - ( -4.021908374204471e-09)) < 1.e-14 );
+assert( fabs(ao_value[26][224] - (  7.175045873560788e-10)) < 1.e-14 );
+
+}
+
+      #+end_src
+
+** Value, gradients, Laplacian
    :PROPERTIES:
    :Name:     qmckl_compute_ao_vgl
    :CRetType: qmckl_exit_code
    :FRetType: qmckl_exit_code
    :END:
-
-** Unoptimized version
+*** Unoptimized version
      #+NAME: qmckl_ao_vgl_args_doc
     | Variable              | Type                              | In/Out | Description                                  |
     |-----------------------+-----------------------------------+--------+----------------------------------------------|
@@ -5445,7 +6323,7 @@ integer function qmckl_compute_ao_vgl_doc_f(context, &
 end function qmckl_compute_ao_vgl_doc_f
      #+end_src
 
-** HPC version
+*** HPC version
      #+NAME: qmckl_ao_vgl_args_hpc_gaussian
     | Variable              | Type                           | In/Out | Description                                  |
     |-----------------------+--------------------------------+--------+----------------------------------------------|
@@ -5786,11 +6664,10 @@ qmckl_compute_ao_vgl_hpc_gaussian (
 #endif
      #+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)
+*** 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,
@@ -5809,6 +6686,7 @@ qmckl_compute_ao_vgl_hpc_gaussian (
           const double* shell_vgl,
           double* const ao_vgl );
      #+end_src
+
      #+begin_src c :tangle (eval h_private_func) :comments org
 #ifdef HAVE_HPC
     qmckl_exit_code qmckl_compute_ao_vgl_hpc_gaussian (
@@ -5894,8 +6772,7 @@ qmckl_compute_ao_vgl_hpc_gaussian (
     end function qmckl_compute_ao_vgl_doc
      #+end_src
 
-
-*** Provide                                                        :noexport:
+**** 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);
@@ -6036,7 +6913,7 @@ qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context)
 }
      #+end_src
 
-*** Test                                                           :noexport:
+**** Test                                                          :noexport:
 
      #+begin_src python :results output :exports none
 import numpy as np

org/qmckl_mo.org

@@ -92,10 +92,12 @@ int main() {
 
   Computed data:
 
-  |---------------+--------------------------+-------------------------------------------------------------------------------------|
+  |-----------------+--------------------------+-------------------------------------------------------------------------------------|
+  | ~mo_value~      | ~[point_num][mo_num]~    | Value of the MOs at point positions                                                 |
+  | ~mo_value_date~ | ~uint64_t~               | Late modification date of the value of the MOs at point positions                   |
   | ~mo_vgl~        | ~[point_num][5][mo_num]~ | Value, gradients, Laplacian of the MOs at point positions                           |
   | ~mo_vgl_date~   | ~uint64_t~               | Late modification date of Value, gradients, Laplacian of the MOs at point positions |
-  |---------------+--------------------------+-------------------------------------------------------------------------------------|
+  |-----------------+--------------------------+-------------------------------------------------------------------------------------|
 
 ** Data structure
 
@@ -106,7 +108,9 @@ typedef struct qmckl_mo_basis_struct {
   double * restrict coefficient_t;
 
   double * restrict mo_vgl;
+  double * restrict mo_value;
   uint64_t  mo_vgl_date;
+  uint64_t  mo_value_date;
 
   int32_t   uninitialized;
   bool      provided;
@@ -418,7 +422,464 @@ qmckl_exit_code qmckl_finalize_mo_basis(qmckl_context context) {
 
 * Computation
 
-** Computation of MOs
+** Computation of MOs: values only
+
+*** Get
+
+    #+begin_src c :comments org :tangle (eval h_func) :noweb yes
+qmckl_exit_code
+qmckl_get_mo_basis_mo_value(qmckl_context context,
+                            double* const mo_value,
+                            const int64_t size_max);
+    #+end_src
+
+    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
+qmckl_exit_code
+qmckl_get_mo_basis_mo_value(qmckl_context context,
+                            double* const mo_value,
+                            const int64_t size_max)
+{
+
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return QMCKL_NULL_CONTEXT;
+  }
+
+  qmckl_exit_code rc;
+
+  rc = qmckl_provide_ao_value(context);
+  if (rc != QMCKL_SUCCESS) return rc;
+
+  rc = qmckl_provide_mo_value(context);
+  if (rc != QMCKL_SUCCESS) return rc;
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
+  assert (ctx != NULL);
+
+  const int64_t sze = ctx->point.num * ctx->mo_basis.mo_num;
+  if (size_max < sze) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_3,
+                           "qmckl_get_mo_basis_mo_value",
+                           "input array too small");
+  }
+  memcpy(mo_value, ctx->mo_basis.mo_value, 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_mo_basis_mo_value (context, &
+          mo_value, size_max) bind(C)
+       use, intrinsic :: iso_c_binding
+       import
+       implicit none
+
+       integer (c_int64_t) , intent(in)  , value :: context
+       double precision,     intent(out)         :: mo_value(*)
+       integer (c_int64_t) , intent(in)  , value :: size_max
+     end function qmckl_get_mo_basis_mo_value
+  end interface
+    #+end_src
+
+    Uses the given array to compute the values.
+
+    #+begin_src c :comments org :tangle (eval h_func) :noweb yes
+qmckl_exit_code
+qmckl_get_mo_basis_mo_value_inplace (qmckl_context context,
+                                     double* const mo_value,
+                                     const int64_t size_max);
+    #+end_src
+
+    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
+qmckl_exit_code
+qmckl_get_mo_basis_mo_value_inplace (qmckl_context context,
+                                     double* const mo_value,
+                                     const int64_t size_max)
+{
+
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_CONTEXT,
+                           "qmckl_get_mo_basis_mo_value",
+                           NULL);
+  }
+
+  qmckl_exit_code rc;
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
+  assert (ctx != NULL);
+
+  const int64_t sze = ctx->mo_basis.mo_num * ctx->point.num;
+  if (size_max < sze) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_3,
+                           "qmckl_get_mo_basis_mo_value",
+                           "input array too small");
+  }
+
+  rc = qmckl_context_touch(context);
+  if (rc != QMCKL_SUCCESS) return rc;
+
+  double* old_array = ctx->mo_basis.mo_value;
+
+  ctx->mo_basis.mo_value = mo_value;
+
+  rc = qmckl_provide_mo_value(context);
+  if (rc != QMCKL_SUCCESS) return rc;
+
+  ctx->mo_basis.mo_value = 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_mo_basis_mo_value_inplace (context, &
+          mo_value, size_max) bind(C)
+       use, intrinsic :: iso_c_binding
+       import
+       implicit none
+       integer (c_int64_t) , intent(in)  , value :: context
+       double precision,     intent(out)         :: mo_value(*)
+       integer (c_int64_t) , intent(in)  , value :: size_max
+     end function qmckl_get_mo_basis_mo_value_inplace
+  end interface
+    #+end_src
+
+*** Provide
+
+    #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
+qmckl_exit_code qmckl_provide_mo_value(qmckl_context context);
+    #+end_src
+
+    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
+qmckl_exit_code qmckl_provide_mo_value(qmckl_context context)
+{
+
+  qmckl_exit_code rc;
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return QMCKL_NULL_CONTEXT;
+  }
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
+  assert (ctx != NULL);
+
+  if (!ctx->ao_basis.provided) {
+    return qmckl_failwith( context,
+                           QMCKL_NOT_PROVIDED,
+                           "qmckl_ao_basis",
+                           NULL);
+  }
+
+  rc = qmckl_provide_ao_value(context);
+  if (rc != QMCKL_SUCCESS) {
+    return qmckl_failwith( context,
+                           QMCKL_NOT_PROVIDED,
+                           "qmckl_ao_value",
+                           NULL);
+  }
+
+  if (!ctx->mo_basis.provided) {
+    return qmckl_failwith( context,
+                           QMCKL_NOT_PROVIDED,
+                           "qmckl_mo_basis",
+                           NULL);
+  }
+
+ /* Compute if necessary */
+  if (ctx->point.date > ctx->mo_basis.mo_value_date) {
+
+ /* Allocate array */
+    if (ctx->mo_basis.mo_value == NULL) {
+
+      qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
+      mem_info.size = ctx->point.num * ctx->mo_basis.mo_num * sizeof(double);
+      double* mo_value = (double*) qmckl_malloc(context, mem_info);
+
+      if (mo_value == NULL) {
+        return qmckl_failwith( context,
+                               QMCKL_ALLOCATION_FAILED,
+                               "qmckl_mo_basis_mo_value",
+                               NULL);
+      }
+      ctx->mo_basis.mo_value = mo_value;
+    }
+
+    if (ctx->mo_basis.mo_vgl_date == ctx->point.date) {
+
+      // mo_vgl has been computed at this step: Just copy the data.
+      
+      double * v = &(ctx->mo_basis.mo_value[0]);
+      double * vgl = &(ctx->mo_basis.mo_vgl[0]);
+      for (int i=0 ; i<ctx->point.num ; ++i) {
+        for (int k=0 ; k<ctx->mo_basis.mo_num ; ++k) {
+          v[k] = vgl[k];
+        }
+        v   += ctx->mo_basis.mo_num;
+        vgl += ctx->mo_basis.mo_num * 5;
+      }
+
+    } else {
+
+      rc = qmckl_compute_mo_basis_mo_value(context,
+                                           ctx->ao_basis.ao_num,
+                                           ctx->mo_basis.mo_num,
+                                           ctx->point.num,
+                                           ctx->mo_basis.coefficient_t,
+                                           ctx->ao_basis.ao_value,
+                                           ctx->mo_basis.mo_value);
+
+      if (rc != QMCKL_SUCCESS) {
+        return rc;
+      }
+
+    }
+
+    ctx->mo_basis.mo_value_date = ctx->date;
+  }
+
+  return QMCKL_SUCCESS;
+}
+    #+end_src
+
+*** Compute
+   :PROPERTIES:
+   :Name:     qmckl_compute_mo_basis_mo_value
+   :CRetType: qmckl_exit_code
+   :FRetType: qmckl_exit_code
+   :END:
+
+    #+NAME: qmckl_mo_basis_mo_value_args
+    | Variable            | Type                        | In/Out | Description                                     |
+    |---------------------+-----------------------------+--------+-------------------------------------------------|
+    | ~context~           | ~qmckl_context~             | in     | Global state                                    |
+    | ~ao_num~            | ~int64_t~                   | in     | Number of AOs                                   |
+    | ~mo_num~            | ~int64_t~                   | in     | Number of MOs                                   |
+    | ~point_num~         | ~int64_t~                   | in     | Number of points                                |
+    | ~coef_normalized_t~ | ~double[mo_num][ao_num]~    | in     | Transpose of the AO to MO transformation matrix |
+    | ~ao_value~          | ~double[point_num][ao_num]~ | in     | Value of the AOs                                |
+    | ~mo_value~          | ~double[point_num][mo_num]~ | out    | Value of the MOs                                |
+
+
+    The matrix of AO values is very sparse, so we use a sparse-dense
+    matrix multiplication instead of a dgemm, as exposed in
+    https://dx.doi.org/10.1007/978-3-642-38718-0_14.
+
+
+
+    #+begin_src f90 :comments org :tangle (eval f) :noweb yes
+integer function qmckl_compute_mo_basis_mo_value_doc_f(context, &
+     ao_num, mo_num, point_num, &
+     coef_normalized_t, ao_value, mo_value) &
+     result(info)
+  use qmckl
+  implicit none
+  integer(qmckl_context), intent(in)  :: context
+  integer*8             , intent(in)  :: ao_num, mo_num
+  integer*8             , intent(in)  :: point_num
+  double precision      , intent(in)  :: ao_value(ao_num,point_num)
+  double precision      , intent(in)  :: coef_normalized_t(mo_num,ao_num)
+  double precision      , intent(out) :: mo_value(mo_num,point_num)
+  integer*8 :: i,j,k
+  double precision :: c1, c2, c3, c4, c5
+
+  integer*8 :: LDA, LDB, LDC
+
+  info = QMCKL_SUCCESS
+  if (.True.)  then    ! fast algorithm
+     do j=1,point_num
+        mo_value(:,j) = 0.d0
+        do k=1,ao_num
+           if (ao_value(k,j) /= 0.d0) then
+              c1 = ao_value(k,j)
+              do i=1,mo_num
+                 mo_value(i,j) = mo_value(i,j) + coef_normalized_t(i,k) * c1
+              end do
+           end if
+        end do
+     end do
+     
+  else ! dgemm
+
+    LDA = size(coef_normalized_t,1)
+    LDB = size(ao_value,1) 
+    LDC = size(mo_value,1)
+
+    info = qmckl_dgemm(context,'N', 'N', mo_num, point_num, ao_num, 1.d0,     &
+                                    coef_normalized_t, LDA, ao_value, LDB, &
+                                    0.d0, mo_value, LDC)
+
+  end if
+
+end function qmckl_compute_mo_basis_mo_value_doc_f
+    #+end_src
+
+    #+CALL: generate_c_header(table=qmckl_mo_basis_mo_value_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_value"))
+
+   #+RESULTS:
+   #+begin_src c :tangle (eval h_func) :comments org
+   qmckl_exit_code qmckl_compute_mo_basis_mo_value (
+         const qmckl_context context,
+         const int64_t ao_num,
+         const int64_t mo_num,
+         const int64_t point_num,
+         const double* coef_normalized_t,
+         const double* ao_value,
+         double* const mo_value );
+   #+end_src
+
+   #+CALL: generate_c_header(table=qmckl_mo_basis_mo_value_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_value_doc"))
+
+   #+RESULTS:
+   #+begin_src c :tangle (eval h_func) :comments org
+   qmckl_exit_code qmckl_compute_mo_basis_mo_value_doc (
+         const qmckl_context context,
+         const int64_t ao_num,
+         const int64_t mo_num,
+         const int64_t point_num,
+         const double* coef_normalized_t,
+         const double* ao_value,
+         double* const mo_value );
+   #+end_src
+
+   #+CALL: generate_c_interface(table=qmckl_mo_basis_mo_value_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_value_doc"))
+
+    #+RESULTS:
+    #+begin_src f90 :tangle (eval f) :comments org :exports none
+    integer(c_int32_t) function qmckl_compute_mo_basis_mo_value_doc &
+        (context, ao_num, mo_num, point_num, coef_normalized_t, ao_value, mo_value) &
+        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 :: mo_num
+      integer (c_int64_t) , intent(in)  , value :: point_num
+      real    (c_double ) , intent(in)          :: coef_normalized_t(ao_num,mo_num)
+      real    (c_double ) , intent(in)          :: ao_value(ao_num,point_num)
+      real    (c_double ) , intent(out)         :: mo_value(mo_num,point_num)
+
+      integer(c_int32_t), external :: qmckl_compute_mo_basis_mo_value_doc_f
+      info = qmckl_compute_mo_basis_mo_value_doc_f &
+             (context, ao_num, mo_num, point_num, coef_normalized_t, ao_value, mo_value)
+
+    end function qmckl_compute_mo_basis_mo_value_doc
+    #+end_src
+
+    #+begin_src c :tangle (eval c) :comments org
+qmckl_exit_code
+qmckl_compute_mo_basis_mo_value (const qmckl_context context,
+                                 const int64_t ao_num,
+                                 const int64_t mo_num,
+                                 const int64_t point_num,
+                                 const double* coef_normalized_t,
+                                 const double* ao_value,
+                                 double* const mo_value )
+{
+#ifdef HAVE_HPC
+  return qmckl_compute_mo_basis_mo_value_hpc (context, ao_num, mo_num, point_num, coef_normalized_t, ao_value, mo_value);
+#else
+  return qmckl_compute_mo_basis_mo_value_doc (context, ao_num, mo_num, point_num, coef_normalized_t, ao_value, mo_value);
+#endif
+}
+    #+end_src
+
+*** HPC version
+
+
+    #+begin_src c :tangle (eval h_func) :comments org
+#ifdef HAVE_HPC
+qmckl_exit_code
+qmckl_compute_mo_basis_mo_value_hpc (const qmckl_context context,
+                                     const int64_t ao_num,
+                                     const int64_t mo_num,
+                                     const int64_t point_num,
+                                     const double* coef_normalized_t,
+                                     const double* ao_value,
+                                     double* const mo_value );
+#endif
+    #+end_src
+
+    #+begin_src c :tangle (eval c) :comments org
+#ifdef HAVE_HPC
+qmckl_exit_code
+qmckl_compute_mo_basis_mo_value_hpc (const qmckl_context context,
+                                     const int64_t ao_num,
+                                     const int64_t mo_num,
+                                     const int64_t point_num,
+                                     const double* restrict coef_normalized_t,
+                                     const double* restrict ao_value,
+                                     double* restrict const mo_value )
+{
+  assert (context != QMCKL_NULL_CONTEXT);
+
+#ifdef HAVE_OPENMP
+  #pragma omp parallel for
+#endif
+  for (int64_t ipoint=0 ; ipoint < point_num ; ++ipoint) {
+    double* restrict const vgl1  = &(mo_value[ipoint*mo_num]);
+    const double* restrict avgl1 = &(ao_value[ipoint*ao_num]);
+
+    for (int64_t i=0 ; i<mo_num ; ++i) {
+      vgl1[i] = 0.;
+    }
+
+    int64_t nidx=0;
+    int64_t idx[ao_num];
+    double  av1[ao_num];
+    for (int64_t k=0 ; k<ao_num ; ++k) {
+      if (avgl1[k] != 0.) {
+        idx[nidx] = k;
+        av1[nidx] = avgl1[k];
+        ++nidx;
+      }
+    }
+
+    int64_t n;
+
+    for (n=0 ; n < nidx-4 ; n+=4) {
+      const double* restrict ck1 = coef_normalized_t + idx[n  ]*mo_num;
+      const double* restrict ck2 = coef_normalized_t + idx[n+1]*mo_num;
+      const double* restrict ck3 = coef_normalized_t + idx[n+2]*mo_num;
+      const double* restrict ck4 = coef_normalized_t + idx[n+3]*mo_num;
+
+      const double a11 = av1[n  ];
+      const double a21 = av1[n+1];
+      const double a31 = av1[n+2];
+      const double a41 = av1[n+3];
+
+#ifdef HAVE_OPENMP
+#pragma omp simd
+#endif
+      for (int64_t i=0 ; i<mo_num ; ++i) {
+        vgl1[i] = vgl1[i] + ck1[i] * a11 + ck2[i] * a21 + ck3[i] * a31 + ck4[i] * a41;
+      }
+    }
+
+    const int64_t n0 = n < 0 ? 0 : n;
+    for (int64_t m=n0 ; m < nidx ; m+=1) {
+      const double* restrict ck = coef_normalized_t + idx[m]*mo_num;
+      const double a1 = av1[m];
+
+#ifdef HAVE_OPENMP
+  #pragma omp simd
+#endif
+      for (int64_t i=0 ; i<mo_num ; ++i) {
+        vgl1[i] += ck[i] * a1;
+      }
+    }
+  }
+  return QMCKL_SUCCESS;
+}
+#endif
+    #+end_src
+
+** Computation of MOs: values, gradient, Laplacian
 
 *** Get
 
@@ -774,7 +1235,6 @@ qmckl_compute_mo_basis_mo_vgl (const qmckl_context context,
 }
     #+end_src
 
-
 *** HPC version
 
 
@@ -918,7 +1378,7 @@ qmckl_compute_mo_basis_mo_vgl_hpc (const qmckl_context context,
 #endif
     #+end_src
 
-*** Test
+** Test
 
    #+begin_src python :results output :exports none
 import numpy as np
@@ -1107,10 +1567,36 @@ assert (rc == QMCKL_SUCCESS);
 
 assert(qmckl_mo_basis_provided(context));
 
+rc = qmckl_context_touch(context);
+assert (rc == QMCKL_SUCCESS);
+
+double mo_value[walk_num*elec_num][chbrclf_mo_num];
+rc = qmckl_get_mo_basis_mo_value(context, &(mo_value[0][0]), walk_num*elec_num*chbrclf_mo_num);
+assert (rc == QMCKL_SUCCESS);
+
 double mo_vgl[walk_num*elec_num][5][chbrclf_mo_num];
 rc = qmckl_get_mo_basis_mo_vgl(context, &(mo_vgl[0][0][0]), walk_num*elec_num*5*chbrclf_mo_num);
 assert (rc == QMCKL_SUCCESS);
 
+for (int i=0 ; i< walk_num*elec_num; ++i) {
+  for (int k=0 ; k< chbrclf_mo_num ; ++k) {
+    assert(fabs(mo_vgl[i][0][k] - mo_value[i][k]) < 1.e-12) ;
+  }
+ }
+
+rc = qmckl_context_touch(context);
+assert (rc == QMCKL_SUCCESS);
+ 
+rc = qmckl_get_mo_basis_mo_value(context, &(mo_value[0][0]), walk_num*elec_num*chbrclf_mo_num);
+assert (rc == QMCKL_SUCCESS);
+
+for (int i=0 ; i< walk_num*elec_num; ++i) {
+  for (int k=0 ; k< chbrclf_mo_num ; ++k) {
+    assert(fabs(mo_vgl[i][0][k] - mo_value[i][k]) < 1.e-12) ;
+  }
+ }
+
+
 // Test overlap of MO
 //double point_x[10];
 //double point_y[10];