From 51e41b2b3ce269b4ee35cf14369dfb89f7ee05c6 Mon Sep 17 00:00:00 2001
From: Emiel Slootman <e.slootman@utwente.nl>
Date: Wed, 5 Feb 2025 12:44:32 +0100
Subject: [PATCH] optimizations

---
 org/qmckl_forces.org               | 198 +++++++++++-----
 org/qmckl_jastrow_champ_single.org | 365 ++++++++++++++++++++++++++++-
 2 files changed, 493 insertions(+), 70 deletions(-)

diff --git a/org/qmckl_forces.org b/org/qmckl_forces.org
index 449f107..8a3c0db 100644
--- a/org/qmckl_forces.org
+++ b/org/qmckl_forces.org
@@ -2142,20 +2142,19 @@ integer(qmckl_exit_code) function qmckl_compute_forces_dtmp_c( &
   if (cord_num <  0)                 info = QMCKL_INVALID_ARG_5
   if (info /= QMCKL_SUCCESS)         return
 
-  do nw=1, walk_num
-    do i=0, cord_num-1
-    do k =1, 4
-      info = qmckl_dgemm(context,'N','N',elec_num*1_8,&
-            nucl_num*(cord_num+1)*4_8, elec_num*1_8, -1.0d0,     &
-            een_rescaled_e_gl(1,k,1,i,nw),4_8*elec_num,              &
-            een_rescaled_n_gl(1,1,1,0,nw),elec_num*1_8,         &
-            0.0d0,                                              &
-            forces_dtmp_c(1,k,1,1,0,i,nw),4_8*elec_num)
-      end do
-    end do
-  end do
+  !do nw=1, walk_num
+  !  do i=0, cord_num-1
+  !  do k =1, 4
+  !    info = qmckl_dgemm(context,'N','N',elec_num*1_8,&
+  !          nucl_num*(cord_num+1)*4_8, elec_num*1_8, -1.0d0,     &
+  !          een_rescaled_e_gl(1,k,1,i,nw),4_8*elec_num,              &
+  !          een_rescaled_n_gl(1,1,1,0,nw),elec_num*1_8,         &
+  !          0.0d0,                                              &
+  !          forces_dtmp_c(1,k,1,1,0,i,nw),4_8*elec_num)
+  !    end do
+  !  end do
+  !end do
 
-  if (.true.) then
   do nw = 1, walk_num
     do l = 0, cord_num-1
       do m = 0, cord_num
@@ -2175,7 +2174,6 @@ integer(qmckl_exit_code) function qmckl_compute_forces_dtmp_c( &
       end do
     end do
   end do
-  end if
 
 end function qmckl_compute_forces_dtmp_c
      #+end_src
@@ -2543,6 +2541,8 @@ integer(qmckl_exit_code) function qmckl_compute_forces_jastrow_een( &
   if (cord_num <  0)                 info = QMCKL_INVALID_ARG_5
   if (info /= QMCKL_SUCCESS)         return
 
+  forces_jastrow_een = 0.d0
+
   do nw =1, walk_num
     do n = 1, dim_c_vector
       l = lkpm_combined_index(n, 1)
@@ -2555,8 +2555,8 @@ integer(qmckl_exit_code) function qmckl_compute_forces_jastrow_een( &
         if(cn == 0.d0) cycle
           do j = 1, elec_num
             do ii = 1, 3
-              accu = een_rescaled_n(j,a,m,nw) * forces_tmp_c(j,ii,a,m+l,k,nw)
-              accu = accu - een_rescaled_n_gl(j,ii,a,m,nw) * tmp_c(j,a,m+l,k,nw)
+              accu = een_rescaled_n(j,a,m,nw) * forces_tmp_c(j,ii,a,m+l,k,nw) & 
+                  - een_rescaled_n_gl(j,ii,a,m,nw) * tmp_c(j,a,m+l,k,nw)
 
               forces_jastrow_een(ii, a, nw) = forces_jastrow_een(ii, a, nw) + accu * cn
             end do
@@ -2790,7 +2790,7 @@ integer(qmckl_exit_code) function qmckl_compute_forces_een_rescaled_n_gl( &
   real    (c_double )   , intent(in)  :: een_rescaled_n_gl(elec_num,4,nucl_num,0:cord_num,walk_num)
   real    (c_double )   , intent(out) :: forces_een_n(elec_num,4,nucl_num,3,0:cord_num,walk_num)
 
-  double precision                    :: x, ria_inv, kappa_l
+  double precision                    :: x, ria_inv, kappa_l, temp
   integer*8                           :: i, a, k, l, nw, m, n
 
 
@@ -2823,6 +2823,34 @@ integer(qmckl_exit_code) function qmckl_compute_forces_een_rescaled_n_gl( &
 
 
 
+  !forces_een_n = 0.0d0
+  !do nw = 1, walk_num
+  !  do l = 0, cord_num
+  !    do a = 1, nucl_num
+  !      kappa_l = dble(l)*rescale_factor_en(type_nucl_vector(a)+1)
+  !      do i = 1, elec_num
+  !        ria_inv = 1.0d0 / en_distance(a,i,nw)
+  !        do n = 1, 3
+  !          do m = 1, 3
+  !            if (m == n) then
+  !              forces_een_n(i,m,a,n,l,nw) = forces_een_n(i,m,a,n,l,nw) + &
+  !                    kappa_l*een_rescaled_n(i,a,l,nw) * ria_inv
+  !            end if
+  !            temp = een_rescaled_n_gl(i,m,a,l,nw) * een_rescaled_n_gl(i,n,a,l,nw) / een_rescaled_n(i,a,l,nw)
+  !            forces_een_n(i,m,a,n,l,nw) = forces_een_n(i,m,a,n,l,nw) - &
+  !               temp * ria_inv / kappa_l - temp
+  !          end do
+  !          forces_een_n(i,4,a,n,l,nw) = forces_een_n(i,4,a,n,l,nw) + &
+  !            (2.0d0 * ria_inv * ria_inv & 
+  !            - een_rescaled_n_gl(i,4,a,l,nw)/een_rescaled_n(i,a,l,nw)) * & 
+  !            een_rescaled_n_gl(i,n,a,l,nw)
+  !        end do
+  !      end do
+  !    end do
+  !  end do
+  !end do
+
+
     forces_een_n = 0.0d0
     do nw = 1, walk_num
     do i = 1, elec_num
@@ -4345,7 +4373,6 @@ integer(qmckl_exit_code) function qmckl_compute_forces_jastrow_delta_p_doc( &
               do a = 1, nucl_num
                   delta_n_gl = een_rescaled_single_n_gl(k,a,l,nw) - een_rescaled_n_gl(num, k, a, l, nw)
                  do j = 1, elec_num
-
                     forces_delta_p(j,k,a,l,m,nw) =  delta_e(j) * (-een_rescaled_single_n_gl(k,a,l,nw)) - &
                          een_rescaled_e(j,num,m,nw) * delta_n_gl
                  end do
@@ -5308,6 +5335,7 @@ qmckl_exit_code qmckl_provide_forces_ao_value(qmckl_context context)
       }
       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);
@@ -5430,6 +5458,8 @@ function qmckl_compute_forces_ao_value_doc(context, &
      end do
   end do
 
+  deallocate(ao_index)
+
 end function qmckl_compute_forces_ao_value_doc
      #+end_src
 
@@ -5729,26 +5759,41 @@ integer(qmckl_exit_code)  function qmckl_compute_forces_mo_value_doc(context, &
   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
+  integer*8 :: i,j,k,a, m
   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
+  !do a=1,nucl_num
+  !  do j=1,point_num
+  !    do k=1,ao_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
+
+    
+  !do a=1,nucl_num
+  !  do m=1,3
+  !    info = qmckl_dgemm(context,'N', 'N', mo_num, point_num, ao_num, &
+  !               1.0d0, coefficient_t, mo_num, &
+  !               forces_ao_value(:, :, m, a), ao_num, &
+  !               0.0d0, forces_mo_value(:, :, m, a), mo_num)
+  !  end do
+  !end do
+
+      info = qmckl_dgemm(context,'N', 'N', mo_num, 3*nucl_num*point_num, ao_num, &
+                 1.0d0, coefficient_t(1,1), mo_num, &
+                 forces_ao_value(1,1,1,1), ao_num, &
+                 0.0d0, forces_mo_value(1,1,1,1), mo_num)
 
 
 end function qmckl_compute_forces_mo_value_doc
@@ -6062,19 +6107,48 @@ integer function qmckl_compute_forces_mo_g_doc(context, &
 
   info = QMCKL_SUCCESS
 
-  do j=1,point_num
-     forces_mo_g(:,:,j,:,:) = 0.d0
-     do k=1,ao_num
-        do i=1,mo_num
-          do a=1, nucl_num
-            do m = 1, 3 
-              do n = 1, 3
-               forces_mo_g(i, m, j, n, a) = forces_mo_g(i, m, j, n, a) - coefficient_t(i,k) * ao_hessian(k, m, j, n, a)
-              end do
-            end do
-          end do
-        end do
-     end do
+
+  forces_mo_g = 0.d0
+
+  !do j=1,point_num
+  !   !forces_mo_g(:,:,j,:,:) = 0.d0
+  !   do k=1,ao_num
+  !      do i=1,mo_num
+  !        do a=1, nucl_num
+  !          do m = 1, 3 
+  !            do n = 1, 3
+  !             forces_mo_g(i, m, j, n, a) = forces_mo_g(i, m, j, n, a) - coefficient_t(i,k) * ao_hessian(k, m, j, n, a)
+  !            end do
+  !          end do
+  !        end do
+  !      end do
+  !   end do
+  !end do
+
+
+
+    ! do k=1,ao_num
+    !    do i=1,mo_num
+    !        do m = 1, 3 
+    !           forces_mo_g(i, m, :,:,:) = forces_mo_g(i, m, :,:,:) - coefficient_t(i,k) * ao_hessian(k, m, :,:,:)
+    !          end do
+    !        end do
+    !      end do
+  !do m = 1, 3
+  !   info = qmckl_dgemm(context,'N', 'N', mo_num, point_num * 3 * nucl_num, ao_num, &
+  !              -1.d0, coefficient_t, mo_num, &
+  !              ao_hessian(:, m, :, :, :), ao_num, &
+  !              1.d0, forces_mo_g(:, m, :, :, :), mo_num)
+  !end do
+  do a=1,nucl_num
+     do n = 1,3
+    do m = 1, 3
+     info = qmckl_dgemm(context,'N', 'N', mo_num, point_num, ao_num, &
+                -1.d0, coefficient_t, mo_num, &
+                ao_hessian(:, m, :, n, a), ao_num, &
+                1.d0, forces_mo_g(:, m, :, n, a), mo_num)
+      end do
+  end do
   end do
 
 
@@ -6392,21 +6466,31 @@ integer function qmckl_compute_forces_mo_l_doc(context, &
   integer*8 :: i,j,k, m, n,a
 
   info = QMCKL_SUCCESS
+  forces_mo_l = 0.d0
+  !do j=1,point_num
+  !   forces_mo_l(:,j,:,:) = 0.d0
+  !   do k=1,ao_num
+  !      do i=1,mo_num
+  !        do a=1, nucl_num
+  !          do n = 1, 3
+  !            forces_mo_l(i, j, n, a) = forces_mo_l(i, j, n, a) - coefficient_t(i,k) * ao_hessian(k, 4, j, n, a)
+  !          end do
+  !        end do
+  !      end do
+  !   end do
+  !end do
 
-  do j=1,point_num
-     forces_mo_l(:,j,:,:) = 0.d0
-     do k=1,ao_num
-        do i=1,mo_num
-          do a=1, nucl_num
-            do n = 1, 3
-              forces_mo_l(i, j, n, a) = forces_mo_l(i, j, n, a) - coefficient_t(i,k) * ao_hessian(k, 4, j, n, a)
-            end do
-          end do
-        end do
-     end do
+  do a=1,nucl_num
+    do m = 1, 3
+     info = qmckl_dgemm(context,'N', 'N', mo_num, point_num, ao_num, &
+                -1.d0, coefficient_t, mo_num, &
+                ao_hessian(:, 4, :, m, a), ao_num, &
+                1.d0, forces_mo_l(:, :, m, a), mo_num)
+    end do
   end do
 
 
+
 end function qmckl_compute_forces_mo_l_doc
     #+end_src
 
diff --git a/org/qmckl_jastrow_champ_single.org b/org/qmckl_jastrow_champ_single.org
index da44d2f..5a0f8be 100644
--- a/org/qmckl_jastrow_champ_single.org
+++ b/org/qmckl_jastrow_champ_single.org
@@ -3291,12 +3291,8 @@ integer(qmckl_exit_code) function qmckl_compute_jastrow_champ_delta_p_gl_doc( &
            do k = 1, 4
               do a = 1, nucl_num
                  een_rescaled_delta_n = een_rescaled_single_n(a,l,nw) - een_rescaled_n(num, a, l, nw)
-                 !tmp = fk(k) * (een_rescaled_n(num,a,l,nw) + een_rescaled_delta_n)
-                 tmp = fk(k) * een_rescaled_single_n(a,l,nw)
-                 do j = 1, elec_num
-                    delta_p_gl(j,a,k,l,m,nw) =  delta_e_gl(j,k) * tmp  + &
-                         een_rescaled_e_gl(j,k,num,m,nw) * een_rescaled_delta_n
-                 end do
+                 delta_p_gl(:,a,k,l,m,nw) = delta_e_gl(:,k) * fk(k) * een_rescaled_single_n(a,l,nw)  + &
+                      een_rescaled_e_gl(:,k,num,m,nw) * een_rescaled_delta_n
                  do j = 1, elec_num
                     delta_p_gl(num,a,k,l,m,nw) = delta_p_gl(num,a,k,l,m,nw) + delta_e_gl(j,k) * een_rescaled_n(j,a,l,nw)
                  end do
@@ -3699,6 +3695,7 @@ integer(qmckl_exit_code) function qmckl_compute_jastrow_champ_factor_single_een_
                       delta_p(i,a,m  ,k,nw) * een_rescaled_n_gl(i,kk,a,m+l,nw) + &
                       delta_p(i,a,m+l,k,nw) * een_rescaled_n_gl(i,kk,a,m  ,nw) ) * cn
               end do
+
               delta_een_gl(num,kk,nw) = delta_een_gl(num,kk,nw) + ( &
                    (dtmp_c(num,kk,a,m  ,k,nw) + delta_p_gl(num,a,kk,m  ,k,nw)) * een_rescaled_delta_n(a,m+l,nw) + &
                    (dtmp_c(num,kk,a,m+l,k,nw) + delta_p_gl(num,a,kk,m+l,k,nw)) * een_rescaled_delta_n(a,m  ,nw) + &
@@ -3903,12 +3900,12 @@ qmckl_get_jastrow_champ_single_een_g(qmckl_context context,
                            "Array is NULL.");
   }
 
-  int64_t sze = 3 * ctx->electron.num * ctx->electron.walker.num;
+  int64_t sze = 4 * ctx->electron.num * ctx->electron.walker.num;
   if (size_max < sze) {
     return qmckl_failwith( context,
                            QMCKL_INVALID_ARG_3,
                            "qmckl_get_jastrow_champ_single_een_g",
-                           "Array too small. Expected 3 * ctx->electron.num * ctx->electron.walker.num");
+                           "Array too small. Expected 4 * ctx->electron.num * ctx->electron.walker.num");
   }
   memcpy(delta_een_g, ctx->single_point.delta_een_g, sze * sizeof(double));
 
@@ -3985,7 +3982,7 @@ qmckl_exit_code qmckl_provide_jastrow_champ_single_een_g(qmckl_context context)
     if (ctx->single_point.delta_een_g == NULL) {
 
       qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
-      mem_info.size = 3 * ctx->electron.num * ctx->electron.walker.num * sizeof(double);
+      mem_info.size = 4 * ctx->electron.num * ctx->electron.walker.num * sizeof(double);
       double* delta_een_g = (double*) qmckl_malloc(context, mem_info);
 
       if (delta_een_g == NULL) {
@@ -4053,7 +4050,7 @@ qmckl_exit_code qmckl_provide_jastrow_champ_single_een_g(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 derivatives                |
      | ~een_rescaled_single_n_gl~ | ~double[walk_num][0:cord_num][nucl_num][4]~                         | in     | Electron-nucleus single rescaled distances derivatives         |
-     | ~delta_een_g~              | ~double[walk_num][3][elec_num]~                                     | out    | Delta electron-electron-nucleus jastrow gradient               |
+     | ~delta_een_g~              | ~double[walk_num][4][elec_num]~                                     | out    | Delta electron-electron-nucleus jastrow gradient               |
      |----------------------------+---------------------------------------------------------------------+--------+----------------------------------------------------------------|
 
      #+begin_src f90 :comments org :tangle (eval f) :noweb yes
@@ -4078,7 +4075,7 @@ integer(qmckl_exit_code) function qmckl_compute_jastrow_champ_factor_single_een_
   real(c_double)        , intent(in)  :: een_rescaled_single_n(nucl_num, 0:cord_num, walk_num)
   real(c_double)        , intent(in)  :: een_rescaled_n_gl(elec_num, 4, nucl_num, 0:cord_num, walk_num)
   real(c_double)        , intent(in)  :: een_rescaled_single_n_gl(4, nucl_num, 0:cord_num, walk_num)
-  real(c_double)        , intent(out) :: delta_een_g(elec_num, 3, walk_num)
+  real(c_double)        , intent(out) :: delta_een_g(elec_num, 4, walk_num)
 
   integer*8 :: i, a, j, l, k, p, m, n, nw, kk, num
   double precision :: accu, accu2, cn
@@ -4222,7 +4219,7 @@ qmckl_compute_jastrow_champ_factor_single_een_g (const qmckl_context context,
 assert(qmckl_jastrow_champ_provided(context));
 
 double een_g_old[walk_num][4][elec_num];
-double delta_een_g[walk_num][3][elec_num];
+double delta_een_g[walk_num][4][elec_num];
 double een_g_new[walk_num][4][elec_num];
 
 
@@ -4240,7 +4237,7 @@ for (int elec = 0; elec < elec_num; elec++){
   rc = qmckl_set_single_point(context, 'N', elec, new_coords, 3*walk_num);
   assert (rc == QMCKL_SUCCESS);
 
-  rc = qmckl_get_jastrow_champ_single_een_g(context, &delta_een_g[0][0][0], walk_num*3*elec_num);
+  rc = qmckl_get_jastrow_champ_single_een_g(context, &delta_een_g[0][0][0], walk_num*4*elec_num);
   assert (rc == QMCKL_SUCCESS);
 
   coords[0][elec][0] = new_coords[0];
@@ -7122,6 +7119,348 @@ qmckl_get_jastrow_champ_single_accept(qmckl_context context);
 
      #+begin_src c :comments org :tangle (eval c) :noweb yes  
 qmckl_exit_code
+qmckl_get_jastrow_champ_single_accept_alt(qmckl_context context)
+{
+
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return QMCKL_NULL_CONTEXT;
+  }
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
+  assert (ctx != NULL);
+
+  qmckl_exit_code rc;
+
+  double metric[4] = {-1.0, -1.0, -1.0, 1.0};
+
+  int old_date, do_update;
+
+  int shift1, shift2, shift3, shift4, shift5, shift6, shift7;
+
+  int shift8, shift9, shift10, shift11, shift12, shift13, shift14;
+
+  old_date = ctx->single_point.date;
+  do_update = 0;
+
+  rc = qmckl_context_touch(context);
+  if (rc != QMCKL_SUCCESS) return rc;
+
+  ctx->single_point.date = ctx->date;
+
+  if (ctx->jastrow_champ.cord_num > 0) {
+    if(old_date == ctx->single_point.delta_een_date) {
+
+      shift1 = (ctx->jastrow_champ.cord_num+1)*ctx->electron.num*ctx->electron.num;
+      shift2 = ctx->electron.num*ctx->electron.num;
+      shift3 = (ctx->jastrow_champ.cord_num+1)*ctx->nucleus.num*ctx->electron.num;
+      shift4 = ctx->nucleus.num*ctx->electron.num;
+      shift5 = (ctx->jastrow_champ.cord_num+1)*ctx->nucleus.num;
+
+      for (int nw = 0; nw < ctx->electron.walker.num; nw++) {
+
+        ctx->jastrow_champ.factor_een[nw] = ctx->jastrow_champ.factor_een[nw] + ctx->single_point.delta_een[nw];
+
+        for (int l = 0; l <= ctx->jastrow_champ.cord_num; l++){
+          for (int i = 0; i < ctx->electron.num; i++) {
+            ctx->jastrow_champ.een_rescaled_e[nw*shift1
+              + l*shift2 
+              + i*ctx->electron.num 
+              + ctx->single_point.num] = 
+            ctx->single_point.een_rescaled_single_e[nw*shift3
+              + l*ctx->electron.num
+              + i];
+
+            ctx->jastrow_champ.een_rescaled_e[nw*shift1
+              + l*shift2
+              + ctx->single_point.num*ctx->electron.num 
+              + i] = 
+            ctx->single_point.een_rescaled_single_e[nw*shift3
+              + l*ctx->electron.num
+              + i];
+
+          }
+          for (int a = 0; a < ctx->nucleus.num; a++){
+            ctx->jastrow_champ.een_rescaled_n[nw*shift3
+              + l*shift4
+              + a*ctx->electron.num
+              + ctx->single_point.num] = 
+            ctx->single_point.een_rescaled_single_n[nw*shift5
+              + l*ctx->nucleus.num
+              + a];
+          }
+        }
+      }
+      for (int i = 0; i < (ctx->electron.walker.num*(ctx->jastrow_champ.cord_num+1)*ctx->nucleus.num*ctx->electron.num*ctx->jastrow_champ.cord_num); i++) {
+        ctx->jastrow_champ.tmp_c[i] = ctx->jastrow_champ.tmp_c[i] + ctx->single_point.delta_p[i];
+      }
+      ctx->jastrow_champ.een_rescaled_e_date = ctx->date;
+      ctx->jastrow_champ.een_rescaled_n_date = ctx->date;
+      ctx->jastrow_champ.factor_een_date = ctx->date;
+      ctx->jastrow_champ.tmp_c_date = ctx->date;
+      
+      ctx->single_point.een_rescaled_single_e_date = ctx->single_point.date;
+      ctx->single_point.een_rescaled_single_n_date = ctx->single_point.date;
+      ctx->single_point.delta_een_date = ctx->single_point.date;
+      ctx->single_point.delta_p_date = ctx->single_point.date;
+    }
+
+    if(old_date == ctx->single_point.delta_een_gl_date) {
+
+
+      for (int i = 0; i < ctx->electron.walker.num * 4 * ctx->electron.num; i++) {
+        ctx->jastrow_champ.factor_een_gl[i] = ctx->jastrow_champ.factor_een_gl[i] + ctx->single_point.delta_een_gl[i];
+      }  
+      
+      ctx->jastrow_champ.factor_een_gl_date = ctx->date;
+      ctx->single_point.delta_een_gl_date = ctx->single_point.date;
+      do_update = 1;
+    } else if (old_date == ctx->single_point.delta_een_g_date) {
+
+      for (int nw = 0; nw < ctx->electron.walker.num; nw++) {
+        for (int k = 0; k < 3; k++){
+          for (int i = 0; i < ctx->electron.num; i++){
+            ctx->jastrow_champ.factor_een_gl[nw*ctx->electron.num*4 + k*ctx->electron.num + i] = 
+            ctx->jastrow_champ.factor_een_gl[nw*ctx->electron.num*4 + k*ctx->electron.num + i] + 
+            ctx->single_point.delta_een_g[nw*ctx->electron.num*3 + k*ctx->electron.num + i];
+          }
+        }
+      }  
+      ctx->jastrow_champ.factor_een_gl_date = ctx->date;
+      ctx->single_point.delta_een_g_date = ctx->single_point.date;
+      do_update = 1;
+    }
+
+    if (do_update == 1) {
+
+      shift1 = (ctx->jastrow_champ.cord_num+1)*ctx->electron.num*4*ctx->electron.num;
+      shift2 = ctx->electron.num*4*ctx->electron.num;
+      shift3 = ctx->electron.num*4;
+      shift4 = (ctx->jastrow_champ.cord_num+1)*ctx->electron.num*4;
+      shift5 = (ctx->jastrow_champ.cord_num+1)*ctx->nucleus.num*4*ctx->electron.num;
+      shift6 = ctx->nucleus.num*4*ctx->electron.num;
+      shift7 = (ctx->jastrow_champ.cord_num+1)*ctx->nucleus.num*4;
+      shift8 = ctx->nucleus.num*4;
+      shift9 = ctx->nucleus.num*ctx->electron.num;
+
+      for (int nw = 0; nw < ctx->electron.walker.num; nw++) {
+        for (int l = 0; l <= ctx->jastrow_champ.cord_num; l++){
+          for (int i = 0; i < ctx->electron.num; i++) {
+            for (int k = 0; k < 4; k++){
+              ctx->jastrow_champ.een_rescaled_e_gl[nw*shift1
+                + l*shift2
+                + i*shift3
+                + k*ctx->electron.num 
+                + ctx->single_point.num] = 
+              ctx->single_point.een_rescaled_single_e_gl[nw*shift4
+                + l*shift3
+                + i*4 
+                + k]; 
+              ctx->jastrow_champ.een_rescaled_e_gl[nw*shift1
+                + l*shift2
+                + ctx->single_point.num*shift3
+                + k*ctx->electron.num 
+                + i] = 
+              metric[k] * ctx->single_point.een_rescaled_single_e_gl[nw*shift4
+                + l*shift3
+                + i*4
+                + k]; 
+            }
+          }
+          for (int a = 0; a < ctx->nucleus.num; a++){
+            for (int k = 0; k < 4; k++){
+              ctx->jastrow_champ.een_rescaled_n_gl[nw*shift5
+                + l*shift6
+                + a*shift3
+                + k*ctx->electron.num
+                + ctx->single_point.num] = 
+              ctx->single_point.een_rescaled_single_n_gl[nw*shift7
+                + l*shift8
+                + a*4
+                + k];
+            }
+          }
+        }
+      }
+      for (int nw = 0; nw < ctx->electron.walker.num*(ctx->jastrow_champ.cord_num+1)*ctx->jastrow_champ.cord_num; nw++) {
+        for (int a = 0; a < ctx->nucleus.num; a++) {
+          for (int k = 0; k < 4; k++){
+            for (int i = 0; i < ctx->electron.num; i++) {
+              ctx->jastrow_champ.dtmp_c[nw*shift6
+                + a*shift3
+                + k*ctx->electron.num
+                + i] =
+              ctx->jastrow_champ.dtmp_c[nw*shift6
+                + a*shift3
+                + k*ctx->electron.num
+                + i] +
+              ctx->single_point.delta_p_gl[nw*shift6
+                + k*shift9
+                + a*ctx->electron.num
+                + i];
+            }
+          }
+        }
+      }
+      ctx->jastrow_champ.dtmp_c_date = ctx->date;
+      ctx->jastrow_champ.een_rescaled_e_gl_date = ctx->date;
+      ctx->jastrow_champ.een_rescaled_n_gl_date = ctx->date;
+
+      ctx->single_point.een_rescaled_single_e_gl_date = ctx->single_point.date;
+      ctx->single_point.een_rescaled_single_n_gl_date = ctx->single_point.date;
+      ctx->single_point.delta_p_gl_date = ctx->single_point.date;
+    }
+  }
+  
+
+    shift1 = (ctx->jastrow_champ.cord_num+1)*ctx->electron.num*ctx->electron.num;
+    shift2 = ctx->electron.num*ctx->electron.num;
+    shift3 = (ctx->jastrow_champ.cord_num+1)*ctx->electron.num;
+    shift4 = (ctx->jastrow_champ.cord_num+1)*ctx->electron.num*4*ctx->electron.num;
+    shift5 = ctx->electron.num*4*ctx->electron.num;
+    shift6 = ctx->electron.num*4;
+    shift7 = (ctx->jastrow_champ.cord_num+1)*ctx->electron.num*4;
+    shift8 = (ctx->jastrow_champ.cord_num+1)*ctx->nucleus.num*ctx->electron.num;
+    shift9 = (ctx->jastrow_champ.cord_num+1)*ctx->nucleus.num;
+    shift10 = (ctx->jastrow_champ.cord_num+1)*ctx->nucleus.num*4*ctx->electron.num;
+    shift11 = ctx->nucleus.num*4*ctx->electron.num;
+    shift12 = (ctx->jastrow_champ.cord_num+1)*ctx->nucleus.num*4;
+    shift13 = ctx->nucleus.num*4;
+    shift14 = ctx->nucleus.num*ctx->electron.num;
+
+
+  for (int nw = 0; nw < ctx->electron.walker.num; nw++) {
+    ctx->jastrow_champ.factor_en[nw] = ctx->jastrow_champ.factor_en[nw] + ctx->single_point.delta_en[nw];
+    ctx->jastrow_champ.factor_ee[nw] = ctx->jastrow_champ.factor_ee[nw] + ctx->single_point.delta_ee[nw];
+    for (int a = 0; a < ctx->nucleus.num; a++) {
+      ctx->electron.en_distance[nw*shift14
+        + ctx->single_point.num*ctx->nucleus.num
+        + a] =
+      ctx->single_point.single_en_distance[nw*ctx->nucleus.num
+        + a];
+
+      ctx->jastrow_champ.en_distance_rescaled[nw*shift14
+        + a*ctx->electron.num
+        + ctx->single_point.num] =
+      ctx->single_point.en_rescaled_single[nw*ctx->nucleus.num
+        + a];
+      for (int k = 0; k < 4; k++){
+        ctx->jastrow_champ.en_distance_rescaled_gl[nw*shift11
+          + a*shift6
+          + ctx->single_point.num*4
+          + k] = 
+        ctx->single_point.en_rescaled_single_gl[nw*shift13
+          + a*4
+          + k];
+      }
+    }
+    for (int i = 0; i < ctx->electron.num; i++) {
+      ctx->jastrow_champ.ee_distance_rescaled[nw*shift2
+        + i*ctx->electron.num
+        + ctx->single_point.num] =
+      ctx->single_point.ee_rescaled_single[nw*ctx->electron.num
+        + i];
+
+      ctx->jastrow_champ.ee_distance_rescaled[nw*shift2
+        + ctx->single_point.num*ctx->electron.num
+        + i] =
+      ctx->single_point.ee_rescaled_single[nw*ctx->electron.num
+        + i];
+
+      ctx->electron.ee_distance[nw*shift2
+        + i*ctx->electron.num
+        + ctx->single_point.num] =
+      ctx->single_point.single_ee_distance[nw*ctx->electron.num
+        + i];
+
+      ctx->electron.ee_distance[nw*shift2
+        + ctx->single_point.num*ctx->electron.num
+        + i] =
+      ctx->single_point.single_ee_distance[nw*ctx->electron.num
+        + i];
+
+      for (int k = 0; k < 4; k++){
+        ctx->jastrow_champ.ee_distance_rescaled_gl[nw*shift5
+          + i*shift6
+          + ctx->single_point.num*4 
+          + k] = 
+        metric[k] * ctx->single_point.ee_rescaled_single_gl[nw*shift6
+          + i*4 
+          + k];
+        ctx->jastrow_champ.ee_distance_rescaled_gl[nw*shift5
+          + ctx->single_point.num*shift6
+          + i*4 
+          + k] =  
+        ctx->single_point.ee_rescaled_single_gl[nw*shift6
+          + i*4 
+          + k];
+      }
+    }
+    for (int k = 0; k < 4; k++){
+      for (int i = 0; i < ctx->electron.num; i++) {
+
+        ctx->jastrow_champ.factor_ee_gl[nw*shift6
+          + k*ctx->electron.num
+          + i] = 
+        ctx->jastrow_champ.factor_ee_gl[nw*shift6
+          + k*ctx->electron.num
+          + i] + 
+        ctx->single_point.delta_ee_gl[nw*shift6
+          + i*4
+          + k];
+
+        ctx->jastrow_champ.factor_en_gl[nw*shift6
+          + k*ctx->electron.num
+          + i] = 
+        ctx->jastrow_champ.factor_en_gl[nw*shift6
+          + k*ctx->electron.num
+          + i] + 
+        ctx->single_point.delta_en_gl[nw*shift6
+          + i*4
+          + k];
+      }
+    }
+  }
+ 
+  for (int nw = 0; nw < ctx->electron.walker.num; nw++) {
+    for (int k = 0; k < 3; k++) {
+      ctx->point.coord.data[nw*3*ctx->electron.num + k*ctx->electron.num + ctx->single_point.num] = ctx->single_point.coord.data[nw*3 + k];
+    }
+  }
+
+
+  ctx->jastrow_champ.ee_distance_rescaled_date = ctx->date;
+  ctx->jastrow_champ.ee_distance_rescaled_gl_date = ctx->date;
+  ctx->jastrow_champ.en_distance_rescaled_date = ctx->date;
+  ctx->jastrow_champ.en_distance_rescaled_gl_date = ctx->date;
+  ctx->jastrow_champ.factor_ee_date = ctx->date;
+  ctx->jastrow_champ.factor_ee_gl_date = ctx->date;
+  ctx->jastrow_champ.factor_en_date = ctx->date;
+  ctx->jastrow_champ.factor_en_gl_date = ctx->date;
+
+  ctx->electron.ee_distance_date = ctx->date;
+  ctx->electron.en_distance_date = ctx->date;
+  
+  ctx->single_point.date = ctx->date;
+
+
+  ctx->single_point.single_ee_distance_date = ctx->single_point.date;
+  ctx->single_point.single_en_distance_date = ctx->single_point.date;
+  ctx->single_point.ee_rescaled_single_date = ctx->single_point.date;
+  ctx->single_point.en_rescaled_single_date = ctx->single_point.date;
+  ctx->single_point.delta_en_date = ctx->single_point.date;
+  ctx->single_point.delta_ee_date = ctx->single_point.date;
+  ctx->single_point.ee_rescaled_single_gl_date = ctx->single_point.date;
+  ctx->single_point.en_rescaled_single_gl_date = ctx->single_point.date;
+  ctx->single_point.delta_en_gl_date = ctx->single_point.date;
+  ctx->single_point.delta_ee_gl_date = ctx->single_point.date;
+
+
+  return QMCKL_SUCCESS;
+}
+     #+end_src
+
+     #+begin_src c :comments org :tangle (eval c) :noweb yes  
+qmckl_exit_code
 qmckl_get_jastrow_champ_single_accept(qmckl_context context)
 {