From 21336e0178cfaa8cc36580cfbc38421bc276c48d Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Tue, 14 Mar 2023 14:59:14 +0100
Subject: [PATCH] Working on e-n cusp

---
 org/qmckl_ao.org           |  21 +-
 org/qmckl_determinant.org  |   2 +-
 org/qmckl_local_energy.org |   2 +-
 org/qmckl_mo.org           | 914 ++++++++++++++++++++++++++++++++++---
 org/qmckl_trexio.org       |   2 +-
 5 files changed, 864 insertions(+), 77 deletions(-)

diff --git a/org/qmckl_ao.org b/org/qmckl_ao.org
index 8099ff9..e56ba7f 100644
--- a/org/qmckl_ao.org
+++ b/org/qmckl_ao.org
@@ -2587,7 +2587,7 @@ qmckl_exit_code qmckl_finalize_basis(qmckl_context context) {
   /* ao_ang_mom */
   {
     qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
-    mem_info.size = ctx->ao_basis.ao_num * sizeof(int64_t);
+    mem_info.size = ctx->ao_basis.ao_num * sizeof(int32_t);
 
     ctx->ao_basis.ao_ang_mom = (int32_t*) qmckl_malloc(context, mem_info);
 
@@ -2598,36 +2598,39 @@ qmckl_exit_code qmckl_finalize_basis(qmckl_context context) {
                              NULL);
     }
 
+    mem_info = qmckl_memory_info_struct_zero;
+    mem_info.size = ctx->ao_basis.ao_num * sizeof(int64_t);
     ctx->ao_basis.ao_nucl = (int64_t*) qmckl_malloc(context, mem_info);
 
-    if (ctx->ao_basis.ao_ang_mom == NULL) {
+    if (ctx->ao_basis.ao_nucl == NULL) {
       return qmckl_failwith( context,
                              QMCKL_ALLOCATION_FAILED,
                              "ao_basis.ao_nucl",
                              NULL);
     }
 
+    int64_t lstart[32];
+    for (int32_t l=0 ; l<32 ; ++l) {
+      lstart[l] = l*(l+1)*(l+2)/6;
+    }
+    
     int64_t ao_idx = 0;
     for (int64_t inucl=0 ; inucl<nucl_num ; ++inucl) {
       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 int l = ctx->ao_basis.shell_ang_mom[ishell];
-        const int mmax = l*(l+1)*(l+2)/6;
-        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) {
-          for (int m=0 ; m < mmax ; m++) {
+        assert (l<32);
+          for (int m=lstart[l] ; m < lstart[l+1]; m++) {
             ctx->ao_basis.ao_ang_mom[ao_idx] = l;
             ctx->ao_basis.ao_nucl[ao_idx] = inucl;
             ++ao_idx;
           }
-        }
       }
     }
+    assert( ao_idx == ctx->ao_basis.ao_num );
   }
 
