fix bug in ee single

org/qmckl_forces.org

@@ -31,6 +31,7 @@
 #include <stdbool.h>
 #include <stdio.h>
 #include "n2.h"
+#include "chbrclf.h"
 #include "qmckl_jastrow_champ_private_func.h"
 #include "qmckl_jastrow_champ_single_private_func.h"
 #include "qmckl_forces_private_func.h"
@@ -97,6 +98,10 @@ typedef struct qmckl_forces_struct{
     uint64_t  forces_jastrow_een_g_date;
     double  * restrict forces_jastrow_een_l;
     uint64_t  forces_jastrow_een_l_date;
+    double  * restrict forces_ao_value;
+    uint64_t  forces_ao_value_date;
+    double  * restrict forces_mo_value;
+    uint64_t  forces_mo_value_date;
 } qmckl_forces_struct;
     #+end_src
 
@@ -3700,6 +3705,909 @@ printf("OK\n");
 
       #+end_src
 
+* Reset test for orbitals
+
+   #+begin_src c :tangle (eval c_test)
+
+  rc = qmckl_context_destroy(context);
+  assert (rc == QMCKL_SUCCESS);
+   context = qmckl_context_create();
+
+   nucl_num      = chbrclf_nucl_num;
+   nucl_charge   = chbrclf_charge;
+   nucl_coord    = &(chbrclf_nucl_coord[0][0]);
+
+  rc = qmckl_set_nucleus_num (context, nucl_num);
+  assert(rc == QMCKL_SUCCESS);
+
+  rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), 3*nucl_num);
+  assert(rc == QMCKL_SUCCESS);
+
+  rc = qmckl_set_nucleus_charge(context, nucl_charge, nucl_num);
+  assert(rc == QMCKL_SUCCESS);
+
+  assert(qmckl_nucleus_provided(context));
+
+
+  const int64_t    shell_num         =  chbrclf_shell_num;
+  const int64_t    prim_num          =  chbrclf_prim_num;
+  const int64_t    ao_num            =  chbrclf_ao_num;
+  const int64_t *  nucleus_index     =  &(chbrclf_basis_nucleus_index[0]);
+  const int64_t *  nucleus_shell_num =  &(chbrclf_basis_nucleus_shell_num[0]);
+  const int32_t *  shell_ang_mom     =  &(chbrclf_basis_shell_ang_mom[0]);
+  const int64_t *  shell_prim_num    =  &(chbrclf_basis_shell_prim_num[0]);
+  const int64_t *  shell_prim_index  =  &(chbrclf_basis_shell_prim_index[0]);
+  const double  *  shell_factor      =  &(chbrclf_basis_shell_factor[0]);
+  const double  *  exponent          =  &(chbrclf_basis_exponent[0]);
+  const double  *  coefficient       =  &(chbrclf_basis_coefficient[0]);
+  const double  *  prim_factor       =  &(chbrclf_basis_prim_factor[0]);
+  const double  *  ao_factor         =  &(chbrclf_basis_ao_factor[0]);
+
+  const char typ = 'G';
+
+  assert(!qmckl_ao_basis_provided(context));
+
+  rc = qmckl_set_ao_basis_type (context, typ);
+  assert(rc == QMCKL_SUCCESS);
+  assert(!qmckl_ao_basis_provided(context));
+
+  rc = qmckl_set_ao_basis_shell_num (context, shell_num);
+  assert(rc == QMCKL_SUCCESS);
+  assert(!qmckl_ao_basis_provided(context));
+
+  rc = qmckl_set_ao_basis_prim_num (context, prim_num);
+  assert(rc == QMCKL_SUCCESS);
+  assert(!qmckl_ao_basis_provided(context));
+
+  rc = qmckl_set_ao_basis_nucleus_index (context, nucleus_index, nucl_num);
+  assert(rc == QMCKL_SUCCESS);
+  assert(!qmckl_ao_basis_provided(context));
+
+  rc = qmckl_set_ao_basis_nucleus_index (context, nucleus_index, nucl_num);
+  assert(rc == QMCKL_ALREADY_SET);
+
+  rc = qmckl_set_ao_basis_nucleus_shell_num (context, nucleus_shell_num, nucl_num);
+  assert(rc == QMCKL_SUCCESS);
+  assert(!qmckl_ao_basis_provided(context));
+
+  rc = qmckl_set_ao_basis_shell_ang_mom (context, shell_ang_mom, shell_num);
+  assert(rc == QMCKL_SUCCESS);
+  assert(!qmckl_ao_basis_provided(context));
+
+  rc = qmckl_set_ao_basis_shell_factor  (context, shell_factor, shell_num);
+  assert(rc == QMCKL_SUCCESS);
+  assert(!qmckl_ao_basis_provided(context));
+
+  rc = qmckl_set_ao_basis_shell_prim_num (context, shell_prim_num, shell_num);
+  assert(rc == QMCKL_SUCCESS);
+  assert(!qmckl_ao_basis_provided(context));
+
+  rc = qmckl_set_ao_basis_shell_prim_index (context, shell_prim_index, shell_num);
+  assert(rc == QMCKL_SUCCESS);
+  assert(!qmckl_ao_basis_provided(context));
+
+  rc = qmckl_set_ao_basis_exponent      (context, exponent, prim_num);
+  assert(rc == QMCKL_SUCCESS);
+  assert(!qmckl_ao_basis_provided(context));
+
+  rc = qmckl_set_ao_basis_coefficient   (context, coefficient, prim_num);
+  assert(rc == QMCKL_SUCCESS);
+  assert(!qmckl_ao_basis_provided(context));
+
+  rc = qmckl_set_ao_basis_prim_factor (context, prim_factor, prim_num);
+  assert(rc == QMCKL_SUCCESS);
+
+  rc = qmckl_set_ao_basis_ao_num(context, ao_num);
+  assert(rc == QMCKL_SUCCESS);
+
+  rc = qmckl_set_ao_basis_ao_factor (context, ao_factor, ao_num);
+  assert(rc == QMCKL_SUCCESS);
+
+  assert(qmckl_ao_basis_provided(context));
+
+  int64_t    shell_num_test        ;
+  int64_t    prim_num_test         ;
+  int64_t    ao_num_test           ;
+  int64_t *  nucleus_index_test    ;
+  int64_t *  nucleus_shell_num_test;
+  int32_t *  shell_ang_mom_test    ;
+  int64_t *  shell_prim_num_test   ;
+  int64_t *  shell_prim_index_test ;
+  double  *  shell_factor_test     ;
+  double  *  exponent_test         ;
+  double  *  coefficient_test      ;
+  double  *  prim_factor_test      ;
+  double  *  ao_factor_test        ;
+  char    typ_test ;
+
+
+  rc = qmckl_get_ao_basis_type (context, &typ_test);
+  assert (rc == QMCKL_SUCCESS);
+  assert(typ == typ_test);
+
+  rc = qmckl_get_ao_basis_shell_num (context, &shell_num_test);
+  assert (rc == QMCKL_SUCCESS);
+  assert(shell_num == shell_num_test);
+
+  rc = qmckl_get_ao_basis_prim_num (context, &prim_num_test);
+  assert (rc == QMCKL_SUCCESS);
+  assert(prim_num == prim_num_test);
+
+  nucleus_index_test = (int64_t*) malloc (nucl_num * sizeof(int64_t));
+  rc = qmckl_get_ao_basis_nucleus_index (context, nucleus_index_test, nucl_num);
+  assert (rc == QMCKL_SUCCESS);
+  for (int64_t i=0 ; i < nucl_num ; ++i) {
+    assert(nucleus_index_test[i] == nucleus_index[i]);
+  }
+  free(nucleus_index_test);
+
+  nucleus_shell_num_test = (int64_t*) malloc ( nucl_num * sizeof(int64_t));
+  rc = qmckl_get_ao_basis_nucleus_shell_num (context, nucleus_shell_num_test, nucl_num);
+  assert (rc == QMCKL_SUCCESS);
+
+  for (int64_t i=0 ; i < nucl_num ; ++i) {
+    assert(nucleus_shell_num_test[i] == nucleus_shell_num[i]);
+  }
+  free(nucleus_shell_num_test);
+
+  shell_ang_mom_test = (int32_t*) malloc ( shell_num * sizeof(int32_t));
+  rc = qmckl_get_ao_basis_shell_ang_mom (context, shell_ang_mom_test, shell_num);
+  assert (rc == QMCKL_SUCCESS);
+
+  for (int64_t i=0 ; i < shell_num ; ++i) {
+    assert(shell_ang_mom_test[i] == shell_ang_mom[i]);
+  }
+  free(shell_ang_mom_test);
+
+  shell_factor_test = (double*) malloc ( shell_num * sizeof(double));
+  rc = qmckl_get_ao_basis_shell_factor  (context, shell_factor_test, shell_num);
+  assert (rc == QMCKL_SUCCESS);
+
+  for (int64_t i=0 ; i < shell_num ; ++i) {
+    assert(shell_factor_test[i] == shell_factor[i]);
+  }
+  free(shell_factor_test);
+
+  shell_prim_num_test = (int64_t*) malloc ( shell_num * sizeof(int64_t));
+  rc = qmckl_get_ao_basis_shell_prim_num (context, shell_prim_num_test, shell_num);
+  assert (rc == QMCKL_SUCCESS);
+
+  for (int64_t i=0 ; i < shell_num ; ++i) {
+    assert(shell_prim_num_test[i] == shell_prim_num[i]);
+  }
+  free(shell_prim_num_test);
+
+  shell_prim_index_test = (int64_t*) malloc ( shell_num * sizeof(int64_t));
+  rc = qmckl_get_ao_basis_shell_prim_index (context, shell_prim_index_test, shell_num);
+  assert (rc == QMCKL_SUCCESS);
+
+  for (int64_t i=0 ; i < shell_num ; ++i) {
+    assert(shell_prim_index_test[i] == shell_prim_index[i]);
+  }
+  free(shell_prim_index_test);
+
+  exponent_test = (double*) malloc ( prim_num * sizeof(double));
+  rc = qmckl_get_ao_basis_exponent(context, exponent_test, prim_num);
+  assert (rc == QMCKL_SUCCESS);
+
+  for (int64_t i=0 ; i < prim_num ; ++i) {
+    assert(exponent_test[i] == exponent[i]);
+  }
+  free(exponent_test);
+
+  coefficient_test = (double*) malloc ( prim_num * sizeof(double));
+  rc = qmckl_get_ao_basis_coefficient(context, coefficient_test, prim_num);
+  assert (rc == QMCKL_SUCCESS);
+
+  for (int64_t i=0 ; i < prim_num ; ++i) {
+    assert(coefficient_test[i] == coefficient[i]);
+  }
+  free(coefficient_test);
+
+  prim_factor_test = (double*) malloc ( prim_num * sizeof(double));
+  rc = qmckl_get_ao_basis_prim_factor (context, prim_factor_test, prim_num);
+  assert (rc == QMCKL_SUCCESS);
+
+  for (int64_t i=0 ; i < prim_num ; ++i) {
+    assert(prim_factor_test[i] == prim_factor[i]);
+  }
+  free(prim_factor_test);
+
+  rc = qmckl_get_ao_basis_ao_num(context, &ao_num_test);
+  assert(ao_num == ao_num_test);
+
+  ao_factor_test = (double*) malloc ( ao_num * sizeof(double));
+  rc = qmckl_get_ao_basis_ao_factor (context, ao_factor_test, ao_num);
+  assert (rc == QMCKL_SUCCESS);
+
+  for (int64_t i=0 ; i < ao_num ; ++i) {
+    assert(ao_factor_test[i] == ao_factor[i]);
+  }
+  free(ao_factor_test);
+
+#define walk_num 1 // chbrclf_walk_num
+#define elec_num chbrclf_elec_num
+#define prim_num chbrclf_prim_num
+
+  elec_up_num   = chbrclf_elec_up_num;
+  elec_dn_num   = chbrclf_elec_dn_num;
+  elec_coord    = &(chbrclf_elec_coord[0][0][0]);
+
+  rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num);
+  assert (rc == QMCKL_SUCCESS);
+
+  assert(qmckl_electron_provided(context));
+
+  int64_t point_num = elec_num;
+
+  rc = qmckl_set_point(context, 'N', point_num, elec_coord, point_num*3);
+  assert(rc == QMCKL_SUCCESS);
+
+  int64_t mo_num = chbrclf_mo_num;
+  rc = qmckl_set_mo_basis_mo_num(context, mo_num);
+  assert (rc == QMCKL_SUCCESS);
+
+  const double  * mo_coefficient          =  &(chbrclf_mo_coef[0]);
+
+  rc = qmckl_set_mo_basis_coefficient(context, mo_coefficient, chbrclf_mo_num*chbrclf_ao_num);
+  assert (rc == QMCKL_SUCCESS);
+
+    #+end_src
+
+
+* Force of AO value
+
+Here we compute the forces of the AO value. 
+
+
+** Get
+     #+begin_src c :comments org :tangle (eval h_func) :noweb yes
+qmckl_exit_code
+qmckl_get_forces_ao_value(qmckl_context context,
+                                     double* const forces_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_forces_ao_value(qmckl_context context,
+                                     double* const forces_ao_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_forces_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->nucleus.num * 3 * ctx->point.num;
+  if (size_max < sze) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_3,
+                           "qmckl_get_forces_ao_value",
+                           "Array too small. Expected walk_num*nucl_num*point_num*3");
+  }
+  memcpy(forces_ao_value, ctx->forces.forces_ao_value, sze*sizeof(double));
+
+  return QMCKL_SUCCESS;
+}
+     #+end_src
+
+
+** Provide                                                       :noexport:
+
+      #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
+qmckl_exit_code qmckl_provide_forces_ao_value(qmckl_context context);
+      #+end_src
+  
+      #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
+qmckl_exit_code qmckl_provide_forces_ao_value(qmckl_context context)
+{
+  qmckl_exit_code rc = QMCKL_SUCCESS;
+
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_CONTEXT,
+                           "qmckl_provide_forces_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_provide_forces_ao_value",
+                           NULL);
+  }
+
+  /* Compute if necessary */
+  if (ctx->point.date > ctx->forces.forces_ao_value_date) {
+
+    qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
+    mem_info.size = ctx->ao_basis.ao_num * 3 * ctx->nucleus.num * ctx->point.num * sizeof(double);
+
+    if (ctx->forces.forces_ao_value != NULL) {
+      qmckl_memory_info_struct mem_info_test = qmckl_memory_info_struct_zero;
+      rc = qmckl_get_malloc_info(context, ctx->forces.forces_ao_value, &mem_info_test);
+
+      /* if rc != QMCKL_SUCCESS, we are maybe in an _inplace function because the
+         memory was not allocated with qmckl_malloc */
+
+      if ((rc == QMCKL_SUCCESS) && (mem_info_test.size != mem_info.size)) {
+        rc = qmckl_free(context, ctx->forces.forces_ao_value);
+        assert (rc == QMCKL_SUCCESS);
+        ctx->forces.forces_ao_value = NULL;
+      }
+    }
+
+    /* Allocate array */
+    if (ctx->forces.forces_ao_value == NULL) {
+
+      double* forces_ao_value = (double*) qmckl_malloc(context, mem_info);
+
+      if (forces_ao_value == NULL) {
+        return qmckl_failwith( context,
+                               QMCKL_ALLOCATION_FAILED,
+                               "qmckl_forces_ao_value",
+                               NULL);
+      }
+      ctx->forces.forces_ao_value = forces_ao_value;
+    }
+    rc = qmckl_provide_ao_basis_ao_vgl(context);
+    if (rc != QMCKL_SUCCESS) {
+      return qmckl_failwith( context, rc, "qmckl_provide_ao_basis_ao_vgl", NULL);
+    }
+
+    rc = qmckl_compute_forces_ao_value_doc(context,
+                                  ctx->ao_basis.ao_num,
+                                  ctx->ao_basis.shell_num,
+                                  ctx->point.num,
+                                  ctx->nucleus.num,
+                                  ctx->ao_basis.nucleus_index,
+                                  ctx->ao_basis.nucleus_shell_num,
+                                  ctx->ao_basis.shell_ang_mom,
+                                  ctx->ao_basis.ao_factor,
+                                  ctx->ao_basis.ao_vgl,
+                                  ctx->forces.forces_ao_value);
+
+    if (rc != QMCKL_SUCCESS) {
+      return rc;
+    }
+
+    ctx->forces.forces_ao_value_date = ctx->date;
+  }
+
+  return QMCKL_SUCCESS;
+}
+     #+end_src
+
+** Compute
+   :PROPERTIES:
+   :Name:     qmckl_compute_forces_ao_value
+   :CRetType: qmckl_exit_code
+   :FRetType: qmckl_exit_code
+   :END:
+
+     #+NAME: qmckl_forces_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                             |
+    | ~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                 |
+    | ~shell_ang_mom~       | ~int32_t[shell_num]~              | in     | Angular momentum of each shell               |
+    | ~ao_factor~           | ~double[ao_num]~                  | in     | Normalization factor of the AOs              |
+    | ~ao_vgl~           | ~double[point_num][5][shell_num]~ | in     | Value, gradients and Laplacian of the shells |
+    | ~forces_ao_value~     | ~double[nucl_num][3][point_num][ao_num]~     | out    | Forces of the AOs              |
+
+     #+begin_src f90 :comments org :tangle (eval f) :noweb yes
+function qmckl_compute_forces_ao_value_doc(context, &
+     ao_num, shell_num, point_num, nucl_num, &
+     nucleus_index, nucleus_shell_num, &
+     shell_ang_mom, ao_factor, ao_vgl, forces_ao_value) &
+     bind(C) result(info)
+  use qmckl_constants
+  use qmckl, only : qmckl_ao_polynomial_vgl, qmckl_get_numprec_precision
+  implicit none
+  integer (qmckl_context), 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
+  integer (c_int64_t) , intent(in)          :: nucleus_index(nucl_num)
+  integer (c_int64_t) , intent(in)          :: nucleus_shell_num(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)          :: ao_vgl(ao_num,5,point_num)
+  real    (c_double ) , intent(out)         :: forces_ao_value(ao_num,point_num,3,nucl_num)
+  integer(qmckl_exit_code) :: info
+
+  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         , allocatable  ::  ao_index(:)
+  allocate(ao_index(ao_num))
+
+  forces_ao_value = 0.d0
+
+  ! 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
+
+
+  do ipoint = 1, point_num
+     do inucl=1,nucl_num
+        ! 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
+
+              forces_ao_value(k,ipoint,1,inucl) = -ao_vgl(k,2,ipoint)
+              forces_ao_value(k,ipoint,2,inucl) = -ao_vgl(k,3,ipoint)
+              forces_ao_value(k,ipoint,3,inucl) = -ao_vgl(k,4,ipoint)
+
+              k = k+1
+           end do
+        end do
+     end do
+  end do
+
+end function qmckl_compute_forces_ao_value_doc
+     #+end_src
+
+  #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
+    qmckl_exit_code qmckl_compute_forces_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 int64_t* nucleus_index,
+          const int64_t* nucleus_shell_num,
+          const int32_t* shell_ang_mom,
+          const double* ao_factor,
+          const double* ao_vgl,
+          double* const forces_ao_value );
+    #+end_src
+
+
+** Test
+
+   #+begin_src c :tangle (eval c_test)
+printf("Forces AO value\n");
+
+rc = qmckl_set_nucleus_coord(context, 'T', &(nucl_coord[0]), 3*nucl_num);
+assert(rc == QMCKL_SUCCESS);
+
+double forces_ao_value[nucl_num][3][point_num][ao_num];
+rc = qmckl_get_forces_ao_value(context, &forces_ao_value[0][0][0][0], 3*nucl_num*ao_num*point_num);
+assert(rc == QMCKL_SUCCESS);
+
+double finite_difference_force_ao_value[3][nucl_num][point_num][ao_num];
+
+nucleus_coord = (double*) malloc(3 * nucl_num * sizeof(double));
+if (nucleus_coord == NULL) {
+  return QMCKL_ALLOCATION_FAILED;
+}
+
+rc = qmckl_get_nucleus_coord(context, 'N', nucleus_coord, 3*nucl_num);
+
+temp_coord = (double*) malloc(3 * nucl_num * sizeof(double));
+if (temp_coord == NULL) {
+  free(nucleus_coord);
+  return QMCKL_ALLOCATION_FAILED;
+}
+double ao_output[point_num][ao_num];
+
+// Copy original coordinates
+for (int i = 0; i < 3 * nucl_num; i++) {
+  temp_coord[i] = nucleus_coord[i];
+}
+
+for (int64_t a = 0; a < nucl_num; a++) {
+  for (int64_t k = 0; k < 3; k++) {
+    for (int64_t m = -4; m <= 4; m++) {
+      
+      // Apply finite difference displacement
+      temp_coord[k+a*3] = nucleus_coord[k+3*a] + (double) m * delta_x;
+
+      // Update coordinates in the context
+      rc = qmckl_set_nucleus_coord(context, 'N', temp_coord, 3*nucl_num);
+      assert(rc == QMCKL_SUCCESS);
+      
+      rc = qmckl_context_touch(context);
+      assert(rc == QMCKL_SUCCESS);
+      
+      // Call the provided function
+      rc = qmckl_get_ao_basis_ao_value(context,&ao_output[0][0], point_num*ao_num);
+      assert(rc == QMCKL_SUCCESS);
+
+      // Accumulate derivative using finite-difference coefficients
+      for (int i = 0; i < point_num; i++) {
+        for (int j = 0; j < ao_num; j++) {
+          if (m == -4) {
+            finite_difference_force_ao_value[k][a][i][j] = 0.0;
+          }
+          finite_difference_force_ao_value[k][a][i][j] += coef[m + 4] * ao_output[i][j]/delta_x;
+        }
+      }
+    }
+    temp_coord[k+a*3] = nucleus_coord[k+3*a];
+  }
+}
+
+// Reset coordinates in the context
+rc = qmckl_set_nucleus_coord(context, 'N', temp_coord, 3*nucl_num);
+assert(rc == QMCKL_SUCCESS);
+
+rc = qmckl_context_touch(context);
+assert(rc == QMCKL_SUCCESS);
+    
+free(nucleus_coord);
+free(temp_coord);
+
+
+for (int j = 0; j < ao_num; j++){
+  for (int i = 0; i < point_num; i++){
+      for (int a = 0; a < nucl_num; a++) {
+        for (int k = 0; k < 3; k++){
+            //printf("k=%i a=%i i=%i j=%i\n", k, a, i, j);
+            //printf("%.10f\t", finite_difference_force_ao_value[k][a][i][j]);
+            //printf("%.10f\n", forces_ao_value[a][k][i][j]);
+            assert(fabs(finite_difference_force_ao_value[k][a][i][j] -  forces_ao_value[a][k][i][j]) < 1.e-10);
+
+        }
+    }
+  }
+}
+
+printf("OK\n");
+
+      #+end_src
+
+* Force of MO value
+
+** Get
+
+    #+begin_src c :comments org :tangle (eval h_func) :noweb yes
+qmckl_exit_code
+qmckl_get_forces_mo_value(qmckl_context context,
+                          double* const forces_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_forces_mo_value(qmckl_context context,
+                          double* const forces_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_forces_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 * 3 * ctx->mo_basis.mo_num * ctx->nucleus.num;
+  if (size_max < sze) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_3,
+                           "qmckl_get_forces_mo_value",
+                           "input array too small");
+  }
+  memcpy(forces_mo_value, ctx->forces.forces_mo_value, sze * sizeof(double));
+
+  return QMCKL_SUCCESS;
+}
+    #+end_src
+
+    #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
+  interface
+     integer(qmckl_exit_code) function qmckl_get_forces_mo_value (context, &
+          forces_mo_value, size_max) bind(C)
+       use, intrinsic :: iso_c_binding
+       import
+       implicit none
+
+       integer (c_int64_t) , intent(in)  , value :: context
+       real(c_double),       intent(out)         :: forces_mo_value(*)
+       integer (c_int64_t) , intent(in)  , value :: size_max
+     end function qmckl_get_forces_mo_value
+  end interface
+    #+end_src
+
+** Provide
+
+#+begin_src c :comments org :tangle (eval c) :noweb yes :export none
+qmckl_exit_code qmckl_provide_forces_mo_value(qmckl_context context)
+{
+
+  qmckl_exit_code rc = QMCKL_SUCCESS;
+
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_CONTEXT,
+                           "qmckl_provide_forces_mo_value",
+                           NULL);
+  }
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
+  assert (ctx != NULL);
+
+  if (!ctx->mo_basis.provided) {
+    return qmckl_failwith( context,
+                           QMCKL_NOT_PROVIDED,
+                           "qmckl_provide_mo_basis_mo_vgl",
+                           NULL);
+  }
+
+  /* Compute if necessary */
+  if (ctx->point.date > ctx->forces.forces_mo_value_date) {
+
+    qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
+    mem_info.size = 3 * ctx->mo_basis.mo_num * ctx->point.num * ctx->nucleus.num * sizeof(double);
+
+    if (ctx->forces.forces_mo_value != NULL) {
+      qmckl_memory_info_struct mem_info_test = qmckl_memory_info_struct_zero;
+      rc = qmckl_get_malloc_info(context, ctx->forces.forces_mo_value, &mem_info_test);
+
+      /* if rc != QMCKL_SUCCESS, we are maybe in an _inplace function because the
+         memory was not allocated with qmckl_malloc */
+
+      if ((rc == QMCKL_SUCCESS) && (mem_info_test.size != mem_info.size)) {
+        rc = qmckl_free(context, ctx->forces.forces_mo_value);
+        assert (rc == QMCKL_SUCCESS);
+        ctx->forces.forces_mo_value = NULL;
+      }
+    }
+
+    /* Allocate array */
+    if (ctx->forces.forces_mo_value == NULL) {
+
+      double* forces_mo_value = (double*) qmckl_malloc(context, mem_info);
+
+      if (forces_mo_value == NULL) {
+        return qmckl_failwith( context,
+                               QMCKL_ALLOCATION_FAILED,
+                               "qmckl_forces_mo_value",
+                               NULL);
+      }
+      ctx->forces.forces_mo_value = forces_mo_value;
+    }
+
+    rc = qmckl_provide_forces_ao_value(context);
+    if (rc != QMCKL_SUCCESS) {
+      return qmckl_failwith( context,
+                             QMCKL_NOT_PROVIDED,
+                             "qmckl_forces_ao_value",
+                             NULL);
+    }
+
+    rc = qmckl_compute_forces_mo_value_doc(context,
+                                    ctx->nucleus.num,
+                                    ctx->ao_basis.ao_num,
+                                    ctx->mo_basis.mo_num,
+                                    ctx->point.num,
+                                    ctx->mo_basis.coefficient_t,
+                                    ctx->forces.forces_ao_value,
+                                    ctx->forces.forces_mo_value);
+  }
+}
+#+end_src
+
+** Compute
+   :PROPERTIES:
+   :Name:     qmckl_compute_forces_mo_value
+   :CRetType: qmckl_exit_code
+   :FRetType: qmckl_exit_code
+   :END:
+
+    #+NAME: qmckl_forces_mo_value_args
+    | Variable            | Type                           | In/Out | Description                                     |
+    |---------------------+--------------------------------+--------+-------------------------------------------------|
+    | ~context~           | ~qmckl_context~                | in     | Global state                                    |
+    | ~nucl_num~            | ~int64_t~                      | in     | Number of AOs                                   |
+    | ~ao_num~            | ~int64_t~                      | in     | Number of AOs                                   |
+    | ~mo_num~            | ~int64_t~                      | in     | Number of MOs                                   |
+    | ~point_num~         | ~int64_t~                      | in     | Number of points                                |
+    | ~coefficient_t~     | ~double[mo_num][ao_num]~       | in     | Transpose of the AO to MO transformation matrix |
+    | ~forces_ao_value~            | ~double[nucl_num][3][point_num][ao_num]~ | in     | Value, gradients and Laplacian of the AOs       |
+    | ~forces_mo_value~            | ~double[nucl_num][3][point_num][mo_num]~ | out    | Value, gradients and Laplacian 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(qmckl_exit_code)  function qmckl_compute_forces_mo_value_doc(context, &
+     nucl_num,ao_num, mo_num, point_num, &
+     coefficient_t, forces_ao_value, forces_mo_value) &
+     result(info) bind(C)
+  use, intrinsic :: iso_c_binding
+  use qmckl
+  implicit none
+  integer(qmckl_context), intent(in)  :: context
+  integer(c_int64_t)    , intent(in), value  ::nucl_num, ao_num, mo_num, point_num
+  real(c_double)      , intent(in)  :: forces_ao_value(ao_num,point_num,3,nucl_num)
+  real(c_double)      , intent(in)  :: coefficient_t(mo_num,ao_num)
+  real(c_double)      , intent(out) :: forces_mo_value(mo_num,point_num,3,nucl_num)
+  integer*8 :: i,j,k,a
+  double precision :: c1, c2, c3
+
+
+  info = QMCKL_SUCCESS
+ forces_mo_value = 0.0d0
+  do j=1,point_num
+     do k=1,ao_num
+      do a=1,nucl_num
+          c1 = forces_ao_value(k,j,1,a)
+          c2 = forces_ao_value(k,j,2,a)
+          c3 = forces_ao_value(k,j,3,a)
+          do i=1,mo_num
+            forces_mo_value(i,j,1,a) = forces_mo_value(i,j,1,a) + coefficient_t(i,k) * c1
+            forces_mo_value(i,j,2,a) = forces_mo_value(i,j,2,a) + coefficient_t(i,k) * c2
+            forces_mo_value(i,j,3,a) = forces_mo_value(i,j,3,a) + coefficient_t(i,k) * c3
+          end do
+      end do
+     end do
+  end do
+
+
+end function qmckl_compute_forces_mo_value_doc
+    #+end_src
+
+
+   #+begin_src c :tangle (eval h_private_func) :comments org
+   qmckl_exit_code qmckl_compute_forces_mo_value_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 double* coefficient_t,
+         const double* forces_ao_value,
+         double* const forces_mo_value );
+   #+end_src
+
+** Test
+
+   #+begin_src c :tangle (eval c_test)
+printf("Forces MO value\n");
+
+rc = qmckl_set_nucleus_coord(context, 'T', &(nucl_coord[0]), 3*nucl_num);
+assert(rc == QMCKL_SUCCESS);
+
+double forces_mo_value[nucl_num][3][point_num][mo_num];
+rc = qmckl_get_forces_mo_value(context, &forces_mo_value[0][0][0][0], 3*nucl_num*mo_num*point_num);
+assert(rc == QMCKL_SUCCESS);
+
+double finite_difference_force_mo_value[3][nucl_num][point_num][mo_num];
+printf("Mo num %i %i\n", mo_num, ao_num);
+
+nucleus_coord = (double*) malloc(3 * nucl_num * sizeof(double));
+if (nucleus_coord == NULL) {
+  return QMCKL_ALLOCATION_FAILED;
+}
+
+rc = qmckl_get_nucleus_coord(context, 'N', nucleus_coord, 3*nucl_num);
+
+temp_coord = (double*) malloc(3 * nucl_num * sizeof(double));
+if (temp_coord == NULL) {
+  free(nucleus_coord);
+  return QMCKL_ALLOCATION_FAILED;
+}
+double mo_output[point_num][mo_num];
+
+// Copy original coordinates
+for (int i = 0; i < 3 * nucl_num; i++) {
+  temp_coord[i] = nucleus_coord[i];
+}
+
+
+
+for (int64_t a = 0; a < nucl_num; a++) {
+  for (int64_t k = 0; k < 3; k++) {
+    for (int64_t m = -4; m <= 4; m++) {
+      
+      // Apply finite difference displacement
+      temp_coord[k+a*3] = nucleus_coord[k+3*a] + (double) m * delta_x;
+
+      // Update coordinates in the context
+      rc = qmckl_set_nucleus_coord(context, 'N', temp_coord, 3*nucl_num);
+      assert(rc == QMCKL_SUCCESS);
+      
+      rc = qmckl_context_touch(context);
+      assert(rc == QMCKL_SUCCESS);
+      
+      // Call the provided function
+      rc = qmckl_get_mo_basis_mo_value(context,&mo_output[0][0], point_num*mo_num);
+      assert(rc == QMCKL_SUCCESS);
+
+      // Accumulate derivative using finite-difference coefficients
+      for (int i = 0; i < point_num; i++) {
+        for (int j = 0; j < mo_num; j++) {
+          if (m == -4) {
+            finite_difference_force_mo_value[k][a][i][j] = 0.0;
+          }
+          finite_difference_force_mo_value[k][a][i][j] += coef[m + 4] * mo_output[i][j]/delta_x;
+        }
+      }
+    }
+    temp_coord[k+a*3] = nucleus_coord[k+3*a];
+  }
+}
+
+// Reset coordinates in the context
+rc = qmckl_set_nucleus_coord(context, 'N', temp_coord, 3*nucl_num);
+assert(rc == QMCKL_SUCCESS);
+
+rc = qmckl_context_touch(context);
+assert(rc == QMCKL_SUCCESS);
+    
+free(nucleus_coord);
+free(temp_coord);
+
+
+for (int j = 0; j < mo_num; j++){
+  for (int i = 0; i < point_num; i++){
+      for (int a = 0; a < nucl_num; a++) {
+        for (int k = 0; k < 3; k++){
+            //printf("k=%i a=%i i=%i j=%i\n", k, a, i, j);
+            printf("%.10f\t", finite_difference_force_mo_value[k][a][i][j]);
+            printf("%.10f\n", forces_mo_value[a][k][i][j]);
+            assert(fabs(finite_difference_force_mo_value[k][a][i][j] -  forces_mo_value[a][k][i][j]) < 1.e-10);
+
+        }
+    }
+  }
+}
+
+printf("OK\n");
+
+      #+end_src
+  
+
 * End of files                                                     :noexport:
 
   #+begin_src c :tangle (eval h_private_type)

