Merge pull request #16 from v1j4y/rescaled_vj

Rescaled Distances

org/qmckl.org

@@ -208,11 +208,11 @@ cppcheck --addon=cert --enable=all *.c &> cppcheck.out
    value in the ~qmckl_numprec_struct~ contained in the context.
    Because of these intricate dependencies, a private header is
    created, containing the ~qmckl_numprec_struct~. This header is
-   included in the private header which defines the type of the
-   context. Headers for private types are suffixed by =_private_type.h=
-   and headers for private functions, =_private_func.h=.
-   Fortran interfaces should also be written in the =*_f_func.f90= file,
-   and the types definitions should be written in the =*_f_type.f90= file.
+   included in the private header file which defines the type of the
+   context. Header files for private types are suffixed by =_private_type.h=
+   and header files for private functions are suffixed by =_private_func.h=.
+   Fortran interfaces should also be written in the =*fh_func.f90= file,
+   and the types definitions should be written in the =*fh_type.f90= file.
 
    | File               | Scope   | Description                  |
    |--------------------+---------+------------------------------|
@@ -245,6 +245,7 @@ cppcheck --addon=cert --enable=all *.c &> cppcheck.out
    the =context= variable.
 
    # TODO : We need an identifier for impure functions
+   # Suggestion (VJ): using *_unsafe_* for impure functions ?
 
 ** Numerical precision
 

org/qmckl_context.org