-
   /* Normalize coefficients */
   {
     qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
diff --git a/org/qmckl_determinant.org b/org/qmckl_determinant.org
index 2f20c8a..3f65628 100644
--- a/org/qmckl_determinant.org
+++ b/org/qmckl_determinant.org
@@ -1225,7 +1225,7 @@ qmckl_check(context, rc);
 
 const double  * mo_coefficient          =  &(chbrclf_mo_coef[0]);
 
-rc = qmckl_set_mo_basis_coefficient(context, mo_coefficient);
+rc = qmckl_set_mo_basis_coefficient(context, mo_coefficient, chbrclf_mo_num*chbrclf_ao_num);
 qmckl_check(context, rc);
 
 assert(qmckl_mo_basis_provided(context));
diff --git a/org/qmckl_local_energy.org b/org/qmckl_local_energy.org
index 5ddeb39..ad583a7 100644
--- a/org/qmckl_local_energy.org
+++ b/org/qmckl_local_energy.org
@@ -697,7 +697,7 @@ assert (rc == QMCKL_SUCCESS);
 
 const double  * mo_coefficient          =  &(chbrclf_mo_coef[0]);
 
-rc = qmckl_set_mo_basis_coefficient(context, mo_coefficient);
+rc = qmckl_set_mo_basis_coefficient(context, mo_coefficient, chbrclf_ao_num*chbrclf_mo_num);
 assert (rc == QMCKL_SUCCESS);
 
 assert(qmckl_mo_basis_provided(context));
diff --git a/org/qmckl_mo.org b/org/qmckl_mo.org
index 417c170..e2df46a 100644
--- a/org/qmckl_mo.org
+++ b/org/qmckl_mo.org
@@ -14,6 +14,53 @@ By default, all the MOs are computed. A subset of MOs can be selected
 for computation, for example to remove computation of the useless
 virtual MOs present in a MO coefficient matrix.
 
+For each nucleus, a radius $r_{\text{cusp}}$ can be given such that if
+the electron-nucleus distance is smaller than $r_{\text{cusp}}$, all
+the MOs are replaced by functions with the correct cusp condition and
+such that the values and the gradients of the MOs are continuous at
+$r_{\text{cusp}}$.
+
+
+
+Molecular orbitals (MOs) are defined in the basis of atomic orbitals
+(AOs) using a coefficient matrix $C$, which determines how the AOs are
+combined to form the MOs.
+
+The equation for the MOs $\phi_i(\mathbf{r})$ is given by:
+
+$$\phi_i(\mathbf{r}) = \sum_{\mu} C_{\mu i} \chi_{\mu}(\mathbf{r})$$
+
+where $i$ labels the MO, $\mu$ labels the AO, $C_{\mu i}$ is the
+coefficient of AO $\mu$ in MO $i$, and $\chi_{\mu}(\mathbf{r})$ is the
+AO itself.
+
+In some cases, it may be desirable to only compute a subset of the
+MOs, for example to exclude virtual MOs that do not contribute
+to the wave function. This can be achieved by selecting a subset of
+columns from the coefficient matrix $C$.
+
+The exact wave function must have a cusp at the positions of the
+nuclei to ensure that the kinetic energy diverges and cancels out the
+divergence of the potential, resulting in a finite total energy.
+To ensure that the cusp condition is satisfied, a modification can be
+made to the molecular orbitals (MOs) when the electron-nucleus
+distance is smaller than a certain radius $r_{\text{cusp}}$. At
+distances smaller than $r_{\text{cusp}}$, the MOs are replaced by
+functions that have the correct electron-nucleus cusp condition and
+that ensure that the values and the gradients of the MOs are
+continuous at $r_{\text{cusp}}$. 
+
+A radius $r_{\text{cusp}\, A}$ is given for each nucleus $A$, default is
+zero. If an electron is closer to the nucleus $A$ than
+$r_{\text{cusp}\, A}$, the MOs are locally replaced as
+\[
+\phi_{\text{cusp}\, i}(\mathbf{r}) = \phi_i(\mathbf{r}) - \phi_{s_A i}(\mathbf{r}) + \sum_{k=0}^{3} f_k\, |\mathbf{r}-\mathbf{R}_A|^k
+\]
+where $\phi_{s_A i}$ are the contributions of the $s$ AOs centered at
+$A$ to MO $i$.
+The coefficients $f_k$ are such that the value and gradient of the
+MO are continuous at $r_{\text{cusp}}$, and the electron-nucleus
+cusp is exact.
 
 * Headers                                                          :noexport:
   #+begin_src elisp :noexport :results none
@@ -43,6 +90,7 @@ virtual MOs present in a MO coefficient matrix.
 #include <stdio.h>
 #include <math.h>
 #include "chbrclf.h"
+#include "qmckl_electron_private_func.h"
 #include "qmckl_ao_private_func.h"
 #include "qmckl_mo_private_func.h"
 
@@ -83,21 +131,20 @@ int main() {
 
   The following arrays are stored in the context:
 
-
-  |-------------------+--------------------+----------------------------------------|
-  | ~mo_num~          |                    | Number of MOs                          |
-  | ~coefficient~     | ~[mo_num][ao_num]~ | MO coefficients                        |
-  | ~coefficient_t~   | ~[ao_num][mo_num]~ | Transposed of the Orbital coefficients |
-  |-------------------+--------------------+----------------------------------------|
+  |-----------------+--------------------+----------------------------------------------|
+  | ~mo_num~        |                    | Number of MOs                                |
+  | ~coefficient~   | ~[mo_num][ao_num]~ | MO coefficients                              |
+  | ~coefficient_t~ | ~[ao_num][mo_num]~ | Transposed of the Orbital coefficients       |
+  | ~r_cusp~        | ~[nucl_num]~       | Radius of the functions for Cusp adjustments |
+  |-----------------+--------------------+----------------------------------------------|
 
   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 |
-  |-----------------+--------------------------+-------------------------------------------------------------------------------------|
+  |-----------------+--------------------------+----------------------------------------------------------------------------------|
+  | ~r_cusp_param~  | ~[nucl_num][mo_num][4]~  | Parameters of the functions for Cusp adjustments                                 |
+  | ~mo_value~      | ~[point_num][mo_num]~    | Value of the MOs at point positions                                              |
+  | ~mo_vgl~        | ~[point_num][5][mo_num]~ | Value, gradients, Laplacian of the MOs at point positions                        |
+  |-----------------+--------------------------+----------------------------------------------------------------------------------|
 
 ** Data structure
 
@@ -106,9 +153,12 @@ typedef struct qmckl_mo_basis_struct {
   int64_t  mo_num;
   double * restrict coefficient;
   double * restrict coefficient_t;
+  double * restrict r_cusp;
 
   double * restrict mo_vgl;
   double * restrict mo_value;
+  double * restrict r_cusp_param;
+
   uint64_t  mo_vgl_date;
   uint64_t  mo_value_date;
 
@@ -151,7 +201,8 @@ qmckl_exit_code qmckl_init_mo_basis(qmckl_context context) {
 
    #+begin_src c :comments org :tangle (eval h_func)
 qmckl_exit_code  qmckl_set_mo_basis_mo_num           (qmckl_context context, const int64_t   mo_num);
-qmckl_exit_code  qmckl_set_mo_basis_coefficient      (qmckl_context context, const double  * coefficient);
+qmckl_exit_code  qmckl_set_mo_basis_coefficient      (qmckl_context context, const double  * coefficient, const int64_t size_max);
+qmckl_exit_code  qmckl_set_mo_basis_r_cusp           (qmckl_context context, const double  * r_cusp, const int64_t size_max);
    #+end_src
 
    #+NAME:pre
@@ -198,10 +249,14 @@ qmckl_exit_code qmckl_set_mo_basis_mo_num(qmckl_context context, const int64_t m
   ctx->mo_basis.mo_num = mo_num;
 
   <<post>>
+
+  return QMCKL_SUCCESS;
 }
 
-qmckl_exit_code  qmckl_set_mo_basis_coefficient(qmckl_context context, const double* coefficient) {
-  
+qmckl_exit_code  qmckl_set_mo_basis_coefficient(qmckl_context context,
+                                                const double* coefficient,
+                                                const int64_t size_max)
+{
   int32_t mask = 1 << 1;
 
   <<pre>>
@@ -215,6 +270,12 @@ qmckl_exit_code  qmckl_set_mo_basis_coefficient(qmckl_context context, const dou
     }
   }
 
+  if (size_max < ctx->ao_basis.ao_num * ctx->mo_basis.mo_num) {
+    return qmckl_failwith( context, QMCKL_INVALID_ARG_3,
+                             "qmckl_set_mo_basis_coefficient",
+                             "Array too small");
+  }
+  
   qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
   mem_info.size = ctx->ao_basis.ao_num * ctx->mo_basis.mo_num * sizeof(double);
   double* new_array = (double*) qmckl_malloc(context, mem_info);
@@ -232,6 +293,54 @@ qmckl_exit_code  qmckl_set_mo_basis_coefficient(qmckl_context context, const dou
   <<post>>
 }
 
+
+qmckl_exit_code
+qmckl_set_mo_basis_r_cusp(qmckl_context context,
+                          const double* r_cusp,
+                          const int64_t size_max)
+{
+  
+  int32_t mask = 0;
+
+  <<pre>>
+
+  if (r_cusp == NULL) {
+      return qmckl_failwith( context, QMCKL_INVALID_ARG_2,
+                             "qmckl_set_mo_basis_r_cusp",
+                             "r_cusp: Null pointer");
+  }
+        
+  if (size_max < ctx->nucleus.num) {
+    return qmckl_failwith( context, QMCKL_INVALID_ARG_3,
+                             "qmckl_set_mo_basis_r_cusp",
+                             "Array too small");
+  }
+  
+
+  if (ctx->mo_basis.r_cusp != NULL) {
+    qmckl_exit_code rc = qmckl_free(context, ctx->mo_basis.r_cusp);
+    if (rc != QMCKL_SUCCESS) {
+      return rc;
+    }
+  }
+
+  qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
+  mem_info.size = ctx->nucleus.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_mo_basis_r_cusp",
+                           NULL);
+  }
+
+  memcpy(new_array, r_cusp, mem_info.size);
+
+  ctx->mo_basis.r_cusp = new_array;
+  return QMCKL_SUCCESS;
+
+}
+
    #+end_src
 
  When the basis set is completely entered, other data structures are
@@ -437,22 +546,39 @@ interface
     use, intrinsic :: iso_c_binding
     import
     implicit none
-    integer (c_int64_t) , intent(in)  , value :: context
-    double precision, intent(out)             :: coefficient(*)
+    integer (c_int64_t) , intent(in), value   :: context
+    double precision    , intent(out)         :: coefficient(*)
     integer (c_int64_t) , intent(in), value   :: size_max
   end function qmckl_get_mo_basis_coefficient
 end interface
 
+interface
+  integer(qmckl_exit_code) function qmckl_set_mo_basis_r_cusp(context, &
+       r_cusp, size_max) bind(C)
+    use, intrinsic :: iso_c_binding
+    import
+    implicit none
+    integer (c_int64_t) , intent(in), value   :: context
+    double precision    , intent(in)          :: r_cusp(*)
+    integer (c_int64_t) , intent(in), value   :: size_max
+  end function qmckl_set_mo_basis_r_cusp
+end interface
+
     #+end_src
 
 ** Update
 
-   Useless MOs can be removed, for instance virtual MOs in a single
-   determinant calculation.
+   It may be desirable to remove certain molecular orbitals (MOs) that
+   do not significantly contribute to the wave function.  In
+   particular, in a single determinant calculation, the virtual MOs
+   can be removed as they do not participate in the ground state
+   configuration.
 
-   To select a subset of MOs that will be kept, create an array of
-   integers of size =mo_num=. If the integer is zero, the MO is dropped,
-   otherwise it is kept.
+   To select a subset of MOs that will be kept, an array of integers of
+   size ~mo_num~ can be created. If the integer corresponding to an MO is
+   zero, that MO is dropped and will not be included in the
+   calculation. If the integer is non-zero, the MO will be kept.
+   
 
    #+begin_src c :comments org :tangle (eval h_func)
 qmckl_exit_code
@@ -548,6 +674,7 @@ end interface
 
 * Computation
 
+** Parameters of the cusp-correction functions
 ** Computation of MOs: values only
 
 *** Get
@@ -769,14 +896,37 @@ qmckl_exit_code qmckl_provide_mo_basis_mo_value(qmckl_context context)
                                NULL);
       }
 
-      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 (ctx->mo_basis.r_cusp == NULL) {
+        /* No cusp correction */
+        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);
+      } else {
+        rc = qmckl_provide_en_distance(context);
+        if (rc != QMCKL_SUCCESS) {
+          return qmckl_failwith( context,
+                                 QMCKL_NOT_PROVIDED,
+                                 "qmckl_electron_en_distance",
+                                 NULL);
+        }
 
+        rc = qmckl_compute_mo_basis_mo_value_cusp(context,
+                                                  ctx->nucleus.num,
+                                                  ctx->ao_basis.ao_num,
+                                                  ctx->mo_basis.mo_num,
+                                                  ctx->point.num,
+                                                  ctx->ao_basis.ao_nucl, 
+                                                  ctx->ao_basis.ao_ang_mom,
+                                                  ctx->electron.en_distance,
+                                                  ctx->mo_basis.r_cusp,
+                                                  ctx->mo_basis.coefficient_t,
+                                                  ctx->ao_basis.ao_value,
+                                                  ctx->mo_basis.mo_value);
+      }
     }
     #+end_src
 