org/qmckl_jastrow_champ_single.org

@@ -4610,7 +4610,7 @@ qmckl_exit_code qmckl_provide_jastrow_champ_single_een_gl(qmckl_context context)
      | ~een_rescaled_single_n~    | ~double[walk_num][0:cord_num][nucl_num]~                            | in     | Electron-nucleus single rescaled distances |
      | ~een_rescaled_n_gl~        | ~double[walk_num][0:cord_num][nucl_num][4][elec_num]~               | in     | Electron-nucleus rescaled distances        |
      | ~een_rescaled_single_n_gl~ | ~double[walk_num][0:cord_num][nucl_num][4]~                         | in     | Electron-nucleus single rescaled distances |
-     | ~delta_een_gl~             | ~double[walk_num][elec_num][4]~                                     | out    | Electron-nucleus jastrow                   |
+     | ~delta_een_gl~             | ~double[walk_num][4][elec_num]~                                     | out    | Electron-nucleus jastrow                   |
 
      #+begin_src f90 :comments org :tangle (eval f) :noweb yes
 integer(qmckl_exit_code) function qmckl_compute_jastrow_champ_factor_single_een_gl_doc( &
@@ -5727,7 +5727,6 @@ function qmckl_compute_jastrow_champ_single_ee_gl_doc( &
 
   do nw =1, walk_num
      delta_ee_gl(:,:,nw) = 0.0d0
-
         do i = 1, elec_num
            if (i == num) cycle
 
@@ -5775,6 +5774,13 @@ function qmckl_compute_jastrow_champ_single_ee_gl_doc( &
                 + f * (dx(4) - 2.d0 * b_vector(2) * grad_c2 * invdenom) &
                 - f_old * (dx_old(4) - 2.d0 * b_vector(2) * grad_c2_old * invdenom_old)
 
+           delta_ee_gl(1,num,nw) = delta_ee_gl(1,num,nw) - f * dx(1) + f_old * dx_old(1)
+           delta_ee_gl(2,num,nw) = delta_ee_gl(2,num,nw) - f * dx(2) + f_old * dx_old(2)
+           delta_ee_gl(3,num,nw) = delta_ee_gl(3,num,nw) - f * dx(3) + f_old * dx_old(3)
+           delta_ee_gl(4,num,nw) = delta_ee_gl(4,num,nw) &
+                + f * (dx(4) - 2.d0 * b_vector(2) * grad_c2 * invdenom) &
+                - f_old * (dx_old(4) - 2.d0 * b_vector(2) * grad_c2_old * invdenom_old)
+
 
            kf = 2.d0
            x1 = x
@@ -5790,6 +5796,13 @@ function qmckl_compute_jastrow_champ_single_ee_gl_doc( &
               delta_ee_gl(4,i,nw) = delta_ee_gl(4,i,nw) &
                    + f * (x1 * dx(4) + (kf-1.d0) * grad_c2) &
                    - f_old * (x1_old * dx_old(4) + (kf-1.d0) * grad_c2_old)
+
+              delta_ee_gl(1,num,nw) = delta_ee_gl(1,num,nw) - f * x1 * dx(1) + f_old * x1_old * dx_old(1)
+              delta_ee_gl(2,num,nw) = delta_ee_gl(2,num,nw) - f * x1 * dx(2) + f_old * x1_old * dx_old(2)
+              delta_ee_gl(3,num,nw) = delta_ee_gl(3,num,nw) - f * x1 * dx(3) + f_old * x1_old * dx_old(3)
+              delta_ee_gl(4,num,nw) = delta_ee_gl(4,num,nw) &
+                   + f * (x1 * dx(4) + (kf-1.d0) * grad_c2) &
+                   - f_old * (x1_old * dx_old(4) + (kf-1.d0) * grad_c2_old)
               x = x*x1
               x_old = x_old*x1_old
               kf = kf + 1.d0
@@ -5909,7 +5922,6 @@ assert (rc == QMCKL_SUCCESS);
 for (int nw = 0; nw < walk_num; nw++) {
   for (int i = 0; i < elec_num; i++) {
       for (int m = 0; m < 4; m++) {
-        if (i == 2) continue;
         //printf("%f\n",(ee_gl_new[nw][m][i] - ee_gl_old[nw][m][i]));
         //printf("%f\n",delta_ee_gl[nw][i][m]);
         assert(fabs((ee_gl_new[nw][m][i] - ee_gl_old[nw][m][i]) - delta_ee_gl[nw][i][m]) < 1.e-12);