@@ -117,7 +117,6 @@ typedef struct qmckl_context_struct {
   qmckl_ao_basis_struct    ao_basis;
 
   /* To be implemented:
-  qmckl_nucleus_struct     nucleus;
   qmckl_mo_struct          mo;
   qmckl_determinant_struct det;
   ,*/
@@ -174,6 +173,7 @@ qmckl_context qmckl_context_check(const qmckl_context context) {
    To create a new context, ~qmckl_context_create()~ should be used.
    - Upon success, it returns a pointer to a new context with the ~qmckl_context~ type
    - It returns ~QMCKL_NULL_CONTEXT~ upon failure to allocate the internal data structure
+   - A new context always has all its members initialized with a NULL value
 
    # Header
    #+begin_src c :comments org :tangle (eval h_func) :exports none
@@ -222,8 +222,8 @@ qmckl_context qmckl_context_create() {
   ctx->numprec.range = QMCKL_DEFAULT_RANGE;
 
   ctx->ao_basis.uninitialized = (1 << 10) - 1;
-  ctx->nucleus.uninitialized  = (1 << 3) - 1;
-  ctx->electron.uninitialized = (1 << 2) - 1;
+  ctx->nucleus.uninitialized  = (1 << 4) - 1;
+  ctx->electron.uninitialized = (1 << 3) - 1;
 
  /* Allocate qmckl_memory_struct */
   {

org/qmckl_distance.org

@@ -520,6 +520,7 @@ integer function qmckl_distance_f(context, transa, transb, m, n, &
      return
   endif
 
+  ! check for LDA
   if (iand(transab,1) == 0 .and. LDA < 3) then
      info = QMCKL_INVALID_ARG_7
      return
@@ -540,6 +541,33 @@ integer function qmckl_distance_f(context, transa, transb, m, n, &
      return
   endif
 
+  ! check for LDB
+  if (iand(transab,1) == 0 .and. LDB < 3) then
+     info = QMCKL_INVALID_ARG_9
+     return
+  endif
+
+  if (iand(transab,1) == 1 .and. LDB < n) then
+     info = QMCKL_INVALID_ARG_9
+     return
+  endif
+
+  if (iand(transab,2) == 0 .and. LDB < 3) then
+     info = QMCKL_INVALID_ARG_9
+     return
+  endif
+
+  if (iand(transab,2) == 2 .and. LDB < n) then
+     info = QMCKL_INVALID_ARG_9
+     return
+  endif
+
+  ! check for LDC
+  if (LDC < m) then
+     info = QMCKL_INVALID_ARG_11
+     return
+  endif
+
 
   select case (transab)
 
@@ -775,6 +803,308 @@ end function test_qmckl_dist
 qmckl_exit_code test_qmckl_dist(qmckl_context context);
 assert(test_qmckl_dist(context) == QMCKL_SUCCESS);
     #+end_src
+
+* Rescaled Distance
+
+** ~qmckl_distance_rescaled~
+   :PROPERTIES:
+   :Name:     qmckl_distance_rescaled
+   :CRetType: qmckl_exit_code
+   :FRetType: qmckl_exit_code
+   :END:
+
+ ~qmckl_distance_rescaled~ computes the matrix of the rescaled distances between all
+   pairs of points in two sets, one point within each set:
+
+   \[
+   C_{ij} = \left( 1 - \exp{-\kappa C_{ij}}\right)/\kappa
+   \]
+
+   If the input array is normal (~'N'~), the xyz coordinates are in
+   the leading dimension: ~[n][3]~ in C and ~(3,n)~ in Fortran.
+
+   #+NAME: qmckl_distance_rescaled_args
+   | qmckl_context | context              | in  | Global state                                  |
+   | char          | transa               | in  | Array ~A~ is ~'N'~: Normal, ~'T'~: Transposed |
+   | char          | transb               | in  | Array ~B~ is ~'N'~: Normal, ~'T'~: Transposed |
+   | int64_t       | m                    | in  | Number of points in the first set             |
+   | int64_t       | n                    | in  | Number of points in the second set            |
+   | double        | A[][lda]             | in  | Array containing the $m \times 3$ matrix $A$  |
+   | int64_t       | lda                  | in  | Leading dimension of array ~A~                |
+   | double        | B[][ldb]             | in  | Array containing the $n \times 3$ matrix $B$  |
+   | int64_t       | ldb                  | in  | Leading dimension of array ~B~                |
+   | double        | C[n][ldc]            | out | Array containing the $m \times n$ matrix $C$  |
+   | int64_t       | ldc                  | in  | Leading dimension of array ~C~                |
+   | double        | rescale_factor_kappa | in  | Factor for calculating rescaled distances     |
+
+*** Requirements
+
+    - ~context~ is not ~QMCKL_NULL_CONTEXT~
+    - ~m > 0~
+    - ~n > 0~
+    - ~lda >= 3~ if ~transa == 'N'~
+    - ~lda >= m~ if ~transa == 'T'~
+    - ~ldb >= 3~ if ~transb == 'N'~
+    - ~ldb >= n~ if ~transb == 'T'~
+    - ~ldc >= m~
+    - ~A~ is allocated with at least $3 \times m \times 8$ bytes
+    - ~B~ is allocated with at least $3 \times n \times 8$ bytes
+    - ~C~ is allocated with at least $m \times n \times 8$ bytes
+
+*** C header
+
+    #+CALL: generate_c_header(table=qmckl_distance_rescaled_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
+
+    #+RESULTS:
+    #+begin_src c :tangle (eval h_func) :comments org
+    qmckl_exit_code qmckl_distance_rescaled (
+          const qmckl_context context,
+          const char transa,
+          const char transb,
+          const int64_t m,
+          const int64_t n,
+          const double* A,
+          const int64_t lda,
+          const double* B,
+          const int64_t ldb,
+          double* const C,
+          const int64_t ldc,
+          const double  rescale_factor_kappa);
+    #+end_src
+
+*** Source
+    #+begin_src f90 :tangle (eval f)
+integer function qmckl_distance_rescaled_f(context, transa, transb, m, n, &
+     A, LDA, B, LDB, C, LDC, rescale_factor_kappa) &
+     result(info)
+  use qmckl
+  implicit none
+  integer(qmckl_context)  , intent(in)  :: context
+  character  , intent(in)  :: transa, transb
+  integer*8  , intent(in)  :: m, n
+  integer*8  , intent(in)  :: lda
+  real*8     , intent(in)  :: A(lda,*)
+  integer*8  , intent(in)  :: ldb
+  real*8     , intent(in)  :: B(ldb,*)
+  integer*8  , intent(in)  :: ldc
+  real*8     , intent(out) :: C(ldc,*)
+  real*8     , intent(in)  :: rescale_factor_kappa
+
+  integer*8 :: i,j
+  real*8    :: x, y, z, dist, rescale_factor_kappa_inv
+  integer   :: transab
+
+  rescale_factor_kappa_inv = 1.0d0/rescale_factor_kappa;
+
+  info = QMCKL_SUCCESS
+
+  if (context == QMCKL_NULL_CONTEXT) then
+     info = QMCKL_INVALID_CONTEXT
+     return
+  endif
+
+  if (m <= 0_8) then
+     info = QMCKL_INVALID_ARG_4
+     return
+  endif
+
+  if (n <= 0_8) then
+     info = QMCKL_INVALID_ARG_5
+     return
+  endif
+
+  if (transa == 'N' .or. transa == 'n') then
+     transab = 0
+  else if (transa == 'T' .or. transa == 't') then
+     transab = 1
+  else
+     transab = -100
+  endif
+
+  if (transb == 'N' .or. transb == 'n') then
+     continue
+  else if (transa == 'T' .or. transa == 't') then
+     transab = transab + 2
+  else
+     transab = -100
+  endif
+
+  ! check for LDA
+  if (transab < 0) then
+     info = QMCKL_INVALID_ARG_1
+     return
+  endif
+
+  if (iand(transab,1) == 0 .and. LDA < 3) then
+     info = QMCKL_INVALID_ARG_7
+     return
+  endif
+
+  if (iand(transab,1) == 1 .and. LDA < m) then
+     info = QMCKL_INVALID_ARG_7
+     return
+  endif
+
+  if (iand(transab,2) == 0 .and. LDA < 3) then
+     info = QMCKL_INVALID_ARG_7
+     return
+  endif
+
+  if (iand(transab,2) == 2 .and. LDA < m) then
+     info = QMCKL_INVALID_ARG_7
+     return
+  endif
+
+  ! check for LDB
+  if (iand(transab,1) == 0 .and. LDB < 3) then
+     info = QMCKL_INVALID_ARG_9
+     return
+  endif
+
+  if (iand(transab,1) == 1 .and. LDB < n) then
+     info = QMCKL_INVALID_ARG_9
+     return
+  endif
+
+  if (iand(transab,2) == 0 .and. LDB < 3) then
+     info = QMCKL_INVALID_ARG_9
+     return
+  endif
+
+  if (iand(transab,2) == 2 .and. LDB < n) then
+     info = QMCKL_INVALID_ARG_9
+     return
+  endif
+
+  ! check for LDC
+  if (LDC < m) then
+     info = QMCKL_INVALID_ARG_11
+     return
+  endif
+
+
+  select case (transab)
+
+  case(0)
+
+     do j=1,n
+        do i=1,m
+           x = A(1,i) - B(1,j)
+           y = A(2,i) - B(2,j)
+           z = A(3,i) - B(3,j)
+           dist = dsqrt(x*x + y*y + z*z)
+           C(i,j) = (1.0d0 - dexp(-rescale_factor_kappa * dist)) * rescale_factor_kappa_inv
+        end do
+     end do
+
+  case(1)
+
+     do j=1,n
+        do i=1,m
+           x = A(i,1) - B(1,j)
+           y = A(i,2) - B(2,j)
+           z = A(i,3) - B(3,j)
+           dist = dsqrt(x*x + y*y + z*z)
+           C(i,j) = (1.0d0 - dexp(-rescale_factor_kappa * dist)) * rescale_factor_kappa_inv
+        end do
+     end do
+
+  case(2)
+
+     do j=1,n
+        do i=1,m
+           x = A(1,i) - B(j,1)
+           y = A(2,i) - B(j,2)
+           z = A(3,i) - B(j,3)
+           dist = dsqrt(x*x + y*y + z*z)
+           C(i,j) = (1.0d0 - dexp(-rescale_factor_kappa * dist)) * rescale_factor_kappa_inv
+        end do
+     end do
+
+  case(3)
+
+     do j=1,n
+        do i=1,m
+           x = A(i,1) - B(j,1)
+           y = A(i,2) - B(j,2)
+           z = A(i,3) - B(j,3)
+           dist = dsqrt(x*x + y*y + z*z)
+           C(i,j) = (1.0d0 - dexp(-rescale_factor_kappa * dist)) * rescale_factor_kappa_inv
+        end do
+     end do
+
+  end select
+
+end function qmckl_distance_rescaled_f
+    #+end_src
+
+*** Performance
+
+    This function is more efficient when ~A~ and ~B~ are transposed.
+
+** C interface                                                     :noexport:
+
+   #+CALL: generate_c_interface(table=qmckl_distance_rescaled_args,rettyp=get_value("FRetType"),fname=get_value("Name"))
+
+   #+RESULTS:
+   #+begin_src f90 :tangle (eval f) :comments org :exports none
+   integer(c_int32_t) function qmckl_distance_rescaled &
+       (context, transa, transb, m, n, A, lda, B, ldb, C, ldc, rescale_factor_kappa) &
+       bind(C) result(info)
+
+     use, intrinsic :: iso_c_binding
+     implicit none
+
+     integer (c_int64_t) , intent(in)  , value :: context
+     character           , intent(in)  , value :: transa
+     character           , intent(in)  , value :: transb
+     integer (c_int64_t) , intent(in)  , value :: m
+     integer (c_int64_t) , intent(in)  , value :: n
+     real    (c_double ) , intent(in)          :: A(lda,*)
+     integer (c_int64_t) , intent(in)  , value :: lda
+     real    (c_double ) , intent(in)          :: B(ldb,*)
+     integer (c_int64_t) , intent(in)  , value :: ldb
+     real    (c_double ) , intent(out)         :: C(ldc,n)
+     integer (c_int64_t) , intent(in)  , value :: ldc
+     real    (c_double ) , intent(in)  , value :: rescale_factor_kappa
+
+     integer(c_int32_t), external :: qmckl_distance_rescaled_f
+     info = qmckl_distance_rescaled_f &
+            (context, transa, transb, m, n, A, lda, B, ldb, C, ldc, rescale_factor_kappa)
+
+   end function qmckl_distance_rescaled
+   #+end_src
+
+   #+CALL: generate_f_interface(table=qmckl_distance_rescaled_args,rettyp=get_value("FRetType"),fname=get_value("Name"))
+
+   #+RESULTS:
+   #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
+   interface
+     integer(c_int32_t) function qmckl_distance_rescaled &
+         (context, transa, transb, m, n, A, lda, B, ldb, C, ldc, rescale_factor_kappa) &
+         bind(C)
+       use, intrinsic :: iso_c_binding
+       import
+       implicit none
+
+       integer (c_int64_t) , intent(in)  , value :: context
+       character           , intent(in)  , value :: transa
+       character           , intent(in)  , value :: transb
+       integer (c_int64_t) , intent(in)  , value :: m
+       integer (c_int64_t) , intent(in)  , value :: n
+       real    (c_double ) , intent(in)          :: A(lda,*)
+       integer (c_int64_t) , intent(in)  , value :: lda
+       real    (c_double ) , intent(in)          :: B(ldb,*)
+       integer (c_int64_t) , intent(in)  , value :: ldb
+       real    (c_double ) , intent(out)         :: C(ldc,n)
+       integer (c_int64_t) , intent(in)  , value :: ldc
+       real    (c_double ) , intent(in)  , value :: rescale_factor_kappa
+
+     end function qmckl_distance_rescaled
+   end interface
+   #+end_src
+
+*** Test                                                           :noexport:
 * End of files                                                     :noexport:
 
    #+begin_src c :comments link :tangle (eval c_test)

org/qmckl_electron.org

@@ -63,19 +63,25 @@ int main() {
 
   The following data stored in the context:
 
-   | ~uninitialized~    | int32_t                         | Keeps bit set for uninitialized data                      |
-   | ~num~              | int64_t                         | Total number of electrons                                 |
-   | ~up_num~           | int64_t                         | Number of up-spin electrons                               |
-   | ~down_num~         | int64_t                         | Number of down-spin electrons                             |
-   | ~walk_num~         | int64_t                         | Number of walkers                                         |
-   | ~provided~         | bool                            | If true, ~electron~ is valid                              |
-   | ~coord_new~        | double[walk_num][3][num]        | New set of electron coordinates                           |
-   | ~coord_old~        | double[walk_num][3][num]        | Old set of electron coordinates                           |
-   | ~coord_new_date~   | uint64_t                        | Last modification date of the coordinates                 |
-   | ~ee_distance~      | double[walk_num][num][num]      | Electron-electron distances                               |
-   | ~ee_distance_date~ | uint64_t                        | Last modification date of the electron-electron distances |
-   | ~en_distance~      | double[walk_num][nucl_num][num] | Electron-nucleus distances                                |
-   | ~en_distance_date~ | uint64_t                        | Last modification date of the electron-electron distances |
+   | ~uninitialized~             | int32_t                         | Keeps bit set for uninitialized data                      |
+   | ~num~                       | int64_t                         | Total number of electrons                                 |
+   | ~up_num~                    | int64_t                         | Number of up-spin electrons                               |
+   | ~down_num~                  | int64_t                         | Number of down-spin electrons                             |
+   | ~walk_num~                  | int64_t                         | Number of walkers                                         |
+   | ~rescale_factor_kappa_ee~   | double                          | The distance scaling factor                               |
+   | ~rescale_factor_kappa_en~   | double                          | The distance scaling factor                               |
+   | ~provided~                  | bool                            | If true, ~electron~ is valid                              |
+   | ~coord_new~                 | double[walk_num][3][num]        | New set of electron coordinates                           |
+   | ~coord_old~                 | double[walk_num][3][num]        | Old set of electron coordinates                           |
+   | ~coord_new_date~            | uint64_t                        | Last modification date of the coordinates                 |
+   | ~ee_distance~               | double[walk_num][num][num]      | Electron-electron distances                               |
+   | ~ee_distance_date~          | uint64_t                        | Last modification date of the electron-electron distances |
+   | ~en_distance~               | double[walk_num][nucl_num][num] | Electron-nucleus distances                                |
+   | ~en_distance_date~          | uint64_t                        | Last modification date of the electron-electron distances |
+   | ~ee_distance_rescaled~      | double[walk_num][num][num]      | Electron-electron distances                               |
+   | ~ee_distance_rescaled_date~ | uint64_t                        | Last modification date of the electron-electron distances |
+   | ~en_distance_rescaled~      | double[walk_num][nucl_num][num] | Electron-nucleus distances                                |
+   | ~en_distance_rescaled_date~ | uint64_t                        | Last modification date of the electron-electron distances |
 
 ** Data structure
 
@@ -85,13 +91,19 @@ typedef struct qmckl_electron_struct {
   int64_t   up_num;
   int64_t   down_num;
   int64_t   walk_num;
+  double    rescale_factor_kappa_ee;
+  double    rescale_factor_kappa_en;
   int64_t   coord_new_date;
   int64_t   ee_distance_date;
   int64_t   en_distance_date;
+  int64_t   ee_distance_rescaled_date;
+  int64_t   en_distance_rescaled_date;
   double*   coord_new;
   double*   coord_old;
   double*   ee_distance;
   double*   en_distance;
+  double*   ee_distance_rescaled;
+  double*   en_distance_rescaled;
   int32_t   uninitialized;
   bool      provided;
 } qmckl_electron_struct;
@@ -165,7 +177,7 @@ qmckl_get_electron_num (const qmckl_context context, int64_t* const num) {
   qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
   assert (ctx != NULL);
 
-  int32_t mask = 1;
+  int32_t mask = 1 << 0;
 
   if ( (ctx->electron.uninitialized & mask) != 0) {
     return QMCKL_NOT_PROVIDED;
@@ -193,7 +205,7 @@ qmckl_get_electron_up_num (const qmckl_context context, int64_t* const up_num) {
   qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
   assert (ctx != NULL);
 
-  int32_t mask = 1;
+  int32_t mask = 1 << 0;
 
   if ( (ctx->electron.uninitialized & mask) != 0) {
     return QMCKL_NOT_PROVIDED;
@@ -221,7 +233,7 @@ qmckl_get_electron_down_num (const qmckl_context context, int64_t* const down_nu
   qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
   assert (ctx != NULL);
 
-  int32_t mask = 1;
+  int32_t mask = 1 << 0;
 
   if ( (ctx->electron.uninitialized & mask) != 0) {
     return QMCKL_NOT_PROVIDED;
@@ -260,7 +272,7 @@ qmckl_get_electron_walk_num (const qmckl_context context, int64_t* const walk_nu
   qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
   assert (ctx != NULL);
 
-  int32_t mask = 2;
+  int32_t mask = 1 << 1;
 
   if ( (ctx->electron.uninitialized & mask) != 0) {
     return QMCKL_NOT_PROVIDED;
@@ -272,6 +284,71 @@ qmckl_get_electron_walk_num (const qmckl_context context, int64_t* const walk_nu
 }
       #+end_src
 
+*** Scaling factors Kappa
+
+      #+begin_src c :comments org :tangle (eval h_func) :exports none
+qmckl_exit_code qmckl_get_kappa_ee   (const qmckl_context context, double* const rescale_factor_kappa_ee);
+qmckl_exit_code qmckl_get_kappa_en   (const qmckl_context context, double* const rescale_factor_kappa_en);
+      #+end_src
+
+      #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
+qmckl_exit_code
+qmckl_get_kappa_ee (const qmckl_context context, double* const rescale_factor_kappa_ee) {
+
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return QMCKL_INVALID_CONTEXT;
+  }
+
+  if (rescale_factor_kappa_ee == NULL) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_2,
+                           "qmckl_get_kappa_ee",
+                           "rescale_factor_kappa_ee is a null pointer");
+  }
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  assert (ctx != NULL);
+
+  int32_t mask = 1 << 2;
+
+  if ( (ctx->electron.uninitialized & mask) != 0) {
+    return QMCKL_NOT_PROVIDED;
+  }
+
+  // TODO: assert (ctx->electron.rescale_factor_kappa_ee > (double) 0);
+  ,*rescale_factor_kappa_ee = ctx->electron.rescale_factor_kappa_ee;
+  return QMCKL_SUCCESS;
+}
+
+
+qmckl_exit_code
+qmckl_get_kappa_en (const qmckl_context context, double* const rescale_factor_kappa_en) {
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return QMCKL_INVALID_CONTEXT;
+  }
+
+  if (rescale_factor_kappa_en == NULL) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_2,
+                           "qmckl_get_kappa_en",
+                           "rescale_factor_kappa_en is a null pointer");
+  }
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  assert (ctx != NULL);
+
+  int32_t mask = 1 << 2;
+
+  if ( (ctx->electron.uninitialized & mask) != 0) {
+    return QMCKL_NOT_PROVIDED;
+  }
+
+  // TODO: assert (ctx->electron.rescale_factor_kappa_en > (double) 0);
+  ,*rescale_factor_kappa_en = ctx->electron.rescale_factor_kappa_en;
+  return QMCKL_SUCCESS;
+}
+      #+end_src
+
 *** Electron coordinates
 
     Returns the current electron coordinates. The pointer is assumed
@@ -359,6 +436,7 @@ qmckl_get_electron_coord (const qmckl_context context, const char transp, double
 
    #+begin_src c :comments org :tangle (eval h_func)
 qmckl_exit_code  qmckl_set_electron_num        (qmckl_context context, const int64_t up_num, const int64_t down_num);
+qmckl_exit_code  qmckl_set_kappa               (qmckl_context context, const double rescale_factor_kappa_ee, const double rescale_factor_kappa_en);
 qmckl_exit_code  qmckl_set_electron_walk_num   (qmckl_context context, const int64_t walk_num);
 qmckl_exit_code  qmckl_set_electron_coord      (qmckl_context context, const char transp, const double* coord);
    #+end_src
@@ -437,7 +515,7 @@ qmckl_set_electron_num(qmckl_context context,
                            "down_num <= 0");
   }
 
-  int32_t mask = 1;
+  int32_t mask = 1 << 0;
 
   ctx->electron.up_num = up_num;
   ctx->electron.down_num = down_num;
@@ -447,6 +525,7 @@ qmckl_set_electron_num(qmckl_context context,
 }
    #+end_src
 
+    The following function sets the number of walkers.
 
    #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
 qmckl_exit_code
@@ -461,13 +540,46 @@ qmckl_set_electron_walk_num(qmckl_context context, const int64_t walk_num) {
                            "walk_num <= 0");
   }
 
-  int32_t mask = 2;
+  int32_t mask = 1 << 1;
   ctx->electron.walk_num = walk_num;
 
   <<post2>>
 }
     #+end_src
 
+   Next we set the rescale parameter for the rescaled distance metric.
+
+   #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
+qmckl_exit_code
+qmckl_set_kappa(qmckl_context context,
+                const double rescale_factor_kappa_ee,
+                const double rescale_factor_kappa_en) {
+  <<pre2>>
+
+  // TODO: Check for 0 values
+  //if (rescale_factor_kappa_ee != 0) {
+  //  return qmckl_failwith( context,
+  //                         QMCKL_INVALID_ARG_2,
+  //                         "qmckl_set_kappa",
+  //                         "rescale_factor_kappa_ee == 0");
+  //}
+
+  //if (rescale_factor_kappa_en <= 0) {
+  //  return qmckl_failwith( context,
+  //                         QMCKL_INVALID_ARG_3,
+  //                         "qmckl_set_kappa",
+  //                         "rescale_factor_kappa_en == 0");
+  //}
+
+  int32_t mask = 1 << 2;
+
+  ctx->electron.rescale_factor_kappa_ee = rescale_factor_kappa_ee;
+  ctx->electron.rescale_factor_kappa_en = rescale_factor_kappa_en;
+
+  <<post2>>
+}
+   #+end_src
+
     The following function sets the electron coordinates of all the
     walkers. When this is done, the pointers to the old and new sets
     of coordinates are swapped, and the new coordinates are
@@ -568,11 +680,14 @@ int64_t walk_num      = chbrclf_walk_num;
 int64_t elec_num      = chbrclf_elec_num;
 int64_t elec_up_num   = chbrclf_elec_up_num;
 int64_t elec_dn_num   = chbrclf_elec_dn_num;
+double  rescale_factor_kappa_ee      = 1.0; // TODO Get rescale_factor_kappa_ee from chbrclf
+double  rescale_factor_kappa_en      = 1.0; // TODO Get rescale_factor_kappa_en from chbrclf
 double* elec_coord    = &(chbrclf_elec_coord[0][0][0]);
 
 int64_t  nucl_num      = chbrclf_nucl_num;
 double*  charge        = chbrclf_charge;
 double*  nucl_coord    = &(chbrclf_nucl_coord[0][0]);
+double   nucl_rescale_factor_kappa    = 1.0; // TODO Change get rescale_factor_kappa from chbrclf example
 
 /* --- */
 
@@ -607,6 +722,26 @@ rc = qmckl_get_electron_num (context, &n);
 assert(rc == QMCKL_SUCCESS);
 assert(n == elec_num);
 
+double k_ee;
+double k_en;
+rc = qmckl_get_kappa_ee (context, &k_ee);
+assert(rc == QMCKL_NOT_PROVIDED);
+
+rc = qmckl_get_kappa_en (context, &k_en);
+assert(rc == QMCKL_NOT_PROVIDED);
+
+rc = qmckl_set_kappa (context, rescale_factor_kappa_ee, rescale_factor_kappa_en);
+assert(rc == QMCKL_SUCCESS);
+assert(!qmckl_electron_provided(context));
+
+rc = qmckl_get_kappa_ee (context, &k_ee);
+assert(rc == QMCKL_SUCCESS);
+assert(k_ee == rescale_factor_kappa_ee);
+
+rc = qmckl_get_kappa_en (context, &k_en);
+assert(rc == QMCKL_SUCCESS);
+assert(k_en == rescale_factor_kappa_en);
+
 
 int64_t w;
 rc = qmckl_get_electron_walk_num (context, &w);
@@ -874,6 +1009,242 @@ assert(fabs(ee_distance[elec_num*elec_num+1]-6.5517646321055665) < 1.e-12);
 
      #+end_src
 
+** Electron-electron rescaled distances
+
+*** Get
+
+    #+begin_src c :comments org :tangle (eval h_func) :noweb yes
+qmckl_exit_code qmckl_get_electron_ee_distance_rescaled(qmckl_context context, double* const distance_rescaled);
+    #+end_src
+
+    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
+qmckl_exit_code qmckl_get_electron_ee_distance_rescaled(qmckl_context context, double* const distance_rescaled)
+{
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return QMCKL_NULL_CONTEXT;
+  }
+
+  qmckl_exit_code rc;
+
+  rc = qmckl_provide_ee_distance_rescaled(context);
+  if (rc != QMCKL_SUCCESS) return rc;
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  assert (ctx != NULL);
+
+  size_t sze = ctx->electron.num * ctx->electron.num * ctx->electron.walk_num;
+  memcpy(distance_rescaled, ctx->electron.ee_distance_rescaled, 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_ee_distance_rescaled(qmckl_context context);
+    #+end_src
+
+    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
+qmckl_exit_code qmckl_provide_ee_distance_rescaled(qmckl_context context)
+{
+
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return QMCKL_NULL_CONTEXT;
+  }
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  assert (ctx != NULL);
+
+
+ /* Compute if necessary */
+  if (ctx->electron.coord_new_date > ctx->electron.ee_distance_rescaled_date) {
+
+ /* Allocate array */
+    if (ctx->electron.ee_distance_rescaled == NULL) {
+
+      qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
+      mem_info.size = ctx->electron.num * ctx->electron.num *
+        ctx->electron.walk_num * sizeof(double);
+      double* ee_distance_rescaled = (double*) qmckl_malloc(context, mem_info);
+
+      if (ee_distance_rescaled == NULL) {
+        return qmckl_failwith( context,
+                               QMCKL_ALLOCATION_FAILED,
+                               "qmckl_ee_distance_rescaled",
+                               NULL);
+      }
+      ctx->electron.ee_distance_rescaled = ee_distance_rescaled;
+    }
+
+    qmckl_exit_code rc =
+      qmckl_compute_ee_distance_rescaled(context,
+                                ctx->electron.num,
+                                ctx->electron.rescale_factor_kappa_en,
+                                ctx->electron.walk_num,
+                                ctx->electron.coord_new,
+                                ctx->electron.ee_distance_rescaled);
+    if (rc != QMCKL_SUCCESS) {
+      return rc;
+    }
+
+    ctx->electron.ee_distance_rescaled_date = ctx->date;
+  }
+
+  return QMCKL_SUCCESS;
+}
+    #+end_src
+
+*** Compute
+   :PROPERTIES:
+   :Name:     qmckl_compute_ee_distance_rescaled
+   :CRetType: qmckl_exit_code
+   :FRetType: qmckl_exit_code
+   :END:
+
+    #+NAME: qmckl_ee_distance_rescaled_args
+   | qmckl_context | context                                   | in  | Global state                   |
+   | int64_t       | elec_num                                  | in  | Number of electrons            |
+   | double        | rescale_factor_kappa_ee                   | in  | Factor to rescale ee distances |
+   | int64_t       | walk_num                                  | in  | Number of walkers              |
+   | double        | coord[walk_num][3][elec_num]              | in  | Electron coordinates           |
+   | double        | ee_distance[walk_num][elec_num][elec_num] | out | Electron-electron distances    |
+
+    #+begin_src f90 :comments org :tangle (eval f) :noweb yes
+integer function qmckl_compute_ee_distance_rescaled_f(context, elec_num, rescale_factor_kappa_ee, walk_num, &
+     coord, ee_distance_rescaled) &
+     result(info)
+  use qmckl
+  implicit none
+  integer(qmckl_context), intent(in)  :: context
+  integer*8             , intent(in)  :: elec_num
+  double precision      , intent(in)  :: rescale_factor_kappa_ee
+  integer*8             , intent(in)  :: walk_num
+  double precision      , intent(in)  :: coord(elec_num,3,walk_num)
+  double precision      , intent(out) :: ee_distance_rescaled(elec_num,elec_num,walk_num)
+
+  integer*8 :: k
+
+  info = QMCKL_SUCCESS
+
+  if (context == QMCKL_NULL_CONTEXT) then
+     info = QMCKL_INVALID_CONTEXT
+     return
+  endif
+
+  if (elec_num <= 0) then
+     info = QMCKL_INVALID_ARG_2
+     return
+  endif
+
+  if (walk_num <= 0) then
+     info = QMCKL_INVALID_ARG_3
+     return
+  endif
+
+  do k=1,walk_num
+     info = qmckl_distance_rescaled(context, 'T', 'T', elec_num, elec_num, &
+          coord(1,1,k), elec_num, &
+          coord(1,1,k), elec_num, &
+          ee_distance_rescaled(1,1,k), elec_num, rescale_factor_kappa_ee)
+     if (info /= QMCKL_SUCCESS) then
+        exit
+     endif
+  end do
+
+end function qmckl_compute_ee_distance_rescaled_f
+    #+end_src
+
+    #+begin_src c :tangle (eval h_private_func) :comments org :exports none
+qmckl_exit_code qmckl_compute_ee_distance_rescaled (
+          const qmckl_context context,
+          const int64_t elec_num,
+          const double rescale_factor_kappa_ee,
+          const int64_t walk_num,
+          const double* coord,
+          double* const ee_distance_rescaled );
+    #+end_src
+
+    #+CALL: generate_c_interface(table=qmckl_ee_distance_rescaled_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
+
+    #+RESULTS:
+    #+begin_src f90 :tangle (eval f) :comments org :exports none
+    integer(c_int32_t) function qmckl_compute_ee_distance_rescaled &
+        (context, elec_num, rescale_factor_kappa_ee, walk_num, coord, ee_distance_rescaled) &
+        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 :: elec_num
+      real    (c_double ) , intent(in)  , value :: rescale_factor_kappa_ee
+      integer (c_int64_t) , intent(in)  , value :: walk_num
+      real    (c_double ) , intent(in)          :: coord(elec_num,3,walk_num)
+      real    (c_double ) , intent(out)         :: ee_distance_rescaled(elec_num,elec_num,walk_num)
+
+      integer(c_int32_t), external :: qmckl_compute_ee_distance_rescaled_f
+      info = qmckl_compute_ee_distance_rescaled_f &
+             (context, elec_num, rescale_factor_kappa_ee, walk_num, coord, ee_distance_rescaled)
+
+    end function qmckl_compute_ee_distance_rescaled
+    #+end_src
+
+*** Test
+
+    #+begin_src python :results output :exports none
+import numpy as np
+
+elec_1_w1 = np.array( [ -2.26995253563, -5.15737533569, -2.22940072417 ])
+elec_2_w1 = np.array( [  3.51983380318, -1.08717381954, -1.19617708027 ])
+elec_1_w2 = np.array( [ -2.34410619736, -3.20016115904, -1.53496759012 ])
+elec_2_w2 = np.array( [  3.17996025085, -1.40260577202,  1.49473607540 ])
+
+print ( "[0][0][0] : ", np.linalg.norm(elec_1_w1-elec_1_w1) )
+print ( "[0][1][0] : ", np.linalg.norm(elec_1_w1-elec_2_w1) )
+print ( "[1][0][0] : ", np.linalg.norm(elec_2_w1-elec_1_w1) )
+print ( "[0][0][1] : ", np.linalg.norm(elec_1_w2-elec_1_w2) )
+print ( "[0][1][1] : ", np.linalg.norm(elec_1_w2-elec_2_w2) )
+print ( "[1][0][1] : ", np.linalg.norm(elec_2_w2-elec_1_w2) )
+    #+end_src
+
+    #+RESULTS:
+    : [0][0][0] :  0.0
+    : [0][1][0] :  7.152322512964209
+    : [1][0][0] :  7.152322512964209
+    : [0][0][1] :  0.0
+    : [0][1][1] :  6.5517646321055665
+    : [1][0][1] :  6.5517646321055665
+
+     #+begin_src c :tangle (eval c_test)
+assert(qmckl_electron_provided(context));
+
+
+double ee_distance_rescaled[walk_num * elec_num * elec_num];
+rc = qmckl_get_electron_ee_distance_rescaled(context, ee_distance);
+
+// TODO: Get exact values
+//// (e1,e2,w)
+//// (0,0,0) == 0.
+//assert(ee_distance[0] == 0.);
+//
+//// (1,0,0) == (0,1,0)
+//assert(ee_distance[1] == ee_distance[elec_num]);
+//
+//// value of (1,0,0)
+//assert(fabs(ee_distance[1]-7.152322512964209) < 1.e-12);
+//
+//// (0,0,1) == 0.
+//assert(ee_distance[elec_num*elec_num] == 0.);
+//
+//// (1,0,1) == (0,1,1)
+//assert(ee_distance[elec_num*elec_num+1] == ee_distance[elec_num*elec_num+elec_num]);
+//
+//// value of (1,0,1)
+//assert(fabs(ee_distance[elec_num*elec_num+1]-6.5517646321055665) < 1.e-12);
+
+     #+end_src
+
 ** Electron-nucleus distances
 
 *** Get
@@ -1104,6 +1475,9 @@ assert(qmckl_electron_provided(context));
 rc = qmckl_set_nucleus_num (context, nucl_num);
 assert(rc == QMCKL_SUCCESS);
 
+rc = qmckl_set_nucleus_kappa (context, nucl_rescale_factor_kappa);
+assert(rc == QMCKL_SUCCESS);
+
 rc = qmckl_set_nucleus_charge (context, charge);
 assert (rc == QMCKL_SUCCESS);
 
@@ -1138,6 +1512,285 @@ assert(fabs(en_distance[1][0][1] - 3.1804527583077356) < 1.e-12);
 
      #+end_src
 
+** Electron-nucleus rescaled distances
+
+*** Get
+
+    #+begin_src c :comments org :tangle (eval h_func) :noweb yes
+qmckl_exit_code qmckl_get_electron_en_distance_rescaled(qmckl_context context, double* distance_rescaled);
+    #+end_src
+
+    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
+qmckl_exit_code qmckl_get_electron_en_distance_rescaled(qmckl_context context, double* distance_rescaled)
+{
+
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return QMCKL_NULL_CONTEXT;
+  }
+
+  qmckl_exit_code rc;
+
+  rc = qmckl_provide_en_distance_rescaled(context);
+  if (rc != QMCKL_SUCCESS) return rc;
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  assert (ctx != NULL);
+
+  size_t sze = ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num;
+  memcpy(distance_rescaled, ctx->electron.en_distance_rescaled, 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_en_distance_rescaled(qmckl_context context);
+    #+end_src
+
+    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
+qmckl_exit_code qmckl_provide_en_distance_rescaled(qmckl_context context)
+{
+
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return QMCKL_NULL_CONTEXT;
+  }
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  assert (ctx != NULL);
+
+  if (!(ctx->nucleus.provided)) {
+    return QMCKL_NOT_PROVIDED;
+  }
+
+ /* Compute if necessary */
+  if (ctx->electron.coord_new_date > ctx->electron.en_distance_rescaled_date) {
+
+ /* Allocate array */
+    if (ctx->electron.en_distance_rescaled == NULL) {
+
+      qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
+      mem_info.size = ctx->electron.num * ctx->nucleus.num *
+        ctx->electron.walk_num * sizeof(double);
+      double* en_distance_rescaled = (double*) qmckl_malloc(context, mem_info);
+
+      if (en_distance_rescaled == NULL) {
+        return qmckl_failwith( context,
+                               QMCKL_ALLOCATION_FAILED,
+                               "qmckl_en_distance_rescaled",
+                               NULL);
+      }
+      ctx->electron.en_distance_rescaled = en_distance_rescaled;
+    }
+
+    qmckl_exit_code rc =
+      qmckl_compute_en_distance_rescaled(context,
+                                ctx->electron.num,
+                                ctx->nucleus.num,
+                                ctx->electron.rescale_factor_kappa_en,
+                                ctx->electron.walk_num,
+                                ctx->electron.coord_new,
+                                ctx->nucleus.coord,
+                                ctx->electron.en_distance_rescaled);
+    if (rc != QMCKL_SUCCESS) {
+      return rc;
+    }
+
+    ctx->electron.en_distance_rescaled_date = ctx->date;
+  }
+
+  return QMCKL_SUCCESS;
+}
+    #+end_src
+
+*** Compute
+   :PROPERTIES:
+   :Name:     qmckl_compute_en_distance_rescaled
+   :CRetType: qmckl_exit_code
+   :FRetType: qmckl_exit_code
+   :END:
+
+    #+NAME: qmckl_en_distance_rescaled_args
+   | qmckl_context | context                                                 | in  | Global state                      |
+   | int64_t       | elec_num                                                | in  | Number of electrons               |
+   | int64_t       | nucl_num                                                | in  | Number of nuclei                  |
+   | double        | rescale_factor_kappa_en                          | in  | The factor for rescaled distances |
+   | int64_t       | walk_num                                                | in  | Number of walkers                 |
+   | double        | elec_coord[walk_num][3][elec_num]                       | in  | Electron coordinates              |
+   | double        | nucl_coord[3][elec_num]                                 | in  | Nuclear coordinates               |
+   | double        | en_distance_rescaled_date[walk_num][nucl_num][elec_num] | out | Electron-nucleus distances        |
+
+    #+begin_src f90 :comments org :tangle (eval f) :noweb yes
+integer function qmckl_compute_en_distance_rescaled_f(context, elec_num, nucl_num, rescale_factor_kappa_en, walk_num, elec_coord, &
+     nucl_coord, en_distance_rescaled) &
+     result(info)
+  use qmckl
+  implicit none
+  integer(qmckl_context), intent(in)  :: context
+  integer*8             , intent(in)  :: elec_num
+  integer*8             , intent(in)  :: nucl_num
+  double precision      , intent(in)  :: rescale_factor_kappa_en
+  integer*8             , intent(in)  :: walk_num
+  double precision      , intent(in)  :: elec_coord(elec_num,3,walk_num)
+  double precision      , intent(in)  :: nucl_coord(nucl_num,3)
+  double precision      , intent(out) :: en_distance_rescaled(elec_num,nucl_num,walk_num)
+
+  integer*8 :: k
+
+  info = QMCKL_SUCCESS
+
+  if (context == QMCKL_NULL_CONTEXT) then
+     info = QMCKL_INVALID_CONTEXT
+     return
+  endif
+
+  if (elec_num <= 0) then
+     info = QMCKL_INVALID_ARG_2
+     return
+  endif
+
+  if (nucl_num <= 0) then
+     info = QMCKL_INVALID_ARG_3
+     return
+  endif
+
+  ! TODO: comparison with 0
+  !if (rescale_factor_kappa_en <= 0) then
+  !   info = QMCKL_INVALID_ARG_4
+  !   return
+  !endif
+
+  if (walk_num <= 0) then
+     info = QMCKL_INVALID_ARG_5
+     return
+  endif
+
+  do k=1,walk_num
+     info = qmckl_distance_rescaled(context, 'T', 'T', elec_num, nucl_num, &
+          elec_coord(1,1,k), elec_num, &
+          nucl_coord, nucl_num, &
+          en_distance_rescaled(1,1,k), elec_num, rescale_factor_kappa_en)
+     if (info /= QMCKL_SUCCESS) then
+        exit
+     endif
+  end do
+
+end function qmckl_compute_en_distance_rescaled_f
+    #+end_src
+
+    #+begin_src c :tangle (eval h_private_func) :comments org :exports none
+qmckl_exit_code qmckl_compute_en_distance_rescaled (
+          const qmckl_context context,
+          const int64_t elec_num,
+          const int64_t nucl_num,
+          const double  rescale_factor_kappa_en,
+          const int64_t walk_num,
+          const double* elec_coord,
+          const double* nucl_coord,
+          double* const en_distance_rescaled );
+    #+end_src
+
+    #+CALL: generate_c_interface(table=qmckl_en_distance_rescaled_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
+
+    #+RESULTS:
+    #+begin_src f90 :tangle (eval f) :comments org :exports none
+    integer(c_int32_t) function qmckl_compute_en_distance_rescaled &
+        (context, elec_num, nucl_num, rescale_factor_kappa_en, walk_num, elec_coord, nucl_coord, en_distance_rescaled) &
+        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 :: elec_num
+      integer (c_int64_t) , intent(in)  , value :: nucl_num
+      real    (c_double ) , intent(in)  , value :: rescale_factor_kappa_en
+      integer (c_int64_t) , intent(in)  , value :: walk_num
+      real    (c_double ) , intent(in)          :: elec_coord(elec_num,3,walk_num)
+      real    (c_double ) , intent(in)          :: nucl_coord(elec_num,3)
+      real    (c_double ) , intent(out)         :: en_distance_rescaled(elec_num,nucl_num,walk_num)
+
+      integer(c_int32_t), external :: qmckl_compute_en_distance_rescaled_f
+      info = qmckl_compute_en_distance_rescaled_f &
+             (context, elec_num, nucl_num, rescale_factor_kappa_en, walk_num, elec_coord, nucl_coord, en_distance_rescaled)
+
+    end function qmckl_compute_en_distance_rescaled
+    #+end_src
+
+*** Test
+
+    #+begin_src python :results output :exports none
+import numpy as np
+
+elec_1_w1 = np.array( [ -2.26995253563, -5.15737533569, -2.22940072417 ])
+elec_2_w1 = np.array( [  3.51983380318, -1.08717381954, -1.19617708027 ])
+elec_1_w2 = np.array( [ -2.34410619736, -3.20016115904, -1.53496759012 ])
+elec_2_w2 = np.array( [  3.17996025085, -1.40260577202,  1.49473607540 ])
+nucl_1    = np.array( [ 1.096243353458458e+00, 8.907054016973815e-01, 7.777092280258892e-01 ] )
+nucl_2    = np.array( [ 1.168459237342663e+00, 1.125660720053393e+00, 2.833370314829343e+00 ] )
+
+print ( "[0][0][0] : ", np.linalg.norm(elec_1_w1-nucl_1) )
+print ( "[0][1][0] : ", np.linalg.norm(elec_1_w1-nucl_2) )
+print ( "[0][0][1] : ", np.linalg.norm(elec_2_w1-nucl_1) )
+print ( "[1][0][0] : ", np.linalg.norm(elec_1_w2-nucl_1) )
+print ( "[1][1][0] : ", np.linalg.norm(elec_1_w2-nucl_2) )
+print ( "[1][0][1] : ", np.linalg.norm(elec_2_w2-nucl_1) )
+    #+end_src
+
+    #+RESULTS:
+    : [0][0][0] :  7.546738741619978
+    : [0][1][0] :  8.77102435246984
+    : [0][0][1] :  3.698922010513608
+    : [1][0][0] :  5.824059436060509
+    : [1][1][0] :  7.080482110317645
+    : [1][0][1] :  3.1804527583077356
+
+     #+begin_src c :tangle (eval c_test)
+
+assert(qmckl_electron_provided(context));
+
+rc = qmckl_set_nucleus_num (context, nucl_num);
+assert(rc == QMCKL_SUCCESS);
+
+rc = qmckl_set_nucleus_kappa (context, nucl_rescale_factor_kappa);
+assert(rc == QMCKL_SUCCESS);
+
+rc = qmckl_set_nucleus_charge (context, charge);
+assert (rc == QMCKL_SUCCESS);
+
+rc = qmckl_set_nucleus_coord (context, 'T', nucl_coord);
+assert (rc == QMCKL_SUCCESS);
+
+assert(qmckl_nucleus_provided(context));
+
+double en_distance_rescaled[walk_num][nucl_num][elec_num];
+
+rc = qmckl_get_electron_en_distance_rescaled(context, &(en_distance[0][0][0]));
+assert (rc == QMCKL_SUCCESS);
+
+// TODO: check exact values
+//// (e,n,w) in Fortran notation
+//// (1,1,1)
+//assert(fabs(en_distance[0][0][0] - 7.546738741619978) < 1.e-12);
+//
+//// (1,2,1)
+//assert(fabs(en_distance[0][1][0] - 8.77102435246984) < 1.e-12);
+//
+//// (2,1,1)
+//assert(fabs(en_distance[0][0][1] - 3.698922010513608) < 1.e-12);
+//
+//// (1,1,2)
+//assert(fabs(en_distance[1][0][0] - 5.824059436060509) < 1.e-12);
+//
+//// (1,2,2)
+//assert(fabs(en_distance[1][1][0] - 7.080482110317645) < 1.e-12);
+//
+//// (2,1,2)
+//assert(fabs(en_distance[1][0][1] - 3.1804527583077356) < 1.e-12);
+
+     #+end_src
+
 * End of files                                                     :noexport:
 
   #+begin_src c :tangle (eval h_private_type)

org/qmckl_error.org

@@ -30,12 +30,18 @@
    #+end_src
 
    #+begin_src c :tangle (eval c_test) :noweb yes
+#include <string.h>
+#include <stdio.h>
 #include "qmckl.h"
 #include "assert.h"
+#include "qmckl_error_private_type.h"
 #ifdef HAVE_CONFIG_H
 #include "config.h"
 #endif
 int main() {
+  qmckl_context context;
+  context = qmckl_context_create();
+
    #+end_src
 
    #+end_src
@@ -82,6 +88,16 @@ typedef int32_t qmckl_exit_code;
    | ~QMCKL_INVALID_ARG_8~       |    8 | 'Invalid argument 8'   |
    | ~QMCKL_INVALID_ARG_9~       |    9 | 'Invalid argument 9'   |
    | ~QMCKL_INVALID_ARG_10~      |   10 | 'Invalid argument 10'  |
+   | ~QMCKL_INVALID_ARG_11~      |   11 | 'Invalid argument 11'  |
+   | ~QMCKL_INVALID_ARG_12~      |   12 | 'Invalid argument 12'  |
+   | ~QMCKL_INVALID_ARG_13~      |   13 | 'Invalid argument 13'  |
+   | ~QMCKL_INVALID_ARG_14~      |   14 | 'Invalid argument 14'  |
+   | ~QMCKL_INVALID_ARG_15~      |   15 | 'Invalid argument 15'  |
+   | ~QMCKL_INVALID_ARG_16~      |   16 | 'Invalid argument 16'  |
+   | ~QMCKL_INVALID_ARG_17~      |   17 | 'Invalid argument 17'  |
+   | ~QMCKL_INVALID_ARG_18~      |   18 | 'Invalid argument 18'  |
+   | ~QMCKL_INVALID_ARG_19~      |   19 | 'Invalid argument 19'  |
+   | ~QMCKL_INVALID_ARG_20~      |   20 | 'Invalid argument 20'  |
    | ~QMCKL_FAILURE~             |  101 | 'Failure'              |
    | ~QMCKL_ERRNO~               |  102 | strerror(errno)        |
    | ~QMCKL_INVALID_CONTEXT~     |  103 | 'Invalid context'      |
@@ -130,6 +146,16 @@ return '\n'.join(result)
    #define  QMCKL_INVALID_ARG_8            ((qmckl_exit_code) 8)
    #define  QMCKL_INVALID_ARG_9            ((qmckl_exit_code) 9)
    #define  QMCKL_INVALID_ARG_10           ((qmckl_exit_code) 10)
+   #define  QMCKL_INVALID_ARG_11           ((qmckl_exit_code) 11)
+   #define  QMCKL_INVALID_ARG_12           ((qmckl_exit_code) 12)
+   #define  QMCKL_INVALID_ARG_13           ((qmckl_exit_code) 13)
+   #define  QMCKL_INVALID_ARG_14           ((qmckl_exit_code) 14)
+   #define  QMCKL_INVALID_ARG_15           ((qmckl_exit_code) 15)
+   #define  QMCKL_INVALID_ARG_16           ((qmckl_exit_code) 16)
+   #define  QMCKL_INVALID_ARG_17           ((qmckl_exit_code) 17)
+   #define  QMCKL_INVALID_ARG_18           ((qmckl_exit_code) 18)
+   #define  QMCKL_INVALID_ARG_19           ((qmckl_exit_code) 19)
+   #define  QMCKL_INVALID_ARG_20           ((qmckl_exit_code) 20)
    #define  QMCKL_FAILURE                  ((qmckl_exit_code) 101)
    #define  QMCKL_ERRNO                    ((qmckl_exit_code) 102)
    #define  QMCKL_INVALID_CONTEXT          ((qmckl_exit_code) 103)
@@ -151,6 +177,16 @@ return '\n'.join(result)
       integer(qmckl_exit_code), parameter :: QMCKL_INVALID_ARG_8            = 8
       integer(qmckl_exit_code), parameter :: QMCKL_INVALID_ARG_9            = 9
       integer(qmckl_exit_code), parameter :: QMCKL_INVALID_ARG_10           = 10
+      integer(qmckl_exit_code), parameter :: QMCKL_INVALID_ARG_11           = 11
+      integer(qmckl_exit_code), parameter :: QMCKL_INVALID_ARG_12           = 12
+      integer(qmckl_exit_code), parameter :: QMCKL_INVALID_ARG_13           = 13
+      integer(qmckl_exit_code), parameter :: QMCKL_INVALID_ARG_14           = 14
+      integer(qmckl_exit_code), parameter :: QMCKL_INVALID_ARG_15           = 15
+      integer(qmckl_exit_code), parameter :: QMCKL_INVALID_ARG_16           = 16
+      integer(qmckl_exit_code), parameter :: QMCKL_INVALID_ARG_17           = 17
+      integer(qmckl_exit_code), parameter :: QMCKL_INVALID_ARG_18           = 18
+      integer(qmckl_exit_code), parameter :: QMCKL_INVALID_ARG_19           = 19
+      integer(qmckl_exit_code), parameter :: QMCKL_INVALID_ARG_20           = 20
       integer(qmckl_exit_code), parameter :: QMCKL_FAILURE                  = 101
       integer(qmckl_exit_code), parameter :: QMCKL_ERRNO                    = 102
       integer(qmckl_exit_code), parameter :: QMCKL_INVALID_CONTEXT          = 103
@@ -198,6 +234,94 @@ return '\n'.join(result)
 
    #+end_src
 
+   #+RESULTS: cases
+   #+begin_example
+   case QMCKL_SUCCESS:
+		 return "Success";
+	       break;
+   case QMCKL_INVALID_ARG_1:
+		 return "Invalid argument 1";
+	       break;
+   case QMCKL_INVALID_ARG_2:
+		 return "Invalid argument 2";
+	       break;
+   case QMCKL_INVALID_ARG_3:
+		 return "Invalid argument 3";
+	       break;
+   case QMCKL_INVALID_ARG_4:
+		 return "Invalid argument 4";
+	       break;
+   case QMCKL_INVALID_ARG_5:
+		 return "Invalid argument 5";
+	       break;
+   case QMCKL_INVALID_ARG_6:
+		 return "Invalid argument 6";
+	       break;
+   case QMCKL_INVALID_ARG_7:
+		 return "Invalid argument 7";
+	       break;
+   case QMCKL_INVALID_ARG_8:
+		 return "Invalid argument 8";
+	       break;
+   case QMCKL_INVALID_ARG_9:
+		 return "Invalid argument 9";
+	       break;
+   case QMCKL_INVALID_ARG_10:
+		 return "Invalid argument 10";
+	       break;
+   case QMCKL_INVALID_ARG_11:
+		 return "Invalid argument 11";
+	       break;
+   case QMCKL_INVALID_ARG_12:
+		 return "Invalid argument 12";
+	       break;
+   case QMCKL_INVALID_ARG_13:
+		 return "Invalid argument 13";
+	       break;
+   case QMCKL_INVALID_ARG_14:
+		 return "Invalid argument 14";
+	       break;
+   case QMCKL_INVALID_ARG_15:
+		 return "Invalid argument 15";
+	       break;
+   case QMCKL_INVALID_ARG_16:
+		 return "Invalid argument 16";
+	       break;
+   case QMCKL_INVALID_ARG_17:
+		 return "Invalid argument 17";
+	       break;
+   case QMCKL_INVALID_ARG_18:
+		 return "Invalid argument 18";
+	       break;
+   case QMCKL_INVALID_ARG_19:
+		 return "Invalid argument 19";
+	       break;
+   case QMCKL_INVALID_ARG_20:
+		 return "Invalid argument 20";
+	       break;
+   case QMCKL_FAILURE:
+		 return "Failure";
+	       break;
+   case QMCKL_ERRNO:
+		 return strerror(errno);
+	       break;
+   case QMCKL_INVALID_CONTEXT:
+		 return "Invalid context";
+	       break;
+   case QMCKL_ALLOCATION_FAILED:
+		 return "Allocation failed";
+	       break;
+   case QMCKL_DEALLOCATION_FAILED:
+		 return "De-allocation failed";
+	       break;
+   case QMCKL_NOT_PROVIDED:
+		 return "Not provided";
+	       break;
+   case QMCKL_INVALID_EXIT_CODE:
+		 return "Invalid exit code";
+	       break;
+   #+end_example
+
    # Source
    #+begin_src c :comments org :tangle (eval c) :noweb yes
 const char* qmckl_string_of_error(const qmckl_exit_code error) {
@@ -293,8 +417,54 @@ qmckl_set_error(qmckl_context context,
 }
     #+end_src
 
-* Failing
+* Get the error
 
+  Upon error, the error type and message can be obtained from the
+  context using ~qmckl_get_error~. The message and function name
+  is returned in the variables provided. Therefore, passing a 
+  function name and message is mandatory.
+
+  # Header
+   #+begin_src c :comments org :tangle (eval h_func) :exports none
+qmckl_exit_code
+qmckl_get_error(qmckl_context context,
+                qmckl_exit_code *exit_code,
+                char* function_name,
+                char* message);
+   #+end_src
+
+  # Source
+  #+begin_src c :tangle (eval c)
+qmckl_exit_code
+qmckl_get_error(qmckl_context context,
+                qmckl_exit_code *exit_code,
+                char* function_name,
+                char* message)
+{
+ /* Passing a function name and a message is mandatory. */
+  assert (function_name != NULL);
+  assert (message != NULL);
+
+ /* The context is assumed to exist. */
+  assert (qmckl_context_check(context) != QMCKL_NULL_CONTEXT);
+
+  qmckl_lock(context);
+  {
+    qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+    assert (ctx != NULL); /* Impossible because the context is valid. */
+
+    strncpy(function_name, ctx->error.function, QMCKL_MAX_FUN_LEN-1);
+    strncpy(message      , ctx->error.message , QMCKL_MAX_MSG_LEN-1);
+    (*exit_code) = ctx->error.exit_code;
+  }
+  qmckl_unlock(context);
+
+  return QMCKL_SUCCESS;
+}
+    #+end_src
+
+* Failing
+  
    To make a function fail, the ~qmckl_failwith~ function should be
    called, such that information about the failure is stored in
    the context. The desired exit code is given as an argument, as
@@ -310,7 +480,7 @@ qmckl_exit_code qmckl_failwith(qmckl_context context,
                                const char* function,
                                const char* message) ;
    #+end_src
-
+   
    #+begin_src c :comments org :tangle (eval c)
 qmckl_exit_code qmckl_failwith(qmckl_context context,
                                const qmckl_exit_code exit_code,
@@ -363,7 +533,22 @@ if (x < 0) {
 
 ** Test
    #+begin_src c :comments link :tangle (eval c_test)
-return 0;
+ /* Initialize the variables */
+  char function_name[QMCKL_MAX_FUN_LEN]="";
+  char message[QMCKL_MAX_MSG_LEN]="";
+
+ /* Set the error code to be different from Success */
+  qmckl_exit_code exit_code;
+  exit_code = 1;
+
+  assert (qmckl_set_error(context, exit_code, "qmckl_transpose", "Success") == QMCKL_SUCCESS);
+
+  assert (qmckl_get_error(context, &exit_code, function_name, message) == QMCKL_SUCCESS);
+  assert (exit_code == 1);
+  assert (strcmp(function_name,"qmckl_transpose") == 0);
+  assert (strcmp(message,"Success") == 0);
+
+  return 0;
 }
    #+end_src
 

org/qmckl_nucleus.org

@@ -62,15 +62,18 @@ int main() {
 
   The following data stored in the context:
 
-   | ~uninitialized~    | int32_t          | Keeps bit set for uninitialized data                |
-   | ~num~              | int64_t          | Total number of nuclei                              |
-   | ~provided~         | bool             | If true, ~nucleus~ is valid                         |
-   | ~charge~           | double[num]      | Nuclear charges                                     |
-   | ~coord~            | double[3][num]   | Nuclear coordinates, in transposed format           |
-   | ~nn_distance~      | double[num][num] | Nucleus-nucleus distances                           |
-   | ~nn_distance_date~ | int64_t          | Date when Nucleus-nucleus distances were computed   |
-   | ~repulsion~        | double           | Nuclear repulsion energy                            |
-   | ~repulsion_date~   | int64_t          | Date when the nuclear repulsion energy was computed |
+   | ~uninitialized~             | int32_t          | Keeps bit set for uninitialized data                       |
+   | ~num~                       | int64_t          | Total number of nuclei                                     |
+   | ~provided~                  | bool             | If true, ~nucleus~ is valid                                |
+   | ~charge~                    | double[num]      | Nuclear charges                                            |
+   | ~coord~                     | double[3][num]   | Nuclear coordinates, in transposed format                  |
+   | ~nn_distance~               | double[num][num] | Nucleus-nucleus distances                                  |
+   | ~nn_distance_date~          | int64_t          | Date when Nucleus-nucleus distances were computed          |
+   | ~nn_distance_rescaled~      | double[num][num] | Nucleus-nucleus rescaled distances                         |
+   | ~nn_distance_rescaled_date~ | int64_t          | Date when Nucleus-nucleus rescaled distances were computed |
+   | ~repulsion~                 | double           | Nuclear repulsion energy                                   |
+   | ~repulsion_date~            | int64_t          | Date when the nuclear repulsion energy was computed        |
+   | ~rescale_factor_kappa~      | double           | The distance scaling factor                                |
 
 ** Data structure
 
@@ -79,25 +82,29 @@ typedef struct qmckl_nucleus_struct {
   int64_t   num;
   int64_t   repulsion_date;
   int64_t   nn_distance_date;
+  int64_t   nn_distance_rescaled_date;
   double*   coord;
   double*   charge;
   double*   nn_distance;
+  double*   nn_distance_rescaled;
   double    repulsion;
+  double    rescale_factor_kappa;
   int32_t   uninitialized;
   bool      provided;
 } qmckl_nucleus_struct;
    #+end_src
 
    The ~uninitialized~ integer contains one bit set to one for each
-   initialization function which has not bee called. It becomes equal
+   initialization function which has not been called. It becomes equal
    to zero after all initialization functions have been called. The
    struct is then initialized and ~provided == true~.
 
 ** Access functions
-
+   
    #+begin_src c :comments org :tangle (eval h_func) :exports none
 qmckl_exit_code qmckl_get_nucleus_num         (const qmckl_context context, int64_t* const num);
 qmckl_exit_code qmckl_get_nucleus_charge      (const qmckl_context context, double* const charge);
+qmckl_exit_code qmckl_get_nucleus_kappa       (const qmckl_context context, double* const rescale_factor_kappa);
 qmckl_exit_code qmckl_get_nucleus_coord       (const qmckl_context context, const char transp, double* const coord);
    #+end_src
 
@@ -179,6 +186,31 @@ qmckl_get_nucleus_charge (const qmckl_context context, double* const charge) {
 }
 
 
+qmckl_exit_code
+qmckl_get_nucleus_kappa (const qmckl_context context, double* const rescale_factor_kappa) {
+
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return QMCKL_INVALID_CONTEXT;
+  }
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  assert (ctx != NULL);
+
+  int32_t mask = 1 << 2;
+
+  if ( (ctx->nucleus.uninitialized & mask) != 0) {
+    return qmckl_failwith( context,
+                           QMCKL_NOT_PROVIDED,
+                           "qmckl_get_nucleus_kappa",
+                           "nucleus data is not provided");
+  }
+
+  (*rescale_factor_kappa) = ctx->nucleus.rescale_factor_kappa;
+
+  return QMCKL_SUCCESS;
+}
+
+
 qmckl_exit_code
 qmckl_get_nucleus_coord (const qmckl_context context, const char transp, double* const coord) {
   
@@ -203,7 +235,7 @@ qmckl_get_nucleus_coord (const qmckl_context context, const char transp, double*
   qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
   assert (ctx != NULL);
 
-  int32_t mask = 1 << 2;
+  int32_t mask = 1 << 3;
 
   if ( (ctx->nucleus.uninitialized & mask) != 0) {
     return qmckl_failwith( context,
@@ -263,6 +295,7 @@ bool qmckl_nucleus_provided(const qmckl_context context) {
    #+begin_src c :comments org :tangle (eval h_func)
 qmckl_exit_code  qmckl_set_nucleus_num    (qmckl_context context, const int64_t num);
 qmckl_exit_code  qmckl_set_nucleus_charge (qmckl_context context, const double* charge);
+qmckl_exit_code  qmckl_set_nucleus_kappa  (qmckl_context context, const double rescale_factor_kappa);
 qmckl_exit_code  qmckl_set_nucleus_coord  (qmckl_context context, const char transp, const double* coord);
    #+end_src
 
@@ -347,6 +380,32 @@ qmckl_set_nucleus_charge(qmckl_context context, const double* charge) {
   ctx->nucleus.charge= memcpy(ctx->nucleus.charge, charge, num*sizeof(double));
   assert (ctx->nucleus.charge != NULL);
 
+  <<post2>>
+}
+    #+end_src
+    
+    The following function sets the rescale parameter for the nuclear distances.
+
+   #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
+qmckl_exit_code
+qmckl_set_nucleus_kappa(qmckl_context context, const double rescale_factor_kappa) {
+  <<pre2>>
+
+  //TODO: Check for small values of kappa
+  //if (rescale_factor_kappa == 0) {
+  //  return qmckl_failwith( context,
+  //                         QMCKL_INVALID_ARG_2,
+  //                         "qmckl_set_nucleus_kappa",
+  //                         "rescale_factor_kappa cannot be 0");
+  //}
+
+  int32_t mask = 1 << 2;
+
+  int64_t num;
+  qmckl_exit_code rc;
+
+  ctx->nucleus.rescale_factor_kappa = rescale_factor_kappa;
+
   <<post2>>
 }
     #+end_src
@@ -362,7 +421,7 @@ qmckl_set_nucleus_coord(qmckl_context context, const char transp, const double*
   int64_t nucl_num = (int64_t) 0;
   qmckl_exit_code rc;
 
-  int32_t mask = 1 << 2;
+  int32_t mask = 1 << 3;
 
   rc = qmckl_get_nucleus_num(context, &nucl_num);
   if (rc != QMCKL_SUCCESS) return rc;
@@ -403,6 +462,7 @@ qmckl_set_nucleus_coord(qmckl_context context, const char transp, const double*
 const int64_t   nucl_num      = chbrclf_nucl_num;
 const double*   nucl_charge   = chbrclf_charge;
 const double*   nucl_coord    = &(chbrclf_nucl_coord[0][0]);
+const double    nucl_rescale_factor_kappa    = 1.0; // TODO Change get rescale_factor_kappa from chbrclf example
 
 /* --- */
 
@@ -423,6 +483,19 @@ rc = qmckl_get_nucleus_num (context, &n);
 assert(rc == QMCKL_SUCCESS);
 assert(n == nucl_num);
 
+double k;
+rc = qmckl_get_nucleus_kappa (context, &k);
+assert(rc == QMCKL_NOT_PROVIDED);
+
+
+rc = qmckl_set_nucleus_kappa (context, nucl_rescale_factor_kappa);
+assert(rc == QMCKL_SUCCESS);
+assert(!qmckl_nucleus_provided(context));
+
+rc = qmckl_get_nucleus_kappa (context, &k);
+assert(rc == QMCKL_SUCCESS);
+assert(k == nucl_rescale_factor_kappa);
+
 double nucl_coord2[3*nucl_num];
 
 rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2);
@@ -474,7 +547,7 @@ assert(qmckl_nucleus_provided(context));
   current date is stored.
 
 ** Nucleus-nucleus distances
-
+   
 *** Get
 
     #+begin_src c :comments org :tangle (eval h_func) :noweb yes
@@ -639,6 +712,178 @@ assert(fabs(distance[1]-2.070304721365169) < 1.e-12);
 
      #+end_src
 
+
+** Nucleus-nucleus rescaled distances
+
+*** Get
+
+    #+begin_src c :comments org :tangle (eval h_func) :noweb yes
+qmckl_exit_code qmckl_get_nucleus_nn_distance_rescaled(qmckl_context context, double* distance_rescaled);
+    #+end_src
+
+    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
+qmckl_exit_code qmckl_get_nucleus_nn_distance_rescaled(qmckl_context context, double* distance_rescaled)
+{
+ /* Check input parameters */
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return (char) 0;
+  }
+
+  qmckl_exit_code rc = qmckl_provide_nn_distance_rescaled(context);
+  if (rc != QMCKL_SUCCESS) return rc;
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  assert (ctx != NULL);
+
+  size_t sze = ctx->nucleus.num * ctx->nucleus.num;
+  memcpy(distance_rescaled, ctx->nucleus.nn_distance_rescaled, 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_nn_distance_rescaled(qmckl_context context);
+    #+end_src
+
+    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
+qmckl_exit_code qmckl_provide_nn_distance_rescaled(qmckl_context context)
+{
+ /* Check input parameters */
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return (char) 0;
+  }
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  assert (ctx != NULL);
+
+  if (!ctx->nucleus.provided) return QMCKL_NOT_PROVIDED;
+
+ /* Allocate array */
+  if (ctx->nucleus.nn_distance_rescaled == NULL) {
+
+    qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
+    mem_info.size = ctx->nucleus.num * ctx->nucleus.num * sizeof(double);
+    double* nn_distance_rescaled = (double*) qmckl_malloc(context, mem_info);
+
+    if (nn_distance_rescaled == NULL) {
+      return qmckl_failwith( context,
+                             QMCKL_ALLOCATION_FAILED,
+                             "qmckl_nn_distance_rescaled",
+                             NULL);
+    }
+    ctx->nucleus.nn_distance_rescaled = nn_distance_rescaled;
+  }
+
+  qmckl_exit_code rc =
+    qmckl_compute_nn_distance_rescaled(context,
+                              ctx->nucleus.num,
+                              ctx->nucleus.rescale_factor_kappa,
+                              ctx->nucleus.coord,
+                              ctx->nucleus.nn_distance_rescaled);
+  if (rc != QMCKL_SUCCESS) {
+    return rc;
+  }
+
+  ctx->nucleus.nn_distance_rescaled_date = ctx->date;
+
+  return QMCKL_SUCCESS;
+}
+    #+end_src
+
+*** Compute
+
+    #+NAME: qmckl_nn_distance_rescaled_args
+   | qmckl_context | context                                  | in  | Global state                            |
+   | int64_t       | nucl_num                                 | in  | Number of nuclei                        |
+   | double        | coord[3][nucl_num]                       | in  | Nuclear coordinates (au)                |
+   | double        | nn_distance_rescaled[nucl_num][nucl_num] | out | Nucleus-nucleus rescaled distances (au) |
+
+    #+begin_src f90 :comments org :tangle (eval f) :noweb yes
+integer function qmckl_compute_nn_distance_rescaled_f(context, nucl_num, rescale_factor_kappa, coord, nn_distance_rescaled) &
+     result(info)
+  use qmckl
+  implicit none
+  integer(qmckl_context), intent(in)  :: context
+  integer*8             , intent(in)  :: nucl_num
+  double precision      , intent(in)  :: rescale_factor_kappa
+  double precision      , intent(in)  :: coord(nucl_num,3)
+  double precision      , intent(out) :: nn_distance_rescaled(nucl_num,nucl_num)
+
+  integer*8 :: k
+
+  info = QMCKL_SUCCESS
+
+  if (context == QMCKL_NULL_CONTEXT) then
+     info = QMCKL_INVALID_CONTEXT
+     return
+  endif
+
+  if (nucl_num <= 0) then
+     info = QMCKL_INVALID_ARG_2
+     return
+  endif
+
+  info = qmckl_distance_rescaled(context, 'T', 'T', nucl_num, nucl_num, &
+          coord, nucl_num, &
+          coord, nucl_num, &
+          nn_distance_rescaled, nucl_num, rescale_factor_kappa)
+
+end function qmckl_compute_nn_distance_rescaled_f
+    #+end_src
+
+    #+begin_src c :tangle (eval h_private_func) :comments org :exports none
+qmckl_exit_code qmckl_compute_nn_distance_rescaled (
+          const qmckl_context context,
+          const int64_t nucl_num,
+          const double  rescale_factor_kappa,
+          const double* coord,
+          double* const nn_distance_rescaled );
+    #+end_src
+
+
+    #+CALL: generate_c_interface(table=qmckl_nn_distance_rescaled_args,rettyp="qmckl_exit_code",fname="qmckl_compute_nn_distance")
+
+    #+RESULTS:
+    #+begin_src f90 :tangle (eval f) :comments org :exports none
+    integer(c_int32_t) function qmckl_compute_nn_distance_rescaled &
+        (context, nucl_num, rescale_factor_kappa, coord, nn_distance_rescaled) &
+        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
+      real    (c_double ) , intent(in)  , value :: rescale_factor_kappa
+      real    (c_double ) , intent(in)          :: coord(nucl_num,3)
+      real    (c_double ) , intent(out)         :: nn_distance_rescaled(nucl_num,nucl_num)
+
+      integer(c_int32_t), external :: qmckl_compute_nn_distance_rescaled_f
+      info = qmckl_compute_nn_distance_rescaled_f &
+             (context, nucl_num, rescale_factor_kappa, coord, nn_distance_rescaled)
+
+    end function qmckl_compute_nn_distance_rescaled
+    #+end_src
+
+*** Test
+
+     #+begin_src c :tangle (eval c_test)
+/* Reference input data */
+/* TODO */
+
+//assert(qmckl_nucleus_provided(context));
+//
+//double distance[nucl_num*nucl_num];
+//rc = qmckl_get_nucleus_nn_distance(context, distance);
+//assert(distance[0] == 0.);
+//assert(distance[1] == distance[nucl_num]);
+//assert(fabs(distance[1]-2.070304721365169) < 1.e-12);
+
+     #+end_src
+
 ** Nuclear repulsion energy
 
    \[
@@ -646,7 +891,7 @@ assert(fabs(distance[1]-2.070304721365169) < 1.e-12);
    \]
 
 *** Get
-
+    
     #+begin_src c :comments org :tangle (eval h_func) :noweb yes
 qmckl_exit_code qmckl_get_nucleus_repulsion(qmckl_context context, double* energy);
     #+end_src
@@ -815,7 +1060,7 @@ assert(rep - 318.2309879436158 < 1.e-10);
 }
   #+end_src
 
-**✸ Compute file names
+*** Compute file names
     #+begin_src emacs-lisp
 ; The following is required to compute the file names
 

tools/init.el

@@ -51,6 +51,9 @@
    (makefile . t)
    ))
 
+;; Use python3 instead of python2.7 
+(setq org-babel-python-command "python3")
+
 (add-hook 'org-babel-after-execute-hook 'org-display-inline-images)
 '(indent-tabs-mode nil)