@@ -835,37 +985,18 @@ integer function qmckl_compute_mo_basis_mo_value_doc_f(context, &
   double precision      , intent(in)  :: ao_value(ao_num,point_num)
   double precision      , intent(in)  :: coefficient_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
+  integer*8 :: j,k
 
   info = QMCKL_SUCCESS
-  if (.True.)  then    ! fast algorithm
-     do j=1,point_num
-        mo_value(:,j) = 0.d0
-        do k=1,ao_num
-           c1 = ao_value(k,j)
-           if (c1 /= 0.d0) then
-              do i=1,mo_num
-                 mo_value(i,j) = mo_value(i,j) + coefficient_t(i,k) * c1
-              end do
-           end if
-        end do
+
+  do j=1,point_num
+     mo_value(:,j) = 0.d0
+     do k=1,ao_num
+        if (ao_value(k,j) == 0.d0) cycle
+        mo_value(:,j) = mo_value(:,j) + coefficient_t(:,k) * ao_value(k,j)
      end do
+  end do
      
-  else ! dgemm for checking
-
-    LDA = size(coefficient_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,     &
-                                    coefficient_t, LDA, ao_value, LDB, &
-                                    0.d0, mo_value, LDC)
-
-  end if
-
 end function qmckl_compute_mo_basis_mo_value_doc_f
     #+end_src
 
@@ -1234,13 +1365,36 @@ qmckl_exit_code qmckl_provide_mo_basis_mo_vgl(qmckl_context context)
                              NULL);
     }
 
-    rc = qmckl_compute_mo_basis_mo_vgl(context,
-                                       ctx->ao_basis.ao_num,
-                                       ctx->mo_basis.mo_num,
-                                       ctx->point.num,
-                                       ctx->mo_basis.coefficient_t,
-                                       ctx->ao_basis.ao_vgl,
-                                       ctx->mo_basis.mo_vgl);
+    if (ctx->mo_basis.r_cusp == NULL) {
+      /* No cusp correction */
+      rc = qmckl_compute_mo_basis_mo_vgl(context,
+                                         ctx->ao_basis.ao_num,
+                                         ctx->mo_basis.mo_num,
+                                         ctx->point.num,
+                                         ctx->mo_basis.coefficient_t,
+                                         ctx->ao_basis.ao_vgl,
+                                         ctx->mo_basis.mo_vgl);
+    } else {
+      rc = qmckl_provide_en_distance(context);
+      if (rc != QMCKL_SUCCESS) {
+        return qmckl_failwith( context,
+                               QMCKL_NOT_PROVIDED,
+                               "qmckl_electron_en_distance",
+                               NULL);
+      }
+      rc = qmckl_compute_mo_basis_mo_vgl_cusp(context,
+                                              ctx->nucleus.num,
+                                              ctx->ao_basis.ao_num,
+                                              ctx->mo_basis.mo_num,
+                                              ctx->point.num,
+                                              ctx->ao_basis.ao_nucl, 
+                                              ctx->ao_basis.ao_ang_mom,
+                                              ctx->electron.en_distance,
+                                              ctx->mo_basis.r_cusp,
+                                              ctx->mo_basis.coefficient_t,
+                                              ctx->ao_basis.ao_vgl,
+                                              ctx->mo_basis.mo_vgl);
+    }
     #+end_src
 
 #+CALL: write_provider_post( group="mo_basis", data="mo_vgl" )
@@ -1299,6 +1453,8 @@ integer function qmckl_compute_mo_basis_mo_vgl_doc_f(context, &
   integer*8 :: i,j,k
   double precision :: c1, c2, c3, c4, c5
 
+  info = QMCKL_SUCCESS
+
   do j=1,point_num
      mo_vgl(:,:,j) = 0.d0
      do k=1,ao_num
@@ -1318,7 +1474,6 @@ integer function qmckl_compute_mo_basis_mo_vgl_doc_f(context, &
         end if
      end do
   end do
-  info = QMCKL_SUCCESS
 
 ! info = qmckl_dgemm(context,'N', 'N', mo_num, point_num, ao_num, 1.d0, &
 !      coefficient_t, int(size(coefficient_t,1),8),      &
@@ -1542,6 +1697,635 @@ qmckl_compute_mo_basis_mo_vgl_hpc (const qmckl_context context,
 #endif
     #+end_src
 
+** Computation of cusp-corrected MOs: values only
+
+*** Compute
+   :PROPERTIES:
+   :Name:     qmckl_compute_mo_basis_mo_value_cusp
+   :CRetType: qmckl_exit_code
+   :FRetType: qmckl_exit_code
+   :END:
+
+    #+NAME: qmckl_mo_basis_mo_value_cusp_args
+    | Variable        | Type                          | In/Out | Description                                     |
+    |-----------------+-------------------------------+--------+-------------------------------------------------|
+    | ~context~       | ~qmckl_context~               | in     | Global state                                    |
+    | ~nucl_num~      | ~int64_t~                     | in     | Number of nuclei                                |
+    | ~ao_num~        | ~int64_t~                     | in     | Number of AOs                                   |
+    | ~mo_num~        | ~int64_t~                     | in     | Number of MOs                                   |
+    | ~point_num~     | ~int64_t~                     | in     | Number of points                                |
+    | ~ao_nucl~       | ~int64_t[ao_num]~             | in     | Nucleus on which the AO is centered             |
+    | ~ao_ang_mom~    | ~int32_t[ao_num]~             | in     | Angular momentum of the shell                   |
+    | ~en_distance~   | ~double[point_num][nucl_num]~ | in     | Electron-nucleus distances                      |
+    | ~r_cusp~        | ~double[nucl_num]~            | in     | Cusp-adjustment radius                          |
+    | ~coefficient_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    | Cusp correction for the values of the MOs       |
+
+
+
+
+    #+begin_src f90 :comments org :tangle (eval f) :noweb yes
+integer function qmckl_compute_mo_basis_mo_value_cusp_doc_f(context, &
+     nucl_num, ao_num, mo_num, point_num, ao_nucl, ao_ang_mom, en_distance, &
+     r_cusp, coefficient_t, ao_value, mo_value) &
+     result(info)
+  use qmckl
+  implicit none
+  integer(qmckl_context), intent(in)  :: context
+  integer*8             , intent(in)  :: nucl_num, ao_num, mo_num, point_num
+  integer*8             , intent(in)  :: ao_nucl(ao_num)
+  integer*4             , intent(in)  :: ao_ang_mom(ao_num)
+  double precision      , intent(in)  :: en_distance(nucl_num, point_num)
+  double precision      , intent(in)  :: r_cusp(nucl_num)
+  double precision      , intent(in)  :: coefficient_t(mo_num,ao_num)
+  double precision      , intent(in)  :: ao_value(ao_num,point_num)
+  double precision      , intent(out) :: mo_value(mo_num,point_num)
+
+  integer*8 :: i,j,k, inucl
+
+  info = QMCKL_SUCCESS
+
+  do j=1,point_num
+     mo_value(:,j) = 0.d0
+     do k=1,ao_num
+        if (ao_value(k,j) == 0.d0) cycle
+        inucl = ao_nucl(k)+1
+        if ( (en_distance(inucl,j) < r_cusp(inucl)) .and. (ao_ang_mom(k) == 0) ) cycle
+        mo_value(:,j) = mo_value(:,j) + coefficient_t(:,k) * ao_value(k,j)
+     end do
+  end do
+     
+end function qmckl_compute_mo_basis_mo_value_cusp_doc_f
+    #+end_src
+
+    #+CALL: generate_c_header(table=qmckl_mo_basis_mo_value_cusp_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_value_cusp"))
+
+   #+RESULTS:
+   #+begin_src c :tangle (eval h_func) :comments org
+   qmckl_exit_code qmckl_compute_mo_basis_mo_value_cusp (
+         const qmckl_context context,
+         const int64_t nucl_num,
+         const int64_t ao_num,
+         const int64_t mo_num,
+         const int64_t point_num,
+         const int64_t* ao_nucl,
+         const int32_t* ao_ang_mom,
+         const double* en_distance,
+         const double* r_cusp,
+         const double* coefficient_t,
+         const double* ao_value,
+         double* const mo_value ); 
+   #+end_src
+
+   #+CALL: generate_c_header(table=qmckl_mo_basis_mo_value_cusp_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_value_cusp_doc"))
+
+   #+RESULTS:
+   #+begin_src c :tangle (eval h_func) :comments org
+   qmckl_exit_code qmckl_compute_mo_basis_mo_value_cusp_doc (
+         const qmckl_context context,
+         const int64_t nucl_num,
+         const int64_t ao_num,
+         const int64_t mo_num,
+         const int64_t point_num,
+         const int64_t* ao_nucl,
+         const int32_t* ao_ang_mom,
+         const double* en_distance,
+         const double* r_cusp,
+         const double* coefficient_t,
+         const double* ao_value,
+         double* const mo_value ); 
+   #+end_src
+
+   #+CALL: generate_c_interface(table=qmckl_mo_basis_mo_value_cusp_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_value_cusp_doc"))
+
+    #+RESULTS:
+    #+begin_src f90 :tangle (eval f) :comments org :exports none
+    integer(c_int32_t) function qmckl_compute_mo_basis_mo_value_cusp_doc &
+        (context, &
+         nucl_num, &
+         ao_num, &
+         mo_num, &
+         point_num, &
+         ao_nucl, &
+         ao_ang_mom, &
+         en_distance, &
+         r_cusp, &
+         coefficient_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 :: nucl_num
+      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
+      integer (c_int64_t) , intent(in)          :: ao_nucl(ao_num)
+      integer (c_int32_t) , intent(in)          :: ao_ang_mom(ao_num)
+      real    (c_double ) , intent(in)          :: en_distance(nucl_num,point_num)
+      real    (c_double ) , intent(in)          :: r_cusp(nucl_num)
+      real    (c_double ) , intent(in)          :: coefficient_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_cusp_doc_f
+      info = qmckl_compute_mo_basis_mo_value_cusp_doc_f &
+             (context, &
+         nucl_num, &
+         ao_num, &
+         mo_num, &
+         point_num, &
+         ao_nucl, &
+         ao_ang_mom, &
+         en_distance, &
+         r_cusp, &
+         coefficient_t, &
+         ao_value, &
+         mo_value)
+
+    end function qmckl_compute_mo_basis_mo_value_cusp_doc
+    #+end_src
+
+    #+begin_src c :tangle (eval c) :comments org
+qmckl_exit_code
+qmckl_compute_mo_basis_mo_value_cusp (const qmckl_context context,
+                                      const int64_t nucl_num,
+                                      const int64_t ao_num,
+                                      const int64_t mo_num,
+                                      const int64_t point_num,
+                                      const int64_t* ao_nucl,
+                                      const int32_t* ao_ang_mom,
+                                      const double* en_distance,
+                                      const double* r_cusp,
+                                      const double* coefficient_t,
+                                      const double* ao_value,
+                                      double* const mo_value )
+{
+#ifdef HAVE_HPC
+  return qmckl_compute_mo_basis_mo_value_cusp_hpc (context, nucl_num, ao_num, mo_num, point_num, 
+                                                   ao_nucl, ao_ang_mom, en_distance, r_cusp,
+                                                   coefficient_t, ao_value, mo_value ); 
+#else
+  return qmckl_compute_mo_basis_mo_value_cusp_doc (context, nucl_num, ao_num, mo_num, point_num, 
+                                                   ao_nucl, ao_ang_mom, en_distance, r_cusp,
+                                                   coefficient_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_cusp_hpc (const qmckl_context context,
+                                          const int64_t nucl_num,
+                                          const int64_t ao_num,
+                                          const int64_t mo_num,
+                                          const int64_t point_num,
+                                          const int64_t* ao_nucl,
+                                          const int32_t* ao_ang_mom,
+                                          const double* en_distance,
+                                          const double* r_cusp,
+                                          const double* coefficient_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_cusp_hpc (const qmckl_context context,
+                                          const int64_t nucl_num,
+                                          const int64_t ao_num,
+                                          const int64_t mo_num,
+                                          const int64_t point_num,
+                                          const int64_t* ao_nucl,
+                                          const int32_t* ao_ang_mom,
+                                          const double* en_distance,
+                                          const double* r_cusp,
+                                          const double* coefficient_t,
+                                          const double* ao_value,
+                                          double* 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]);
+    const double* restrict ria   = &(en_distance[ipoint*nucl_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.) {
+        const int64_t inucl = ao_nucl[k];
+        if (ria[inucl] > r_cusp[inucl] || ao_ang_mom[k] > 0) {
+          idx[nidx] = k;
+          av1[nidx] = avgl1[k];
+          ++nidx;
+        }
+      }
+    }
+
+    int64_t n=0;
+
+    for (n=0 ; n < nidx-4 ; n+=4) {
+      const double* restrict ck1 = coefficient_t + idx[n  ]*mo_num;
+      const double* restrict ck2 = coefficient_t + idx[n+1]*mo_num;
+      const double* restrict ck3 = coefficient_t + idx[n+2]*mo_num;
+      const double* restrict ck4 = coefficient_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;
+      }
+    }
+
+    for (int64_t m=n ; m < nidx ; m+=1) {
+      const double* restrict ck = coefficient_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 cusp-corrected MOs: values, gradient, Laplacian
+
+*** Compute
+   :PROPERTIES:
+   :Name:     qmckl_compute_mo_basis_mo_vgl
+   :CRetType: qmckl_exit_code
+   :FRetType: qmckl_exit_code
+   :END:
+
+    #+NAME: qmckl_mo_basis_mo_vgl_cusp_args
+    | Variable        | Type                           | In/Out | Description                                     |
+    |-----------------+--------------------------------+--------+-------------------------------------------------|
+    | ~context~       | ~qmckl_context~                | in     | Global state                                    |
+    | ~nucl_num~      | ~int64_t~                      | in     | Number of nuclei                                |
+    | ~ao_num~        | ~int64_t~                      | in     | Number of AOs                                   |
+    | ~mo_num~        | ~int64_t~                      | in     | Number of MOs                                   |
+    | ~point_num~     | ~int64_t~                      | in     | Number of points                                |
+    | ~ao_nucl~       | ~int64_t[ao_num]~              | in     | Nucleus on which the AO is centered             |
+    | ~ao_ang_mom~    | ~int32_t[ao_num]~              | in     | Angular momentum of the shell                   |
+    | ~en_distance~   | ~double[point_num][nucl_num]~  | in     | Electron-nucleus distances                      |
+    | ~r_cusp~        | ~double[nucl_num]~             | in     | Cusp-adjustment radius                          |
+    | ~coefficient_t~ | ~double[mo_num][ao_num]~       | in     | Transpose of the AO to MO transformation matrix |
+    | ~ao_vgl~        | ~double[point_num][5][ao_num]~ | in     | Value, gradients and Laplacian of the AOs       |
+    | ~mo_vgl~        | ~double[point_num][5][mo_num]~ | out    | Value, gradients and Laplacian of the MOs       |
+
+
+
+    #+begin_src f90 :comments org :tangle (eval f) :noweb yes
+integer function qmckl_compute_mo_basis_mo_vgl_cusp_doc_f(context, &
+     nucl_num, ao_num, mo_num, point_num, ao_nucl, ao_ang_mom, en_distance, &
+     r_cusp, coefficient_t, ao_vgl, mo_vgl) &
+     result(info)
+  use qmckl
+  implicit none
+  integer(qmckl_context), intent(in)  :: context
+  integer*8             , intent(in)  :: nucl_num, ao_num, mo_num, point_num
+  integer*8             , intent(in)  :: ao_nucl(ao_num)
+  integer*4             , intent(in)  :: ao_ang_mom(ao_num)
+  double precision      , intent(in)  :: en_distance(nucl_num, point_num)
+  double precision      , intent(in)  :: r_cusp(nucl_num)
+  double precision      , intent(in)  :: coefficient_t(mo_num,ao_num)
+  double precision      , intent(in)  :: ao_vgl(ao_num,5,point_num)
+  double precision      , intent(out) :: mo_vgl(mo_num,5,point_num)
+  integer*8 :: i,j,k, inucl
+  double precision :: c1, c2, c3, c4, c5
+
+  do j=1,point_num
+     mo_vgl(:,:,j) = 0.d0
+     do k=1,ao_num
+        if (ao_vgl(k,1,j) == 0.d0) cycle
+        inucl = ao_nucl(k)+1
+        if ( (en_distance(inucl,j) < r_cusp(inucl)) .and. (ao_ang_mom(k) == 0) ) cycle
+        c1 = ao_vgl(k,1,j)
+        c2 = ao_vgl(k,2,j)
+        c3 = ao_vgl(k,3,j)
+        c4 = ao_vgl(k,4,j)
+        c5 = ao_vgl(k,5,j)
+        do i=1,mo_num
+           mo_vgl(i,1,j) = mo_vgl(i,1,j) + coefficient_t(i,k) * c1
+           mo_vgl(i,2,j) = mo_vgl(i,2,j) + coefficient_t(i,k) * c2
+           mo_vgl(i,3,j) = mo_vgl(i,3,j) + coefficient_t(i,k) * c3
+           mo_vgl(i,4,j) = mo_vgl(i,4,j) + coefficient_t(i,k) * c4
+           mo_vgl(i,5,j) = mo_vgl(i,5,j) + coefficient_t(i,k) * c5
+        end do
+     end do
+  end do
+  info = QMCKL_SUCCESS
+
+end function qmckl_compute_mo_basis_mo_vgl_cusp_doc_f
+    #+end_src
+
+    #+CALL: generate_c_header(table=qmckl_mo_basis_mo_vgl_cusp_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_vgl_cusp"))
+
+   #+RESULTS:
+   #+begin_src c :tangle (eval h_func) :comments org
+   qmckl_exit_code qmckl_compute_mo_basis_mo_vgl_cusp (
+         const qmckl_context context,
+         const int64_t nucl_num,
+         const int64_t ao_num,
+         const int64_t mo_num,
+         const int64_t point_num,
+         const int64_t* ao_nucl,
+         const int32_t* ao_ang_mom,
+         const double* en_distance,
+         const double* r_cusp,
+         const double* coefficient_t,
+         const double* ao_vgl,
+         double* const mo_vgl ); 
+   #+end_src
+
+   #+CALL: generate_c_header(table=qmckl_mo_basis_mo_vgl_cusp_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_vgl_cusp_doc"))
+
+   #+RESULTS:
+   #+begin_src c :tangle (eval h_func) :comments org
+   qmckl_exit_code qmckl_compute_mo_basis_mo_vgl_cusp_doc (
+         const qmckl_context context,
+         const int64_t nucl_num,
+         const int64_t ao_num,
+         const int64_t mo_num,
+         const int64_t point_num,
+         const int64_t* ao_nucl,
+         const int32_t* ao_ang_mom,
+         const double* en_distance,
+         const double* r_cusp,
+         const double* coefficient_t,
+         const double* ao_vgl,
+         double* const mo_vgl ); 
+   #+end_src
+
+   #+CALL: generate_c_interface(table=qmckl_mo_basis_mo_vgl_cusp_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_vgl_cusp_doc"))
+
+    #+RESULTS:
+    #+begin_src f90 :tangle (eval f) :comments org :exports none
+    integer(c_int32_t) function qmckl_compute_mo_basis_mo_vgl_cusp_doc &
+        (context, &
+         nucl_num, &
+         ao_num, &
+         mo_num, &
+         point_num, &
+         ao_nucl, &
+         ao_ang_mom, &
+         en_distance, &
+         r_cusp, &
+         coefficient_t, &
+         ao_vgl, &
+         mo_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 :: nucl_num
+      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
+      integer (c_int64_t) , intent(in)          :: ao_nucl(ao_num)
+      integer (c_int32_t) , intent(in)          :: ao_ang_mom(ao_num)
+      real    (c_double ) , intent(in)          :: en_distance(nucl_num,point_num)
+      real    (c_double ) , intent(in)          :: r_cusp(nucl_num)
+      real    (c_double ) , intent(in)          :: coefficient_t(ao_num,mo_num)
+      real    (c_double ) , intent(in)          :: ao_vgl(ao_num,5,point_num)
+      real    (c_double ) , intent(out)         :: mo_vgl(mo_num,5,point_num)
+
+      integer(c_int32_t), external :: qmckl_compute_mo_basis_mo_vgl_cusp_doc_f
+      info = qmckl_compute_mo_basis_mo_vgl_cusp_doc_f &
+             (context, &
+         nucl_num, &
+         ao_num, &
+         mo_num, &
+         point_num, &
+         ao_nucl, &
+         ao_ang_mom, &
+         en_distance, &
+         r_cusp, &
+         coefficient_t, &
+         ao_vgl, &
+         mo_vgl)
+
+    end function qmckl_compute_mo_basis_mo_vgl_cusp_doc
+    #+end_src
+
+    #+begin_src c :tangle (eval c) :comments org
+qmckl_exit_code
+qmckl_compute_mo_basis_mo_vgl_cusp (const qmckl_context context,
+                                    const int64_t nucl_num,
+                                    const int64_t ao_num,
+                                    const int64_t mo_num,
+                                    const int64_t point_num,
+                                    const int64_t* ao_nucl,
+                                    const int32_t* ao_ang_mom,
+                                    const double* en_distance,
+                                    const double* r_cusp,
+                                    const double* coefficient_t,
+                                    const double* ao_vgl,
+                                    double* const mo_vgl )
+{
+#ifdef HAVE_HPC
+  return qmckl_compute_mo_basis_mo_vgl_cusp_hpc (context, nucl_num, ao_num, mo_num, point_num,
+                                                 ao_nucl, ao_ang_mom, en_distance, r_cusp,
+                                                 coefficient_t, ao_vgl, mo_vgl );
+#else
+  return qmckl_compute_mo_basis_mo_vgl_cusp_doc (context, nucl_num, ao_num, mo_num, point_num,
+                                                 ao_nucl, ao_ang_mom, en_distance, r_cusp,
+                                                 coefficient_t, ao_vgl, mo_vgl );
+#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_vgl_cusp_hpc (const qmckl_context context,
+                                        const int64_t nucl_num,
+                                        const int64_t ao_num,
+                                        const int64_t mo_num,
+                                        const int64_t point_num,
+                                        const int64_t* ao_nucl,
+                                        const int32_t* ao_ang_mom,
+                                        const double* en_distance,
+                                        const double* r_cusp,
+                                        const double* coefficient_t,
+                                        const double* ao_vgl,
+                                        double* const mo_vgl ); 
+#endif
+    #+end_src
+
+    #+begin_src c :tangle (eval c) :comments org
+#ifdef HAVE_HPC
+qmckl_exit_code
+qmckl_compute_mo_basis_mo_vgl_cusp_hpc (const qmckl_context context,
+                                        const int64_t nucl_num,
+                                        const int64_t ao_num,
+                                        const int64_t mo_num,
+                                        const int64_t point_num,
+                                        const int64_t* ao_nucl,
+                                        const int32_t* ao_ang_mom,
+                                        const double* en_distance,
+                                        const double* r_cusp,
+                                        const double* coefficient_t,
+                                        const double* ao_vgl,
+                                        double* const mo_vgl )
+{
+  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_vgl[ipoint*5*mo_num]);
+    double* restrict const vgl2 =  vgl1 + mo_num;
+    double* restrict const vgl3 =  vgl1 + (mo_num << 1);
+    double* restrict const vgl4 =  vgl1 + (mo_num << 1) + mo_num;
+    double* restrict const vgl5 =  vgl1 + (mo_num << 2);
+
+    const double* restrict avgl1 = &(ao_vgl[ipoint*5*ao_num]);
+    const double* restrict avgl2 = avgl1 + ao_num;
+    const double* restrict avgl3 = avgl1 + (ao_num << 1);
+    const double* restrict avgl4 = avgl1 + (ao_num << 1) + ao_num;
+    const double* restrict avgl5 = avgl1 + (ao_num << 2);
+
+    for (int64_t i=0 ; i<mo_num ; ++i) {
+      vgl1[i] = 0.;
+      vgl2[i] = 0.;
+      vgl3[i] = 0.;
+      vgl4[i] = 0.;
+      vgl5[i] = 0.;
+    }
+
+    const double* restrict ria   = &(en_distance[ipoint*nucl_num]);
+
+    int64_t nidx=0;
+    int64_t idx[ao_num];
+    double  av1[ao_num];
+    double  av2[ao_num];
+    double  av3[ao_num];
+    double  av4[ao_num];
+    double  av5[ao_num];
+    for (int64_t k=0 ; k<ao_num ; ++k) {
+      if (avgl1[k] != 0.) {
+        const int64_t inucl = ao_nucl[k];
+        if (ria[inucl] > r_cusp[inucl] || ao_ang_mom[k] > 0) {
+          idx[nidx] = k;
+          av1[nidx] = avgl1[k];
+          av2[nidx] = avgl2[k];
+          av3[nidx] = avgl3[k];
+          av4[nidx] = avgl4[k];
+          av5[nidx] = avgl5[k];
+          ++nidx;
+        }
+      }
+    }
+
+    int64_t n=0;
+
+    for (n=0 ; n < nidx-4 ; n+=4) {
+      const double* restrict ck1 = coefficient_t + idx[n  ]*mo_num;
+      const double* restrict ck2 = coefficient_t + idx[n+1]*mo_num;
+      const double* restrict ck3 = coefficient_t + idx[n+2]*mo_num;
+      const double* restrict ck4 = coefficient_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];
+
+      const double a12 = av2[n  ];
+      const double a22 = av2[n+1];
+      const double a32 = av2[n+2];
+      const double a42 = av2[n+3];
+
+      const double a13 = av3[n  ];
+      const double a23 = av3[n+1];
+      const double a33 = av3[n+2];
+      const double a43 = av3[n+3];
+
+      const double a14 = av4[n  ];
+      const double a24 = av4[n+1];
+      const double a34 = av4[n+2];
+      const double a44 = av4[n+3];
+
+      const double a15 = av5[n  ];
+      const double a25 = av5[n+1];
+      const double a35 = av5[n+2];
+      const double a45 = av5[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;
+        vgl2[i] = vgl2[i] + ck1[i] * a12 + ck2[i] * a22 + ck3[i] * a32 + ck4[i] * a42;
+        vgl3[i] = vgl3[i] + ck1[i] * a13 + ck2[i] * a23 + ck3[i] * a33 + ck4[i] * a43;
+        vgl4[i] = vgl4[i] + ck1[i] * a14 + ck2[i] * a24 + ck3[i] * a34 + ck4[i] * a44;
+        vgl5[i] = vgl5[i] + ck1[i] * a15 + ck2[i] * a25 + ck3[i] * a35 + ck4[i] * a45;
+      }
+    }
+
+    for (int64_t m=n ; m < nidx ; m+=1) {
+      const double* restrict ck = coefficient_t + idx[m]*mo_num;
+      const double a1 = av1[m];
+      const double a2 = av2[m];
+      const double a3 = av3[m];
+      const double a4 = av4[m];
+      const double a5 = av5[m];
+
+#ifdef HAVE_OPENMP
+  #pragma omp simd
+#endif
+      for (int64_t i=0 ; i<mo_num ; ++i) {
+        vgl1[i] += ck[i] * a1;
+        vgl2[i] += ck[i] * a2;
+        vgl3[i] += ck[i] * a3;
+        vgl4[i] += ck[i] * a4;
+        vgl5[i] += ck[i] * a5;
+      }
+    }
+  }
+  return QMCKL_SUCCESS;
+}
+#endif
+    #+end_src
+
 ** Rescaling of MO coefficients
 
    When evaluating Slater determinants, the value of the determinants
@@ -1799,7 +2583,7 @@ assert (rc == QMCKL_SUCCESS);
 
 const double  * mo_coefficient          =  &(chbrclf_mo_coef[0]);
 
-rc = qmckl_set_mo_basis_coefficient(context, mo_coefficient);
+rc = qmckl_set_mo_basis_coefficient(context, mo_coefficient, chbrclf_mo_num*chbrclf_ao_num);
 assert (rc == QMCKL_SUCCESS);
 
 assert(qmckl_mo_basis_provided(context));
diff --git a/org/qmckl_trexio.org b/org/qmckl_trexio.org
index 7761799..0deb13e 100644
--- a/org/qmckl_trexio.org
+++ b/org/qmckl_trexio.org
@@ -1030,7 +1030,7 @@ qmckl_trexio_read_mo_X(qmckl_context context, trexio_t* const file)
                              trexio_string_of_error(rcio));
     }
 
-    rc = qmckl_set_mo_basis_coefficient(context, mo_coef);
+    rc = qmckl_set_mo_basis_coefficient(context, mo_coef, ao_num*mo_num);
 
     qmckl_free(context, mo_coef);
     mo_coef = NULL;