qmckl/org/qmckl_jastrow.org

#+TITLE: Jastrow Factor
#+SETUPFILE: ../tools/theme.setup
#+INCLUDE: ../tools/lib.org

Functions for the calculation of the Jastrow factor \(f_{ee}, f_{en}, f_{een}\).
These are stored in the ~factor_ee~, ~factor_en~, and ~factor_een~ variables.
The ~jastrow~ structure contains all the information required to build
these factors along with their derivatives.

* Headers                                                          :noexport:
  #+begin_src elisp :noexport :results none
(org-babel-lob-ingest "../tools/lib.org")
#+end_src


  #+begin_src c :tangle (eval h_private_type)
#ifndef QMCKL_JASTROW_HPT
#define QMCKL_JASTROW_HPT
#include <stdbool.h>
  #+end_src

  #+begin_src c :tangle (eval c_test) :noweb yes
#include "qmckl.h"
#include <assert.h>
#include <math.h>
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include <stdio.h>
#include "n2.h"

int main() {
  qmckl_context context;
  context = qmckl_context_create();
  #+end_src

  #+begin_src c :tangle (eval c)
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#ifdef HAVE_STDINT_H
#include <stdint.h>
#elif HAVE_INTTYPES_H
#include <inttypes.h>
#endif

#include <stdlib.h>
#include <string.h>
#include <stdbool.h>
#include <assert.h>
#include <math.h>

#include <stdio.h>

#include "qmckl.h"
#include "qmckl_context_private_type.h"
#include "qmckl_memory_private_type.h"
#include "qmckl_memory_private_func.h"
#include "qmckl_jastrow_private_func.h"
#include "qmckl_jastrow_private_type.h"
  #+end_src

* Context
   :PROPERTIES:
   :Name:     qmckl_jastrow
   :CRetType: qmckl_exit_code
   :FRetType: qmckl_exit_code
   :END:

  The following data stored in the context:

   #+NAME: qmckl_jastrow_args
  |-----------+--------------------------------------------+-----+-------------------------------------------------------------------|
  | ~int32_t~ | ~uninitialized~                            | in  | Keeps bit set for uninitialized data                              |
  | ~int64_t~ | ~aord_num~                                 | in  | The number of a coeffecients                                      |
  | ~int64_t~ | ~bord_num~                                 | in  | The number of b coeffecients                                      |
  | ~int64_t~ | ~cord_num~                                 | in  | The number of c coeffecients                                      |
  | ~int64_t~ | ~type_nucl_num~                            | in  | Number of Nucleii types                                           |
  | ~int64_t~ | ~type_nucl_vector[nucl_num]~               | in  | IDs of types of Nucleii                                           |
  | ~double~  | ~aord_vector[aord_num + 1][type_nucl_num]~ | in  | Order of a polynomial coefficients                                |
  | ~double~  | ~bord_vector[bord_num + 1]~                | in  | Order of b polynomial coefficients                                |
  | ~double~  | ~cord_vector[cord_num][type_nucl_num]~     | in  | Order of c polynomial coefficients                                |
  | ~double~  | ~factor_ee[walk_num]~                      | out | Jastrow factor: electron-electron part                            |
  | ~double~  | ~factor_ee_date~                           | out | Jastrow factor: electron-electron part                            |
  | ~double~  | ~factor_en[walk_num]~                      | out | Jastrow factor: electron-nucleus  part                            |
  | ~double~  | ~factor_en_date~                           | out | Jastrow factor: electron-nucleus  part                            |
  | ~double~  | ~factor_een[walk_num]~                     | out | Jastrow factor: electron-electron-nucleus  part                   |
  | ~double~  | ~factor_een_date~                          | out | Jastrow factor: electron-electron-nucleus  part                   |
  | ~double~  | ~factor_ee_deriv_e[4][nelec][walk_num]~    | out | Derivative of the Jastrow factor: electron-electron-nucleus  part |
  | ~double~  | ~factor_ee_deriv_e_date~                   | out | Keep track of the date for the derivative                         |
  | ~double~  | ~factor_en_deriv_e[4][nelec][walk_num]~    | out | Derivative of the Jastrow factor: electron-electron-nucleus  part |
  | ~double~  | ~factor_en_deriv_e_date~                   | out | Keep track of the date for the en derivative                      |
  | ~double~  | ~factor_een_deriv_e[4][nelec][walk_num]~   | out | Derivative of the Jastrow factor: electron-electron-nucleus  part |
  | ~double~  | ~factor_een_deriv_e_date~                  | out | Keep track of the date for the een derivative                     |

  computed data:

  |-----------+-------------------------------------------------------------+-------------------------------------------------|
  | ~int64_t~ | ~dim_cord_vec~                                              | Number of unique C coefficients                 |
  | ~double~  | ~asymp_jasb[2]~                                             | Asymptotic component                            |
  | ~int64_t~ | ~asymp_jasb_date~                                           | Asymptotic component                            |
  | ~double~  | ~coord_vect_full[dim_cord_vec][nucl_num]~                   | vector of non-zero coefficients                 |
  | ~int64_t~ | ~lkpm_of_cindex[4][dim_cord_vec]~                           | Transform l,k,p, and m into consecutive indices |
  | ~double~  | ~tmp_c[elec_num][nucl_num][ncord + 1][ncord][walk_num]~     | vector of non-zero coefficients                 |
  | ~double~  | ~dtmp_c[elec_num][4][nucl_num][ncord + 1][ncord][walk_num]~ | vector of non-zero coefficients                 |

  For H2O we have the following data:

   #+NAME: jastrow_data
   #+BEGIN_SRC python :results output
import numpy as np

elec_num     = 10
nucl_num     = 2
up_num       = 5
down_num     = 5
nucl_coord = np.array([ [0.000000,  0.000000 ],
   [0.000000,  0.000000 ],
   [0.000000,  2.059801 ] ])

elec_coord =    [[[-0.250655104764153      ,  0.503070975550133      ,  -0.166554344502303],
     [-0.587812193472177      , -0.128751981129274      ,   0.187773606533075], 
     [ 1.61335569047166       , -0.615556732874863      ,  -1.43165470979934 ],
     [-4.901239896295210E-003 , -1.120440036458986E-002 ,   1.99761909330422 ],
     [ 0.766647499681200      , -0.293515395797937      ,   3.66454589201239 ],
     [-0.127732483187947      , -0.138975497694196      ,  -8.669850480215846E-002],
     [-0.232271834949124      , -1.059321673434182E-002 ,  -0.504862241464867],
     [ 1.09360863531826       , -2.036103063808752E-003 ,  -2.702796910818986E-002],
     [-0.108090166832043      ,  0.189161729653261      ,   2.15398313919894],
     [ 0.397978144318712      , -0.254277292595981      ,   2.54553335476344]]];

ee_distance_rescaled = [
[ 0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000,
  0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000,
  0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000,
  0.000000000000000E+000],
[ 0.550227800352402      ,0.000000000000000E+000 ,0.000000000000000E+000,
  0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000,
  0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000,
  0.000000000000000E+000],
[ 0.919155060185168      ,0.937695909123175      ,0.000000000000000E+000,
  0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000,
  0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000,
  0.000000000000000E+000],
[ 0.893325429242815      ,0.851181978173561      ,0.978501685226877     ,
  0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000,
  0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000,
  0.000000000000000E+000],
[ 0.982457268305353      ,0.976125002619471      ,0.994349933143149     ,
  0.844077311588328      ,0.000000000000000E+000 ,0.000000000000000E+000,
  0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000,
  0.000000000000000E+000],
[ 0.482407528408731      ,0.414816073699124      ,0.894716035479343     ,
  0.876540187084407      ,0.978921170036895      ,0.000000000000000E+000,
  0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000,
  0.000000000000000E+000],
[ 0.459541909660400      ,0.545007215761510      ,0.883752955884551     ,
  0.918958134888791      ,0.986386936267237      ,0.362209822236419     ,
  0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000,
  0.000000000000000E+000],
[ 0.763732576854455      ,0.817282762358449      ,0.801802919535959     ,
  0.900089095449775      ,0.975704636491453      ,0.707836537586060     ,
  0.755705808346586      ,0.000000000000000E+000 ,0.000000000000000E+000,
  0.000000000000000E+000],
[ 0.904249454052971      ,0.871097965261373      ,0.982717262706270     ,
  0.239901207363622      ,0.836519456769083      ,0.896135326270534     ,
  0.930694340243023      ,0.917708540815567      ,0.000000000000000E+000,
  0.000000000000000E+000],
[ 0.944400908070716      ,0.922589018494961      ,0.984615718580670     ,
  0.514328661540623      ,0.692362267147064      ,0.931894098453677     ,
  0.956034127544344      ,0.931221472309472      ,0.540903688625053     ,
  0.000000000000000E+000]]

en_distance_rescaled = np.transpose(np.array([
[ 0.443570948411811  ,    0.467602196999105    ,  0.893870160799932 ,   
  0.864347190364447  ,    0.976608182392358    ,  0.187563183468210 ,   
  0.426404699872689  ,    0.665107090128166    ,  0.885246991424583 ,   
  0.924902909715270     ],
[ 0.899360150637444     , 0.860035135365386   ,   0.979659405613798 ,   
  6.140678415933776E-002, 0.835118398056681   ,   0.884071658981068 ,   
  0.923860000907362     , 0.905203414522289   ,   0.211286300932359 ,   
  0.492104840907350 ]]))

# symmetrize it
for i in range(elec_num):
  for j in range(elec_num):
    ee_distance_rescaled[i][j] = ee_distance_rescaled[j][i]

type_nucl_num = 1
aord_num     = 5
bord_num     = 5
cord_num     = 23
dim_cord_vec = 23
type_nucl_vector = [ 1, 1]
aord_vector = [ 
[0.000000000000000E+000],
[0.000000000000000E+000],
[-0.380512000000000E+000],
[-0.157996000000000E+000],
[-3.155800000000000E-002],
[2.151200000000000E-002]]

bord_vector = [ 0.500000000000000E-000,  0.153660000000000E-000,  6.722620000000000E-002,
  2.157000000000000E-002,  7.309600000000000E-003,  2.866000000000000E-003]
cord_vector = [ 0.571702000000000E-000, -0.514253000000000E-000, -0.513043000000000E-000, 
  9.486000000000000E-003, -4.205000000000000E-003,  0.426325800000000E-000,
  8.288150000000000E-002,  5.118600000000000E-003, -2.997800000000000E-003,
 -5.270400000000000E-003, -7.499999999999999E-005, -8.301649999999999E-002,
  1.454340000000000E-002,  5.143510000000000E-002,  9.250000000000000E-004,
 -4.099100000000000E-003,  4.327600000000000E-003, -1.654470000000000E-003,
  2.614000000000000E-003, -1.477000000000000E-003, -1.137000000000000E-003,
 -4.010475000000000E-002,  6.106710000000000E-003 ]
cord_vector_full = [
[ 0.571702000000000E-000, -0.514253000000000E-000, -0.513043000000000E-000, 
  9.486000000000000E-003, -4.205000000000000E-003,  0.426325800000000E-000,
  8.288150000000000E-002,  5.118600000000000E-003, -2.997800000000000E-003,
 -5.270400000000000E-003, -7.499999999999999E-005, -8.301649999999999E-002,
  1.454340000000000E-002,  5.143510000000000E-002,  9.250000000000000E-004,
 -4.099100000000000E-003,  4.327600000000000E-003, -1.654470000000000E-003,
  2.614000000000000E-003, -1.477000000000000E-003, -1.137000000000000E-003,
 -4.010475000000000E-002,  6.106710000000000E-003 ],
[ 0.571702000000000E-000, -0.514253000000000E-000, -0.513043000000000E-000, 
  9.486000000000000E-003, -4.205000000000000E-003,  0.426325800000000E-000,
  8.288150000000000E-002,  5.118600000000000E-003, -2.997800000000000E-003,
 -5.270400000000000E-003, -7.499999999999999E-005, -8.301649999999999E-002,
  1.454340000000000E-002,  5.143510000000000E-002,  9.250000000000000E-004,
 -4.099100000000000E-003,  4.327600000000000E-003, -1.654470000000000E-003,
  2.614000000000000E-003, -1.477000000000000E-003, -1.137000000000000E-003,
 -4.010475000000000E-002,  6.106710000000000E-003 ],
]
lkpm_of_cindex = [[1 , 1 , 2 , 0], 
          [0 , 0 , 2 , 1],
          [1 , 2 , 3 , 0],
          [2 , 1 , 3 , 0],
          [0 , 1 , 3 , 1],
          [1 , 0 , 3 , 1],
          [1 , 3 , 4 , 0],
          [2 , 2 , 4 , 0],
          [0 , 2 , 4 , 1],
          [3 , 1 , 4 , 0],
          [1 , 1 , 4 , 1],
          [2 , 0 , 4 , 1],
          [0 , 0 , 4 , 2],
          [1 , 4 , 5 , 0],
          [2 , 3 , 5 , 0],
          [0 , 3 , 5 , 1],
          [3 , 2 , 5 , 0],
          [1 , 2 , 5 , 1],
          [4 , 1 , 5 , 0],
          [2 , 1 , 5 , 1],
          [0 , 1 , 5 , 2],
          [3 , 0 , 5 , 1],
          [1 , 0 , 5 , 2]]

kappa     = 1.0
kappa_inv = 1.0/kappa
   #+END_SRC
   
   #+RESULTS: jastrow_data

** Data structure

   #+begin_src c :comments org :tangle (eval h_private_type)
typedef struct qmckl_jastrow_struct{
  int32_t  uninitialized;
  int64_t  aord_num;
  int64_t  bord_num;
  int64_t  cord_num;
  int64_t  type_nucl_num;
  int64_t  asymp_jasb_date;
  int64_t  tmp_c_date;
  int64_t  dtmp_c_date;
  int64_t  factor_ee_date;
  int64_t  factor_en_date;
  int64_t  factor_een_date;
  int64_t  factor_ee_deriv_e_date;
  int64_t  factor_en_deriv_e_date;
  int64_t  factor_een_deriv_e_date;
  int64_t* type_nucl_vector;
  double * aord_vector;
  double * bord_vector;
  double * cord_vector;
  double * asymp_jasb;
  double * factor_ee;
  double * factor_en;
  double * factor_een;
  double * factor_ee_deriv_e;
  double * factor_en_deriv_e;
  double * factor_een_deriv_e;
  int64_t  dim_cord_vec;
  double * coord_vect_full;
  double * tmp_c;
  double * dtmp_c;
  bool     provided;
  char *   type;
} qmckl_jastrow_struct;
    #+end_src


   The ~uninitialized~ integer contains one bit set to one for each
   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~.
   Some values are initialized by default, and are not concerned by
   this mechanism.

   #+begin_src c :comments org :tangle (eval h_func) 
qmckl_exit_code qmckl_init_jastrow(qmckl_context context);
   #+end_src
   
   #+begin_src c :comments org :tangle (eval c)
qmckl_exit_code qmckl_init_jastrow(qmckl_context context) {

  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return false;
  }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

  ctx->jastrow.uninitialized = (1 << 6) - 1;

 /* Default values */

  return QMCKL_SUCCESS;
}
   #+end_src
   
** Access functions

   #+begin_src c :comments org :tangle (eval h_func) :exports none
qmckl_exit_code  qmckl_get_jastrow_aord_num          (qmckl_context context, int64_t* const aord_num);
qmckl_exit_code  qmckl_get_jastrow_bord_num          (qmckl_context context, int64_t* const bord_num);
qmckl_exit_code  qmckl_get_jastrow_cord_num          (qmckl_context context, int64_t* const bord_num);
qmckl_exit_code  qmckl_get_jastrow_type_nucl_num     (qmckl_context context, int64_t* const type_nucl_num);
qmckl_exit_code  qmckl_get_jastrow_type_nucl_vector  (qmckl_context context, int64_t* const type_nucl_num);
qmckl_exit_code  qmckl_get_jastrow_aord_vector       (qmckl_context context, double * const aord_vector);
qmckl_exit_code  qmckl_get_jastrow_bord_vector       (qmckl_context context, double * const bord_vector);
qmckl_exit_code  qmckl_get_jastrow_cord_vector       (qmckl_context context, double * const cord_vector);
   #+end_src

   Along with these core functions, calculation of the jastrow factor
   requires the following additional information to be set:


   When all the data for the AOs have been provided, the following
   function returns ~true~.

   #+begin_src c :comments org :tangle (eval h_func)
bool      qmckl_jastrow_provided           (const qmckl_context context);
   #+end_src

      #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
bool qmckl_jastrow_provided(const qmckl_context context) {

  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return false;
  }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

  return ctx->jastrow.provided;
}
      #+end_src

   #+NAME:post
   #+begin_src c :exports none
if ( (ctx->jastrow.uninitialized & mask) != 0) {
  return NULL;
}
   #+end_src

   #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_jastrow_aord_num (const qmckl_context context, int64_t* const aord_num) {

  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return (char) 0;
  }
  
  if (aord_num == NULL) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_get_jastrow_aord_num",
                           "aord_num is a null pointer");
  }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

  int32_t mask = 1 << 0;

  if ( (ctx->jastrow.uninitialized & mask) != 0) {
    return QMCKL_NOT_PROVIDED;
  }

  assert (ctx->jastrow.aord_num > 0);
  *aord_num = ctx->jastrow.aord_num;
  return QMCKL_SUCCESS;
}

qmckl_exit_code qmckl_get_jastrow_bord_num (const qmckl_context context, int64_t* const bord_num) {

  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return (char) 0;
  }
  
  if (bord_num == NULL) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_get_jastrow_bord_num",
                           "aord_num is a null pointer");
  }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

  int32_t mask = 1 << 0;

  if ( (ctx->jastrow.uninitialized & mask) != 0) {
    return QMCKL_NOT_PROVIDED;
  }

  assert (ctx->jastrow.bord_num > 0);
  *bord_num = ctx->jastrow.bord_num;
  return QMCKL_SUCCESS;
}

qmckl_exit_code qmckl_get_jastrow_cord_num (const qmckl_context context, int64_t* const cord_num) {

  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return (char) 0;
  }
  
  if (cord_num == NULL) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_get_jastrow_cord_num",
                           "aord_num is a null pointer");
  }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

  int32_t mask = 1 << 0;

  if ( (ctx->jastrow.uninitialized & mask) != 0) {
    return QMCKL_NOT_PROVIDED;
  }

  assert (ctx->jastrow.cord_num > 0);
  *cord_num = ctx->jastrow.cord_num;
  return QMCKL_SUCCESS;
}

qmckl_exit_code qmckl_get_jastrow_type_nucl_num (const qmckl_context context, int64_t* const type_nucl_num) {

  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return (char) 0;
  }
  
  if (type_nucl_num == NULL) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_get_jastrow_type_nucl_num",
                           "type_nucl_num is a null pointer");
  }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

  int32_t mask = 1 << 1;

  if ( (ctx->jastrow.uninitialized & mask) != 0) {
    return QMCKL_NOT_PROVIDED;
  }

  assert (ctx->jastrow.type_nucl_num > 0);
  *type_nucl_num = ctx->jastrow.type_nucl_num;
  return QMCKL_SUCCESS;
}

qmckl_exit_code qmckl_get_jastrow_type_nucl_vector (const qmckl_context context, int64_t * const type_nucl_vector) {

  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return (char) 0;
  }

  if (type_nucl_vector == NULL) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_get_jastrow_type_nucl_vector",
                           "type_nucl_vector is a null pointer");
  }
  
  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

  int32_t mask = 1 << 2;

  if ( (ctx->jastrow.uninitialized & mask) != 0) {
    return QMCKL_NOT_PROVIDED;
  }

  assert (ctx->jastrow.type_nucl_vector != NULL);
  memcpy(type_nucl_vector, ctx->jastrow.type_nucl_vector, ctx->jastrow.type_nucl_num*sizeof(int64_t));
  return QMCKL_SUCCESS;
}

qmckl_exit_code qmckl_get_jastrow_aord_vector (const qmckl_context context, double * const aord_vector) {

  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return (char) 0;
  }

  if (aord_vector == NULL) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_get_jastrow_aord_vector",
                           "aord_vector is a null pointer");
  }
  
  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

  int32_t mask = 1 << 3;

  if ( (ctx->jastrow.uninitialized & mask) != 0) {
    return QMCKL_NOT_PROVIDED;
  }

  assert (ctx->jastrow.aord_vector != NULL);
  memcpy(aord_vector, ctx->jastrow.aord_vector, ctx->jastrow.aord_num*sizeof(double));
  return QMCKL_SUCCESS;
}

qmckl_exit_code qmckl_get_jastrow_bord_vector (const qmckl_context context, double * const bord_vector) {

  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return (char) 0;
  }

  if (bord_vector == NULL) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_get_jastrow_bord_vector",
                           "bord_vector is a null pointer");
  }
  
  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

  int32_t mask = 1 << 4;

  if ( (ctx->jastrow.uninitialized & mask) != 0) {
    return QMCKL_NOT_PROVIDED;
  }

  assert (ctx->jastrow.bord_vector != NULL);
  memcpy(bord_vector, ctx->jastrow.bord_vector, ctx->jastrow.bord_num*sizeof(double));
  return QMCKL_SUCCESS;
}

qmckl_exit_code qmckl_get_jastrow_cord_vector (const qmckl_context context, double * const cord_vector) {

  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return (char) 0;
  }

  if (cord_vector == NULL) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_get_jastrow_cord_vector",
                           "cord_vector is a null pointer");
  }
  
  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

  int32_t mask = 1 << 5;

  if ( (ctx->jastrow.uninitialized & mask) != 0) {
    return QMCKL_NOT_PROVIDED;
  }

  assert (ctx->jastrow.cord_vector != NULL);
  memcpy(cord_vector, ctx->jastrow.cord_vector, ctx->jastrow.cord_num*sizeof(double));
  return QMCKL_SUCCESS;
}

   #+end_src

** Initialization functions

   To prepare for the Jastrow and its derivative, all the following functions need to be
   called.

   #+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code  qmckl_set_jastrow_ord_num           (qmckl_context context, const int64_t aord_num, const int64_t bord_num, const int64_t cord_num);
qmckl_exit_code  qmckl_set_jastrow_type_nucl_num     (qmckl_context context, const int64_t type_nucl_num);
qmckl_exit_code  qmckl_set_jastrow_type_nucl_vector  (qmckl_context context, const int64_t* type_nucl_vector, const int64_t nucl_num);
qmckl_exit_code  qmckl_set_jastrow_aord_vector       (qmckl_context context, const double * aord_vector);
qmckl_exit_code  qmckl_set_jastrow_bord_vector       (qmckl_context context, const double * bord_vector);
qmckl_exit_code  qmckl_set_jastrow_cord_vector       (qmckl_context context, const double * cord_vector);
qmckl_exit_code  qmckl_set_jastrow_dependencies      (qmckl_context context);
   #+end_src

   #+NAME:pre2
   #+begin_src c  :exports none
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
  return QMCKL_NULL_CONTEXT;
 }

qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
   #+end_src

   #+NAME:post2
   #+begin_src c  :exports none
ctx->jastrow.uninitialized &= ~mask;
ctx->jastrow.provided = (ctx->jastrow.uninitialized == 0);
if (ctx->jastrow.provided) {
  //qmckl_exit_code rc_ = qmckl_set_jastrow_dependencies(context);
  //if (rc_ != QMCKL_SUCCESS) return rc_;
 }

return QMCKL_SUCCESS;
   #+end_src


   #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_set_jastrow_ord_num(qmckl_context context, const int64_t aord_num, const int64_t bord_num, const int64_t cord_num) {
<<pre2>>

  if (aord_num <= 0) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_set_jastrow_ord_num",
                           "aord_num <= 0");
  }

  if (bord_num <= 0) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_set_jastrow_ord_num",
                           "bord_num <= 0");
  }

  if (cord_num <= 0) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_set_jastrow_ord_num",
                           "cord_num <= 0");
  }

  int32_t mask = 1 << 0;
  ctx->jastrow.aord_num = aord_num;
  ctx->jastrow.bord_num = bord_num;
  ctx->jastrow.cord_num = cord_num;

  <<post2>>
}

qmckl_exit_code qmckl_set_jastrow_type_nucl_num(qmckl_context context, const int64_t type_nucl_num) {
<<pre2>>

  if (type_nucl_num <= 0) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_set_jastrow_type_nucl_num",
                           "type_nucl_num < 0");
  }

  int32_t mask = 1 << 1;
  ctx->jastrow.type_nucl_num = type_nucl_num;

  <<post2>>
}

qmckl_exit_code qmckl_set_jastrow_type_nucl_vector(qmckl_context context, int64_t const * type_nucl_vector, const int64_t nucl_num) {
<<pre2>>

  int32_t mask = 1 << 2;

  int64_t type_nucl_num;
  qmckl_exit_code rc = qmckl_get_jastrow_type_nucl_num(context, &type_nucl_num);
  if (rc != QMCKL_SUCCESS) return rc;
 
  if (type_nucl_num == 0) {
    return qmckl_failwith( context,
                           QMCKL_FAILURE,
                           "qmckl_set_jastrow_type_nucl_vector",
                           "type_nucl_num is not set");
  }

  if (type_nucl_vector == NULL) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_set_jastrow_type_nucl_vector",
                           "type_nucl_vector = NULL");
  }
  
  if (ctx->jastrow.type_nucl_vector != NULL) {
    qmckl_exit_code rc = qmckl_free(context, ctx->jastrow.type_nucl_vector);
    if (rc != QMCKL_SUCCESS) {
      return qmckl_failwith( context, rc,
                             "qmckl_set_type_nucl_vector",
                             NULL);
    }
  }

  qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
  mem_info.size = nucl_num * sizeof(int64_t);
  int64_t* new_array = (int64_t*) qmckl_malloc(context, mem_info);

  if(new_array == NULL) {
    return qmckl_failwith( context,
                           QMCKL_ALLOCATION_FAILED,
                           "qmckl_set_jastrow_type_nucl_vector",
                           NULL);
  }

  memcpy(new_array, type_nucl_vector, mem_info.size);

  ctx->jastrow.type_nucl_vector = new_array;

  <<post2>>
}

qmckl_exit_code qmckl_set_jastrow_aord_vector(qmckl_context context, double const * aord_vector) {
<<pre2>>

  int32_t mask = 1 << 3;

  int64_t aord_num;
  qmckl_exit_code rc = qmckl_get_jastrow_aord_num(context, &aord_num);
  if (rc != QMCKL_SUCCESS) return rc;
 
  int64_t type_nucl_num;
  rc = qmckl_get_jastrow_type_nucl_num(context, &type_nucl_num);
  if (rc != QMCKL_SUCCESS) return rc;
 
  if (aord_num == 0) {
    return qmckl_failwith( context,
                           QMCKL_FAILURE,
                           "qmckl_set_jastrow_coefficient",
                           "aord_num is not set");
  }

  if (aord_vector == NULL) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_set_jastrow_aord_vector",
                           "aord_vector = NULL");
  }
  
  if (ctx->jastrow.aord_vector != NULL) {
    qmckl_exit_code rc = qmckl_free(context, ctx->jastrow.aord_vector);
    if (rc != QMCKL_SUCCESS) {
      return qmckl_failwith( context, rc,
                             "qmckl_set_ord_vector",
                             NULL);
    }
  }

  qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
  mem_info.size = (aord_num + 1) * type_nucl_num * sizeof(double);
  double* new_array = (double*) qmckl_malloc(context, mem_info);

  if(new_array == NULL) {
    return qmckl_failwith( context,
                           QMCKL_ALLOCATION_FAILED,
                           "qmckl_set_jastrow_coefficient",
                           NULL);
  }

  memcpy(new_array, aord_vector, mem_info.size);

  ctx->jastrow.aord_vector = new_array;

  <<post2>>
}

qmckl_exit_code qmckl_set_jastrow_bord_vector(qmckl_context context, double const * bord_vector) {
<<pre2>>

  int32_t mask = 1 << 4;

  int64_t bord_num;
  qmckl_exit_code rc = qmckl_get_jastrow_bord_num(context, &bord_num);
  if (rc != QMCKL_SUCCESS) return rc;

  if (bord_num == 0) {
    return qmckl_failwith( context,
                           QMCKL_FAILURE,
                           "qmckl_set_jastrow_coefficient",
                           "bord_num is not set");
  }

  if (bord_vector == NULL) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_set_jastrow_bord_vector",
                           "bord_vector = NULL");
  }
  
  if (ctx->jastrow.bord_vector != NULL) {
    qmckl_exit_code rc = qmckl_free(context, ctx->jastrow.bord_vector);
    if (rc != QMCKL_SUCCESS) {
      return qmckl_failwith( context, rc,
                             "qmckl_set_ord_vector",
                             NULL);
    }
  }

  qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
  mem_info.size = (bord_num + 1) * sizeof(double);
  double* new_array = (double*) qmckl_malloc(context, mem_info);

  if(new_array == NULL) {
    return qmckl_failwith( context,
                           QMCKL_ALLOCATION_FAILED,
                           "qmckl_set_jastrow_coefficient",
                           NULL);
  }

  memcpy(new_array, bord_vector, mem_info.size);

  ctx->jastrow.bord_vector = new_array;

  <<post2>>
}

qmckl_exit_code qmckl_set_jastrow_cord_vector(qmckl_context context, double const * cord_vector) {
<<pre2>>

  int32_t mask = 1 << 5;

  int64_t cord_num;
  qmckl_exit_code rc = qmckl_get_jastrow_cord_num(context, &cord_num);
  if (rc != QMCKL_SUCCESS) return rc;
 
  int64_t type_nucl_num;
  rc = qmckl_get_jastrow_type_nucl_num(context, &type_nucl_num);
  if (rc != QMCKL_SUCCESS) return rc;

  if (cord_num == 0) {
    return qmckl_failwith( context,
                           QMCKL_FAILURE,
                           "qmckl_set_jastrow_coefficient",
                           "cord_num is not set");
  }

  if (cord_vector == NULL) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_set_jastrow_cord_vector",
                           "cord_vector = NULL");
  }
  
  if (ctx->jastrow.cord_vector != NULL) {
    qmckl_exit_code rc = qmckl_free(context, ctx->jastrow.cord_vector);
    if (rc != QMCKL_SUCCESS) {
      return qmckl_failwith( context, rc,
                             "qmckl_set_ord_vector",
                             NULL);
    }
  }

  qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
  mem_info.size = cord_num * type_nucl_num * sizeof(double);
  double* new_array = (double*) qmckl_malloc(context, mem_info);

  if(new_array == NULL) {
    return qmckl_failwith( context,
                           QMCKL_ALLOCATION_FAILED,
                           "qmckl_set_jastrow_coefficient",
                           NULL);
  }

  memcpy(new_array, cord_vector, mem_info.size);

  ctx->jastrow.cord_vector = new_array;

  <<post2>>
}

qmckl_exit_code qmckl_set_jastrow_dependencies(qmckl_context context) {
<<pre2>>

 /* Check for electron data */
  if (!(ctx->electron.provided)) {
    return qmckl_failwith( context,
                           QMCKL_NOT_PROVIDED,
                           "qmckl_provide_ee_distance",
                           NULL);
  }

 /* Check for nucleus data */
  if (!(ctx->nucleus.provided)) {
    return qmckl_failwith( context,
                           QMCKL_NOT_PROVIDED,
                           "qmckl_provide_en_distance",
                           NULL);
  }

  int32_t mask = 1 << 6;

  <<post2>>
}

   #+end_src

 When the required information is completely entered, other data structures are
 computed to accelerate the calculations. The intermediates factors
 are precontracted using BLAS LEVEL 3 operations for an optimal FLOP count.

   #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_finalize_jastrow(qmckl_context context);
   #+end_src

   #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_finalize_jastrow(qmckl_context context) {

  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);

 /* ----------------------------------- */
 /* Check for the necessary information */
 /* ----------------------------------- */

  /* Check for the electron data
     1. elec_num
     2. ee_distances_rescaled
  */
  if (!(ctx->electron.provided)) {
    return qmckl_failwith( context,
                           QMCKL_NOT_PROVIDED,
                           "qmckl_electron",
                           NULL);
  }

  /* Check for the nucleus data
     1. nucl_num
     2. en_distances_rescaled
  */
  if (!(ctx->nucleus.provided)) {
    return qmckl_failwith( context,
                           QMCKL_NOT_PROVIDED,
                           "qmckl_nucleus",
                           NULL);
  }

  qmckl_exit_code rc = QMCKL_FAILURE;


 /* ----------------------------------- */
 /* Start calculation of data           */
 /* ----------------------------------- */
  

}
   #+end_src

** Test
   #+begin_src c :tangle (eval c_test)
/* Reference input data */
int64_t walk_num      = n2_walk_num;
int64_t elec_num      = n2_elec_num;
int64_t elec_up_num   = n2_elec_up_num;
int64_t elec_dn_num   = n2_elec_dn_num;
double  rescale_factor_kappa_ee   = 1.0;
double  rescale_factor_kappa_en   = 1.0;
double  nucl_rescale_factor_kappa = 1.0;
double* elec_coord    = &(n2_elec_coord[0][0][0]);

const double*   nucl_charge   = n2_charge;
int64_t  nucl_num      = n2_nucl_num;
double*  charge        = n2_charge;
double*  nucl_coord    = &(n2_nucl_coord[0][0]);

/* Provide Electron data */

qmckl_exit_code rc;

assert(!qmckl_electron_provided(context));

int64_t n;
rc = qmckl_get_electron_num (context, &n);
assert(rc == QMCKL_NOT_PROVIDED);

rc = qmckl_get_electron_up_num (context, &n);
assert(rc == QMCKL_NOT_PROVIDED);

rc = qmckl_get_electron_down_num (context, &n);
assert(rc == QMCKL_NOT_PROVIDED);


rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_electron_provided(context));

rc = qmckl_get_electron_up_num (context, &n);
assert(rc == QMCKL_SUCCESS);
assert(n == elec_up_num);

rc = qmckl_get_electron_down_num (context, &n);
assert(rc == QMCKL_SUCCESS);
assert(n == elec_dn_num);

rc = qmckl_get_electron_num (context, &n);
assert(rc == QMCKL_SUCCESS);
assert(n == elec_num);

double k_ee = 0.;
double k_en = 0.;
rc = qmckl_get_electron_rescale_factor_ee (context, &k_ee);
assert(rc == QMCKL_SUCCESS);
assert(k_ee == 1.0);

rc = qmckl_get_electron_rescale_factor_en (context, &k_en);
assert(rc == QMCKL_SUCCESS);
assert(k_en == 1.0);

rc = qmckl_set_electron_rescale_factor_en(context, rescale_factor_kappa_en);
assert(rc == QMCKL_SUCCESS);

rc = qmckl_set_electron_rescale_factor_ee(context, rescale_factor_kappa_ee);
assert(rc == QMCKL_SUCCESS);

rc = qmckl_get_electron_rescale_factor_ee (context, &k_ee);
assert(rc == QMCKL_SUCCESS);
assert(k_ee == rescale_factor_kappa_ee);

rc = qmckl_get_electron_rescale_factor_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);
assert(rc == QMCKL_NOT_PROVIDED);


rc = qmckl_set_electron_walk_num (context, walk_num);
assert(rc == QMCKL_SUCCESS);

rc = qmckl_get_electron_walk_num (context, &w);
assert(rc == QMCKL_SUCCESS);
assert(w == walk_num);

assert(qmckl_electron_provided(context));

rc = qmckl_set_electron_coord (context, 'N', elec_coord);
assert(rc == QMCKL_SUCCESS);

double elec_coord2[walk_num*3*elec_num];

rc = qmckl_get_electron_coord (context, 'N', elec_coord2);
assert(rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i<3*elec_num ; ++i) {
  assert( elec_coord[i] == elec_coord2[i] );
 }


/* Provide Nucleus data */

assert(!qmckl_nucleus_provided(context));

rc = qmckl_get_nucleus_num (context, &n);
assert(rc == QMCKL_NOT_PROVIDED);


rc = qmckl_set_nucleus_num (context, nucl_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_nucleus_provided(context));

rc = qmckl_get_nucleus_num (context, &n);
assert(rc == QMCKL_SUCCESS);
assert(n == nucl_num);

double k;
rc = qmckl_get_nucleus_rescale_factor (context, &k);
assert(rc == QMCKL_SUCCESS);
assert(k == 1.0);


rc = qmckl_set_nucleus_rescale_factor (context, nucl_rescale_factor_kappa);
assert(rc == QMCKL_SUCCESS);

rc = qmckl_get_nucleus_rescale_factor (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);
assert(rc == QMCKL_NOT_PROVIDED);

rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]));
assert(rc == QMCKL_SUCCESS);

assert(!qmckl_nucleus_provided(context));

rc = qmckl_get_nucleus_coord (context, 'N', nucl_coord2);
assert(rc == QMCKL_SUCCESS);
for (size_t k=0 ; k<3 ; ++k) {
  for (size_t i=0 ; i<nucl_num ; ++i) {
    assert( nucl_coord[nucl_num*k+i] == nucl_coord2[3*i+k] );
  }
}

rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2);
assert(rc == QMCKL_SUCCESS);
for (size_t i=0 ; i<3*nucl_num ; ++i) {
  assert( nucl_coord[i] == nucl_coord2[i] );
}

double nucl_charge2[nucl_num];

rc = qmckl_get_nucleus_charge(context, nucl_charge2);
assert(rc == QMCKL_NOT_PROVIDED);

rc = qmckl_set_nucleus_charge(context, nucl_charge);
assert(rc == QMCKL_SUCCESS);

rc = qmckl_get_nucleus_charge(context, nucl_charge2);
assert(rc == QMCKL_SUCCESS);
for (size_t i=0 ; i<nucl_num ; ++i) {
  assert( nucl_charge[i] == nucl_charge2[i] );
 }
assert(qmckl_nucleus_provided(context));

    #+end_src

* Computation

  The computed data is stored in the context so that it can be reused
  by different kernels. To ensure that the data is valid, for each
  computed data the date of the context is stored when it is computed.
  To know if some data needs to be recomputed, we check if the date of
  the dependencies are more recent than the date of the data to
  compute. If it is the case, then the data is recomputed and the
  current date is stored.

** Asymptotic component for \(f_{ee}\)

  Calculate the asymptotic component ~asymp_jasb~ to be substracted from the final
  electron-electron jastrow factor \(f_{ee}\). The asymptotic componenet is calculated
  via the ~bord_vector~ and the electron-electron rescale factor ~rescale_factor_kappa~.

  \[
  J_{asymp} = \frac{b_1 \kappa^-1}{1 + b_2  \kappa^-1}
  \]
  
*** Get
    #+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_jastrow_asymp_jasb(qmckl_context context, double* const asymp_jasb);
    #+end_src

    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
qmckl_exit_code qmckl_get_jastrow_asymp_jasb(qmckl_context context, double* const asymp_jasb)
{
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return QMCKL_NULL_CONTEXT;
  }

  qmckl_exit_code rc;

  rc = qmckl_provide_asymp_jasb(context);
  if (rc != QMCKL_SUCCESS) return rc;

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

  size_t sze = 2;
  memcpy(asymp_jasb, ctx->jastrow.asymp_jasb, 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_asymp_jasb(qmckl_context context);
    #+end_src
    
    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
qmckl_exit_code qmckl_provide_asymp_jasb(qmckl_context context)
{

  qmckl_exit_code rc;

  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return QMCKL_NULL_CONTEXT;
  }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

 /* Check if ee kappa is provided */
  double rescale_factor_kappa_ee;
  rc = qmckl_get_electron_rescale_factor_ee(context, &rescale_factor_kappa_ee);
  if(rc != QMCKL_SUCCESS) return rc;

 /* Compute if necessary */
  if (ctx->date > ctx->jastrow.asymp_jasb_date) {

 /* Allocate array */
    if (ctx->jastrow.asymp_jasb == NULL) {

      qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
      mem_info.size = 2 * sizeof(double);
      double* asymp_jasb = (double*) qmckl_malloc(context, mem_info);

      if (asymp_jasb == NULL) {
        return qmckl_failwith( context,
                               QMCKL_ALLOCATION_FAILED,
                               "qmckl_asymp_jasb",
                               NULL);
      }
      ctx->jastrow.asymp_jasb = asymp_jasb;
    }

    qmckl_exit_code rc =
      qmckl_compute_asymp_jasb(context,
                                ctx->jastrow.bord_num,
                                ctx->jastrow.bord_vector,
                                rescale_factor_kappa_ee,
                                ctx->jastrow.asymp_jasb);
    if (rc != QMCKL_SUCCESS) {
      return rc;
    }

    ctx->jastrow.asymp_jasb_date = ctx->date;
  }

  return QMCKL_SUCCESS;
}
    #+end_src

*** Compute
   :PROPERTIES:
   :Name:     qmckl_compute_asymp_jasb
   :CRetType: qmckl_exit_code
   :FRetType: qmckl_exit_code
   :END:

    #+NAME: qmckl_asymp_jasb_args
   | qmckl_context | context                   | in  | Global state                |
   | int64_t       | bord_num                  | in  | Number of electrons         |
   | double        | bord_vector[bord_num + 1] | in  | Number of walkers           |
   | double        | rescale_factor_kappa_ee   | in  | Electron coordinates        |
   | double        | asymp_jasb[2]             | out | Electron-electron distances |

    #+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_asymp_jasb_f(context, bord_num, bord_vector, rescale_factor_kappa_ee, asymp_jasb) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer*8             , intent(in)  :: bord_num
  double precision      , intent(in)  :: bord_vector(bord_num)
  double precision      , intent(in)  :: rescale_factor_kappa_ee
  double precision      , intent(out) :: asymp_jasb(2)

  integer*8 :: i, p
  double precision   :: kappa_inv, x, asym_one
  kappa_inv = 1.0d0 / rescale_factor_kappa_ee

  info = QMCKL_SUCCESS

  if (context == QMCKL_NULL_CONTEXT) then
     info = QMCKL_INVALID_CONTEXT
     return
  endif

  if (bord_num <= 0) then
     info = QMCKL_INVALID_ARG_2
     return
  endif

  asym_one = bord_vector(1) * kappa_inv / (1.0d0 + bord_vector(2) * kappa_inv)
  asymp_jasb(:) = (/asym_one, 0.5d0 * asym_one/)
  
  do i = 1, 2
     x = kappa_inv
     do p = 2, bord_num
        x = x * kappa_inv
        asymp_jasb(i) = asymp_jasb(i) + bord_vector(p + 1) * x 
     end do
  end do

end function qmckl_compute_asymp_jasb_f
    #+end_src

    #+CALL: generate_c_header(table=qmckl_asymp_jasb_args,rettyp=get_value("CRetType"),fname=get_value("Name"))

    #+RESULTS:
    #+begin_src c :tangle (eval h_private_func) :comments org
    qmckl_exit_code qmckl_compute_asymp_jasb (
	  const qmckl_context context,
	  const int64_t bord_num,
	  const double* bord_vector,
	  const double rescale_factor_kappa_ee,
	  double* const asymp_jasb ); 
    #+end_src


    #+CALL: generate_c_interface(table=qmckl_asymp_jasb_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_asymp_jasb &
	(context, bord_num, bord_vector, rescale_factor_kappa_ee, asymp_jasb) &
	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 :: bord_num
      real    (c_double ) , intent(in)          :: bord_vector(bord_num + 1)
      real    (c_double ) , intent(in)  , value :: rescale_factor_kappa_ee
      real    (c_double ) , intent(out)         :: asymp_jasb(2)

      integer(c_int32_t), external :: qmckl_compute_asymp_jasb_f
      info = qmckl_compute_asymp_jasb_f &
	     (context, bord_num, bord_vector, rescale_factor_kappa_ee, asymp_jasb)

    end function qmckl_compute_asymp_jasb
    #+end_src

*** Test
    #+name: asymp_jasb
    #+begin_src python :results output :exports none :noweb yes
import numpy as np

<<jastrow_data>>

asym_one = bord_vector[0] * kappa_inv / (1.0 + bord_vector[1]*kappa_inv)
asymp_jasb = np.array([asym_one, 0.5 * asym_one])

for i in range(2):
  x = kappa_inv
  for p in range(1,bord_num):
    x = x * kappa_inv
    asymp_jasb[i] += bord_vector[p + 1] * x

print("asym_one         : ", asym_one)
print("asymp_jasb[0]    : ", asymp_jasb[0])
print("asymp_jasb[1]    : ", asymp_jasb[1])
    #+end_src

    #+RESULTS: asymp_jasb
    : asym_one         :  0.6634291325000664
    : asymp_jasb[0]    :  1.043287918508297
    : asymp_jasb[1]    :  0.7115733522582638
    
    #+RESULTS:
    : asym_one         :  0.43340325572525706
    : asymp_jasb[0]    :  0.5323750557252571
    : asymp_jasb[1]    :  0.31567342786262853

     #+begin_src c :tangle (eval c_test)
assert(qmckl_electron_provided(context));

int64_t type_nucl_num = n2_type_nucl_num;
int64_t* type_nucl_vector = &(n2_type_nucl_vector[0]);
int64_t aord_num = n2_aord_num;
int64_t bord_num = n2_bord_num;
int64_t cord_num = n2_cord_num;
double* aord_vector = &(n2_aord_vector[0][0]);
double* bord_vector = &(n2_bord_vector[0]);
double* cord_vector = &(n2_cord_vector[0][0]);

/* Initialize the Jastrow data */
rc = qmckl_init_jastrow(context);
assert(!qmckl_jastrow_provided(context));

/* Set the data */
rc = qmckl_set_jastrow_ord_num(context, aord_num, bord_num, cord_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_jastrow_type_nucl_num(context, type_nucl_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_jastrow_type_nucl_vector(context, type_nucl_vector, nucl_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_jastrow_aord_vector(context, aord_vector);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_jastrow_bord_vector(context, bord_vector);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_jastrow_cord_vector(context, cord_vector);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_jastrow_dependencies(context);
assert(rc == QMCKL_SUCCESS);

/* Check if Jastrow is properly initialized */
assert(qmckl_jastrow_provided(context));

double asymp_jasb[2];
rc = qmckl_get_jastrow_asymp_jasb(context, asymp_jasb);

// calculate asymp_jasb
assert(fabs(asymp_jasb[0]-0.5323750557252571) < 1.e-12);
assert(fabs(asymp_jasb[1]-0.31567342786262853) < 1.e-12);

     #+end_src


** Electron-electron component \(f_{ee}\)

  Calculate the electron-electron jastrow component ~factor_ee~ using the ~asymp_jasb~
  componenet and the electron-electron rescaled distances ~ee_distance_rescaled~. 

  \[
f_{ee} = \sum_{i,j<i} \left\{ \frac{ \eta B_0 C_{ij}}{1 - B_1 C_{ij}}  - J_{asymp} + \sum^{nord}_{k}B_k C_{ij}^k \right\}
  \]

  
*** Get
    #+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_jastrow_factor_ee(qmckl_context context, double* const factor_ee);
    #+end_src

    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
qmckl_exit_code qmckl_get_jastrow_factor_ee(qmckl_context context, double* const factor_ee)
{
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return QMCKL_NULL_CONTEXT;
  }

  qmckl_exit_code rc;

  rc = qmckl_provide_factor_ee(context);
  if (rc != QMCKL_SUCCESS) return rc;

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

  memcpy(factor_ee, ctx->jastrow.factor_ee, ctx->electron.walk_num*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_factor_ee(qmckl_context context);
    #+end_src
    
    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
qmckl_exit_code qmckl_provide_factor_ee(qmckl_context context)
{

  qmckl_exit_code rc;

  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return QMCKL_NULL_CONTEXT;
  }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

 /* Check if ee rescaled distance is provided */
  rc = qmckl_provide_ee_distance_rescaled(context);
  if(rc != QMCKL_SUCCESS) return rc;

 /* Compute if necessary */
  if (ctx->date > ctx->jastrow.factor_ee_date) {

 /* Allocate array */
    if (ctx->jastrow.factor_ee == NULL) {

      qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
      mem_info.size = ctx->electron.walk_num * sizeof(double);
      double* factor_ee = (double*) qmckl_malloc(context, mem_info);

      if (factor_ee == NULL) {
        return qmckl_failwith( context,
                               QMCKL_ALLOCATION_FAILED,
                               "qmckl_provide_factor_ee",
                               NULL);
      }
      ctx->jastrow.factor_ee = factor_ee;
    }
    
    qmckl_exit_code rc =
      qmckl_compute_factor_ee(context,
                              ctx->electron.walk_num,
                              ctx->electron.num,
                              ctx->electron.up_num,
                              ctx->jastrow.bord_num,
                              ctx->jastrow.bord_vector,
                              ctx->electron.ee_distance_rescaled,
                              ctx->jastrow.asymp_jasb,
                              ctx->jastrow.factor_ee);
    if (rc != QMCKL_SUCCESS) {
      return rc;
    }

    ctx->jastrow.factor_ee_date = ctx->date;
  }

  return QMCKL_SUCCESS;
}
    #+end_src

*** Compute
   :PROPERTIES:
   :Name:     qmckl_compute_factor_ee
   :CRetType: qmckl_exit_code
   :FRetType: qmckl_exit_code
   :END:

    #+NAME: qmckl_factor_ee_args
   | qmckl_context | context                                            | in  | Global state                |
   | int64_t       | walk_num                                           | in  | Number of walkers           |
   | int64_t       | elec_num                                           | in  | Number of electrons         |
   | int64_t       | up_num                                             | in  | Number of alpha electrons   |
   | int64_t       | bord_num                                           | in  | Number of coefficients      |
   | double        | bord_vector[bord_num + 1]                          | in  | List of coefficients        |
   | double        | ee_distance_rescaled[walk_num][elec_num][elec_num] | in  | Electron-electron distances |
   | double        | asymp_jasb[2]                                      | in  | Electron-electron distances |
   | double        | factor_ee[walk_num]                                | out | Electron-electron distances |

    #+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_factor_ee_f(context, walk_num, elec_num, up_num, bord_num,            &
                                           bord_vector, ee_distance_rescaled, asymp_jasb, factor_ee) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer*8             , intent(in)  :: walk_num, elec_num, bord_num, up_num
  double precision      , intent(in)  :: bord_vector(bord_num)
  double precision      , intent(in)  :: ee_distance_rescaled(walk_num, elec_num, elec_num)
  double precision      , intent(in)  :: asymp_jasb(2)
  double precision      , intent(out) :: factor_ee(walk_num)

  integer*8 :: i, j, p, ipar, nw
  double precision   :: pow_ser, x, spin_fact, power_ser

  info = QMCKL_SUCCESS

  if (context == QMCKL_NULL_CONTEXT) then
     info = QMCKL_INVALID_CONTEXT
     return
  endif

  if (walk_num <= 0) then
     info = QMCKL_INVALID_ARG_2
     return
  endif

  if (elec_num <= 0) then
     info = QMCKL_INVALID_ARG_3
     return
  endif

  if (bord_num <= 0) then
     info = QMCKL_INVALID_ARG_4
     return
  endif

  factor_ee = 0.0d0
  
  do nw =1, walk_num
  do j = 1, elec_num
     do i = 1, j - 1
        x = ee_distance_rescaled(nw,i,j)
        power_ser = 0.0d0
        spin_fact = 1.0d0
        ipar = 1

        do p = 2, bord_num
          x = x * ee_distance_rescaled(nw,i,j)
          power_ser = power_ser + bord_vector(p + 1) * x
        end do

        if(j .LE. up_num .OR. i .GT. up_num) then
          spin_fact = 0.5d0
          ipar = 2
        endif
        
        factor_ee(nw) = factor_ee(nw) + spin_fact * bord_vector(1)  * &
                                ee_distance_rescaled(nw,i,j) / &
                                (1.0d0 + bord_vector(2) *   &
                                ee_distance_rescaled(nw,i,j))  &
                               -asymp_jasb(ipar) + power_ser

     end do
  end do
  end do

end function qmckl_compute_factor_ee_f
    #+end_src

    #+CALL: generate_c_header(table=qmckl_factor_ee_args,rettyp=get_value("CRetType"),fname=get_value("Name"))

    #+RESULTS:
    #+begin_src c :tangle (eval h_func) :comments org
    qmckl_exit_code qmckl_compute_factor_ee (
	  const qmckl_context context,
	  const int64_t walk_num,
	  const int64_t elec_num,
	  const int64_t up_num,
	  const int64_t bord_num,
	  const double* bord_vector,
	  const double* ee_distance_rescaled,
	  const double* asymp_jasb,
	  double* const factor_ee ); 
    #+end_src


    #+CALL: generate_c_interface(table=qmckl_factor_ee_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_factor_ee &
	(context, &
		 walk_num, &
		 elec_num, &
		 up_num, &
		 bord_num, &
		 bord_vector, &
		 ee_distance_rescaled, &
		 asymp_jasb, &
		 factor_ee) &
	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 :: walk_num
      integer (c_int64_t) , intent(in)  , value :: elec_num
      integer (c_int64_t) , intent(in)  , value :: up_num
      integer (c_int64_t) , intent(in)  , value :: bord_num
      real    (c_double ) , intent(in)          :: bord_vector(bord_num + 1)
      real    (c_double ) , intent(in)          :: ee_distance_rescaled(elec_num,elec_num,walk_num)
      real    (c_double ) , intent(in)          :: asymp_jasb(2)
      real    (c_double ) , intent(out)         :: factor_ee(walk_num)

      integer(c_int32_t), external :: qmckl_compute_factor_ee_f
      info = qmckl_compute_factor_ee_f &
	     (context, &
		 walk_num, &
		 elec_num, &
		 up_num, &
		 bord_num, &
		 bord_vector, &
		 ee_distance_rescaled, &
		 asymp_jasb, &
		 factor_ee)

    end function qmckl_compute_factor_ee
    #+end_src

*** Test
    #+begin_src python :results output :exports none :noweb yes
import numpy as np

<<jastrow_data>>

<<asymp_jasb>>

factor_ee = 0.0
for i in range(0,elec_num):
  for j in range(0,i):
    x = ee_distance_rescaled[i][j]
    pow_ser = 0.0
    spin_fact = 1.0
    ipar = 0
    
    for p in range(1,bord_num):
      x = x * ee_distance_rescaled[i][j]
      pow_ser = pow_ser + bord_vector[p + 1] * x

    if(i < up_num or j >= up_num):
      spin_fact = 0.5
      ipar = 1
    
    factor_ee = factor_ee + spin_fact * bord_vector[0] * ee_distance_rescaled[i][j] \
                          / (1.0 + bord_vector[1] * ee_distance_rescaled[i][j])     \
                          - asymp_jasb[ipar] + pow_ser
print("factor_ee        :",factor_ee)

    #+end_src

    #+RESULTS:
    : asym_one         :  0.43340325572525706
    : asymp_jasb[0]    :  0.5323750557252571
    : asymp_jasb[1]    :  0.31567342786262853
    : factor_ee        : -4.282760865958113
    

     #+begin_src c :tangle (eval c_test)
/* Check if Jastrow is properly initialized */
assert(qmckl_jastrow_provided(context));

double factor_ee[walk_num];
rc = qmckl_get_jastrow_factor_ee(context, factor_ee);

// calculate factor_ee
assert(fabs(factor_ee[0]+4.282760865958113) < 1.e-12);

     #+end_src
     
** Electron-electron component derivative \(f'_{ee}\)

  Calculate the derivative of the ~factor_ee~ using the ~ee_distance_rescaled~ and
  the electron-electron rescaled distances derivatives ~ee_distance_rescaled_deriv_e~.
  There are four components, the gradient which has 3 components in the \(x, y, z\) 
  directions and the laplacian as the last component.

  TODO: Add equation

  
*** Get
    #+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_jastrow_factor_ee_deriv_e(qmckl_context context, double* const factor_ee_deriv_e);
    #+end_src

    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
qmckl_exit_code qmckl_get_jastrow_factor_ee_deriv_e(qmckl_context context, double* const factor_ee_deriv_e)
{
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return QMCKL_NULL_CONTEXT;
  }

  qmckl_exit_code rc;

  rc = qmckl_provide_factor_ee_deriv_e(context);
  if (rc != QMCKL_SUCCESS) return rc;

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

  int64_t sze = ctx->electron.walk_num * 4 * ctx->electron.num;
  memcpy(factor_ee_deriv_e, ctx->jastrow.factor_ee_deriv_e, 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_factor_ee_deriv_e(qmckl_context context);
    #+end_src
    
    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
qmckl_exit_code qmckl_provide_factor_ee_deriv_e(qmckl_context context)
{

  qmckl_exit_code rc;

  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return QMCKL_NULL_CONTEXT;
  }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

 /* Check if ee rescaled distance is provided */
  rc = qmckl_provide_ee_distance_rescaled(context);
  if(rc != QMCKL_SUCCESS) return rc;

 /* Check if ee rescaled distance deriv e is provided */
  rc = qmckl_provide_ee_distance_rescaled_deriv_e(context);
  if(rc != QMCKL_SUCCESS) return rc;

 /* Compute if necessary */
  if (ctx->date > ctx->jastrow.factor_ee_deriv_e_date) {

 /* Allocate array */
    if (ctx->jastrow.factor_ee_deriv_e == NULL) {

      qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
      mem_info.size = ctx->electron.walk_num * 4 * ctx->electron.num * sizeof(double);
      double* factor_ee_deriv_e = (double*) qmckl_malloc(context, mem_info);

      if (factor_ee_deriv_e == NULL) {
        return qmckl_failwith( context,
                               QMCKL_ALLOCATION_FAILED,
                               "qmckl_provide_factor_ee_deriv_e",
                               NULL);
      }
      ctx->jastrow.factor_ee_deriv_e = factor_ee_deriv_e;
    }
    
    qmckl_exit_code rc =
      qmckl_compute_factor_ee_deriv_e(context,
                                      ctx->electron.walk_num,
                                      ctx->electron.num,
                                      ctx->electron.up_num,
                                      ctx->jastrow.bord_num,
                                      ctx->jastrow.bord_vector,
                                      ctx->electron.ee_distance_rescaled,
                                      ctx->electron.ee_distance_rescaled_deriv_e,
                                      ctx->jastrow.asymp_jasb,
                                      ctx->jastrow.factor_ee_deriv_e);
    if (rc != QMCKL_SUCCESS) {
      return rc;
    }

    ctx->jastrow.factor_ee_date = ctx->date;
  }

  return QMCKL_SUCCESS;
}
    #+end_src

*** Compute
   :PROPERTIES:
   :Name:     qmckl_compute_factor_ee_deriv_e
   :CRetType: qmckl_exit_code
   :FRetType: qmckl_exit_code
   :END:

    #+NAME: qmckl_factor_ee_deriv_e_args
   | qmckl_context | context                                                       | in  | Global state                |
   | int64_t       | walk_num                                                      | in  | Number of walkers           |
   | int64_t       | elec_num                                                      | in  | Number of electrons         |
   | int64_t       | up_num                                                        | in  | Number of alpha electrons   |
   | int64_t       | bord_num                                                      | in  | Number of coefficients      |
   | double        | bord_vector[bord_num + 1]                                     | in  | List of coefficients        |
   | double        | ee_distance_rescaled[walk_num][elec_num][elec_num]            | in  | Electron-electron distances |
   | double        | ee_distance_rescaled_deriv_e[walk_num][4][elec_num][elec_num] | in  | Electron-electron distances |
   | double        | asymp_jasb[2]                                                 | in  | Electron-electron distances |
   | double        | factor_ee_deriv_e[walk_num][4][elec_num]                      | out | Electron-electron distances |

    #+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_factor_ee_deriv_e_f(context, walk_num, elec_num, up_num, bord_num,            &
                                           bord_vector, ee_distance_rescaled, ee_distance_rescaled_deriv_e,  &
					   asymp_jasb, factor_ee_deriv_e) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer*8             , intent(in)  :: walk_num, elec_num, bord_num, up_num
  double precision      , intent(in)  :: bord_vector(bord_num)
  double precision      , intent(in)  :: ee_distance_rescaled(walk_num, elec_num, elec_num)
  double precision      , intent(in)  :: ee_distance_rescaled_deriv_e(walk_num, 4, elec_num, elec_num)
  double precision      , intent(in)  :: asymp_jasb(2)
  double precision      , intent(out) :: factor_ee_deriv_e(elec_num,4,walk_num)

  integer*8 :: i, j, p, ipar, nw, ii
  double precision   :: x, spin_fact, y
  double precision   :: den, invden, invden2, invden3, xinv
  double precision   :: lap1, lap2, lap3, third
  double precision, dimension(3) :: pow_ser_g
  double precision, dimension(4) :: dx

  info = QMCKL_SUCCESS

  if (context == QMCKL_NULL_CONTEXT) then
     info = QMCKL_INVALID_CONTEXT
     return
  endif

  if (walk_num <= 0) then
     info = QMCKL_INVALID_ARG_2
     return
  endif

  if (elec_num <= 0) then
     info = QMCKL_INVALID_ARG_3
     return
  endif

  if (bord_num <= 0) then
     info = QMCKL_INVALID_ARG_4
     return
  endif

  factor_ee_deriv_e = 0.0d0
  third = 1.0d0 / 3.0d0
  
  do nw =1, walk_num
  do j = 1, elec_num
     do i = 1, elec_num
        x = ee_distance_rescaled(nw, i, j)
        if(abs(x) < 1.0d-18) cycle
        pow_ser_g   = 0.0d0
        spin_fact   = 1.0d0
        den         = 1.0d0 + bord_vector(2) * x
        invden      = 1.0d0 / den
        invden2     = invden * invden
        invden3     = invden2 * invden
        xinv        = 1.0d0 / (x + 1.0d-18)
        ipar        = 1

        do ii = 1, 4
          dx(ii) = ee_distance_rescaled_deriv_e(nw, ii, i, j)
        end do

        if((i .LE. up_num .AND. j .LE. up_num ) .OR.  &
           (i .GT. up_num .AND. j .GT. up_num)) then
          spin_fact = 0.5d0
        endif
        
        lap1 = 0.0d0
        lap2 = 0.0d0
        lap3 = 0.0d0
        do ii = 1, 3
          x = ee_distance_rescaled(nw, i, j)
          if(abs(x) < 1.0d-18) cycle
          do p = 2, bord_num
            y = p * bord_vector(p + 1) * x
            pow_ser_g(ii) = pow_ser_g(ii) + y * dx(ii)
            lap1 = lap1 + (p - 1) * y * xinv * dx(ii) * dx(ii)
            lap2 = lap2 + y
            x = x * ee_distance_rescaled(nw, i, j)
          end do
          
          lap3 = lap3 - 2.0d0 * bord_vector(2) * dx(ii) * dx(ii)

          factor_ee_deriv_e( j, ii, nw) = factor_ee_deriv_e( j, ii, nw) + spin_fact * bord_vector(1)  * &
                                dx(ii) * invden2 + pow_ser_g(ii)
        end do

        ii = 4
        lap2 = lap2 * dx(ii) * third
        lap3 = lap3 + den * dx(ii)
        lap3 = lap3 * (spin_fact * bord_vector(1) * invden3)
        factor_ee_deriv_e( j, ii, nw) = factor_ee_deriv_e( j, ii, nw) + lap1 + lap2 + lap3

     end do
  end do
  end do

end function qmckl_compute_factor_ee_deriv_e_f
    #+end_src

    #+CALL: generate_c_header(table=qmckl_factor_ee_deriv_e_args,rettyp=get_value("CRetType"),fname=get_value("Name"))

    #+RESULTS:
    #+begin_src c :tangle (eval h_func) :comments org
    qmckl_exit_code qmckl_compute_factor_ee_deriv_e (
	  const qmckl_context context,
	  const int64_t walk_num,
	  const int64_t elec_num,
	  const int64_t up_num,
	  const int64_t bord_num,
	  const double* bord_vector,
	  const double* ee_distance_rescaled,
	  const double* ee_distance_rescaled_deriv_e,
	  const double* asymp_jasb,
	  double* const factor_ee_deriv_e ); 
    #+end_src


    #+CALL: generate_c_interface(table=qmckl_factor_ee_deriv_e_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_factor_ee_deriv_e &
	(context, &
		 walk_num, &
		 elec_num, &
		 up_num, &
		 bord_num, &
		 bord_vector, &
		 ee_distance_rescaled, &
		 ee_distance_rescaled_deriv_e, &
		 asymp_jasb, &
		 factor_ee_deriv_e) &
	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 :: walk_num
      integer (c_int64_t) , intent(in)  , value :: elec_num
      integer (c_int64_t) , intent(in)  , value :: up_num
      integer (c_int64_t) , intent(in)  , value :: bord_num
      real    (c_double ) , intent(in)          :: bord_vector(bord_num + 1)
      real    (c_double ) , intent(in)          :: ee_distance_rescaled(elec_num,elec_num,walk_num)
      real    (c_double ) , intent(in)          :: ee_distance_rescaled_deriv_e(elec_num,elec_num,4,walk_num)
      real    (c_double ) , intent(in)          :: asymp_jasb(2)
      real    (c_double ) , intent(out)         :: factor_ee_deriv_e(elec_num,4,walk_num)

      integer(c_int32_t), external :: qmckl_compute_factor_ee_deriv_e_f
      info = qmckl_compute_factor_ee_deriv_e_f &
	     (context, &
		 walk_num, &
		 elec_num, &
		 up_num, &
		 bord_num, &
		 bord_vector, &
		 ee_distance_rescaled, &
		 ee_distance_rescaled_deriv_e, &
		 asymp_jasb, &
		 factor_ee_deriv_e)

    end function qmckl_compute_factor_ee_deriv_e
    #+end_src

*** Test
    #+begin_src python :results output :exports none :noweb yes
import numpy as np

<<jastrow_data>>

<<asymp_jasb>>

kappa = 1.0

elec_coord = np.array(elec_coord)[0]
elec_dist = np.zeros(shape=(elec_num, elec_num),dtype=float)
for i in range(elec_num):
  for j in range(elec_num):
    elec_dist[i, j] = np.linalg.norm(elec_coord[i] - elec_coord[j])

elec_dist_deriv_e = np.zeros(shape=(4,elec_num, elec_num),dtype=float)
for j in range(elec_num):
  for i in range(elec_num):
    rij_inv = 1.0 / elec_dist[i, j]
    for ii in range(3):
      elec_dist_deriv_e[ii, i, j] = (elec_coord[j][ii] - elec_coord[i][ii]) * rij_inv
    elec_dist_deriv_e[3, i, j] = 2.0 * rij_inv
  elec_dist_deriv_e[:, j, j] = 0.0

ee_distance_rescaled_deriv_e = np.zeros(shape=(4,elec_num,elec_num),dtype=float)
for j in range(elec_num):
  for i in range(elec_num):
    f = 1.0 - kappa * ee_distance_rescaled[i][j]
    for ii in range(4):
      ee_distance_rescaled_deriv_e[ii][i][j] = elec_dist_deriv_e[ii][i][j]
    ee_distance_rescaled_deriv_e[3][i][j] = ee_distance_rescaled_deriv_e[3][i][j] + \
                              (-kappa * ee_distance_rescaled_deriv_e[0][i][j] * ee_distance_rescaled_deriv_e[0][i][j]) + \
                              (-kappa * ee_distance_rescaled_deriv_e[1][i][j] * ee_distance_rescaled_deriv_e[1][i][j]) + \
                              (-kappa * ee_distance_rescaled_deriv_e[2][i][j] * ee_distance_rescaled_deriv_e[2][i][j])
    for ii in range(4):
      ee_distance_rescaled_deriv_e[ii][i][j] = ee_distance_rescaled_deriv_e[ii][i][j] * f

third = 1.0 / 3.0
factor_ee_deriv_e = np.zeros(shape=(4,elec_num),dtype=float)
dx = np.zeros(shape=(4),dtype=float)
pow_ser_g = np.zeros(shape=(4),dtype=float)
for j in range(elec_num):
   for i in range(elec_num):
      x = ee_distance_rescaled[j][i]
      if abs(x) < 1e-18:
        continue
      pow_ser_g = np.zeros(shape=(4),dtype=float)
      spin_fact   = 1.0
      den         = 1.0 + bord_vector[1] * ee_distance_rescaled[j][i]
      invden      = 1.0 / den
      invden2     = invden * invden
      invden3     = invden2 * invden
      xinv        = 1.0 / (ee_distance_rescaled[j][i] + 1.0E-18)
      ipar        = 1

      for ii in range(4):
        dx[ii] = ee_distance_rescaled_deriv_e[ii][j][i]

      if((i <= (up_num-1) and j <= (up_num-1) ) or  \
         (i >  (up_num-1) and j >  (up_num-1))):
        spin_fact = 0.5
      
      lap1 = 0.0
      lap2 = 0.0
      lap3 = 0.0
      for ii in range(3):
        x = ee_distance_rescaled[j][i]
        if x < 1e-18:
          continue
        for p in range(2,bord_num+1):
          y = p * bord_vector[(p-1) + 1] * x
          pow_ser_g[ii] = pow_ser_g[ii] + y * dx[ii]
          lap1 = lap1 + (p - 1) * y * xinv * dx[ii] * dx[ii]
          lap2 = lap2 + y
          x = x * ee_distance_rescaled[j][i]
        
        lap3 = lap3 - 2.0 * bord_vector[1] * dx[ii] * dx[ii]

        factor_ee_deriv_e[ii][j] = factor_ee_deriv_e[ii][j] + spin_fact * bord_vector[0]  * \
                              dx[ii] * invden2 + pow_ser_g[ii]

      ii = 3
      lap2 = lap2 * dx[ii] * third
      lap3 = lap3 + den * dx[ii]
      lap3 = lap3 * (spin_fact * bord_vector[0] * invden3)
      factor_ee_deriv_e[ii][j] = factor_ee_deriv_e[ii][j] + lap1 + lap2 + lap3

print("factor_ee_deriv_e[0][0]:",factor_ee_deriv_e[0][0])
print("factor_ee_deriv_e[1][0]:",factor_ee_deriv_e[1][0])
print("factor_ee_deriv_e[2][0]:",factor_ee_deriv_e[2][0])
print("factor_ee_deriv_e[3][0]:",factor_ee_deriv_e[3][0])
print(factor_ee_deriv_e)

    #+end_src

    #+RESULTS:
    #+begin_example
    asym_one         :  0.43340325572525706
    asymp_jasb[0]    :  0.5323750557252571
    asymp_jasb[1]    :  0.31567342786262853
    factor_ee_deriv_e[0][0]: 0.16364894652107934
    factor_ee_deriv_e[1][0]: -0.6927548119830084
    factor_ee_deriv_e[2][0]: 0.073267755223968
    factor_ee_deriv_e[3][0]: 1.5111672803213185
    [[ 0.16364895  0.60354957 -0.19825547  0.02359797 -0.13123153 -0.18789233
       0.07762515 -0.42459184  0.27920265 -0.2056531 ]
     [-0.69275481  0.15690393  0.09831069  0.18490587  0.04361723  0.3250686
       0.12657961 -0.01736522 -0.40149005  0.17622416]
     [ 0.07326776 -0.27532276  0.22396943  0.18771633 -0.34506246  0.07298062
       0.63302352 -0.00910198 -0.30238713 -0.25908332]
     [ 1.51116728  1.5033247   0.00325003  2.89377255  0.1338393   2.15893795
       1.74732003  0.23561147  2.67455607  0.82810434]]
    #+end_example
    

     #+begin_src c :tangle (eval c_test)
/* Check if Jastrow is properly initialized */
assert(qmckl_jastrow_provided(context));

// calculate factor_ee_deriv_e
double factor_ee_deriv_e[walk_num][4][elec_num];
rc = qmckl_get_jastrow_factor_ee_deriv_e(context, &(factor_ee_deriv_e[0][0][0]));

// check factor_ee_deriv_e
assert(fabs(factor_ee_deriv_e[0][0][0]-0.16364894652107934) < 1.e-12);
assert(fabs(factor_ee_deriv_e[0][1][0]+0.6927548119830084 ) < 1.e-12);
assert(fabs(factor_ee_deriv_e[0][2][0]-0.073267755223968  ) < 1.e-12);
assert(fabs(factor_ee_deriv_e[0][3][0]-1.5111672803213185 ) < 1.e-12);

     #+end_src

** Electron-nucleus component \(f_{en}\)

  Calculate the electron-electron jastrow component ~factor_en~ using the ~aord_vector~
  coeffecients and the electron-nucleus rescaled distances ~en_distance_rescaled~. 

  \[
f_{en} = \sum_{i,j<i} \left\{ \frac{ A_0 C_{ij}}{1 - A_1 C_{ij}} + \sum^{nord}_{k}A_k C_{ij}^k \right\}
  \]

  
*** Get
    #+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_jastrow_factor_en(qmckl_context context, double* const factor_en);
    #+end_src

    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
qmckl_exit_code qmckl_get_jastrow_factor_en(qmckl_context context, double* const factor_en)
{
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return QMCKL_NULL_CONTEXT;
  }

  qmckl_exit_code rc;

  rc = qmckl_provide_factor_en(context);
  if (rc != QMCKL_SUCCESS) return rc;

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

  memcpy(factor_en, ctx->jastrow.factor_en, ctx->electron.walk_num*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_factor_en(qmckl_context context);
    #+end_src
    
    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
qmckl_exit_code qmckl_provide_factor_en(qmckl_context context)
{

  qmckl_exit_code rc;

  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return QMCKL_NULL_CONTEXT;
  }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

 /* Check if en rescaled distance is provided */
  rc = qmckl_provide_en_distance_rescaled(context);
  if(rc != QMCKL_SUCCESS) return rc;

 /* Compute if necessary */
  if (ctx->date > ctx->jastrow.factor_en_date) {

 /* Allocate array */
    if (ctx->jastrow.factor_en == NULL) {

      qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
      mem_info.size = ctx->electron.walk_num * sizeof(double);
      double* factor_en = (double*) qmckl_malloc(context, mem_info);

      if (factor_en == NULL) {
        return qmckl_failwith( context,
                               QMCKL_ALLOCATION_FAILED,
                               "qmckl_provide_factor_en",
                               NULL);
      }
      ctx->jastrow.factor_en = factor_en;
    }
    
    qmckl_exit_code rc =
      qmckl_compute_factor_en(context,
                              ctx->electron.walk_num,
                              ctx->electron.num,
                              ctx->nucleus.num,
                              ctx->jastrow.type_nucl_num,
                              ctx->jastrow.type_nucl_vector,
                              ctx->jastrow.aord_num,
                              ctx->jastrow.aord_vector,
                              ctx->electron.en_distance_rescaled,
                              ctx->jastrow.factor_en);
    if (rc != QMCKL_SUCCESS) {
      return rc;
    }

    ctx->jastrow.factor_en_date = ctx->date;
  }

  return QMCKL_SUCCESS;
}
    #+end_src

*** Compute
   :PROPERTIES:
   :Name:     qmckl_compute_factor_en
   :CRetType: qmckl_exit_code
   :FRetType: qmckl_exit_code
   :END:

    #+NAME: qmckl_factor_en_args
   | qmckl_context | context                                            | in  | Global state               |
   | int64_t       | walk_num                                           | in  | Number of walkers          |
   | int64_t       | elec_num                                           | in  | Number of electrons        |
   | int64_t       | nucl_num                                           | in  | Number of nucleii          |
   | int64_t       | type_nucl_num                                      | in  | Number of unique nuclei    |
   | int64_t       | type_nucl_vector[type_nucl_num]                    | in  | IDs of unique nucleii      |
   | int64_t       | aord_num                                           | in  | Number of coefficients     |
   | double        | aord_vector[aord_num + 1][type_nucl_num]           | in  | List of coefficients       |
   | double        | en_distance_rescaled[walk_num][nucl_num][elec_num] | in  | Electron-nucleus distances |
   | double        | factor_en[walk_num]                                | out | Electron-nucleus jastrow   |

    #+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_factor_en_f(context, walk_num, elec_num, nucl_num, type_nucl_num,            &
                                           type_nucl_vector, aord_num, aord_vector,                         &
                                           en_distance_rescaled, factor_en) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer*8             , intent(in)  :: walk_num, elec_num, aord_num, nucl_num, type_nucl_num
  integer*8             , intent(in)  :: type_nucl_vector(type_nucl_num)
  double precision      , intent(in)  :: aord_vector(aord_num, nucl_num)
  double precision      , intent(in)  :: en_distance_rescaled(walk_num, elec_num, nucl_num)
  double precision      , intent(out) :: factor_en(walk_num)

  integer*8 :: i, a, p, ipar, nw
  double precision   :: x, spin_fact, power_ser

  info = QMCKL_SUCCESS

  if (context == QMCKL_NULL_CONTEXT) then
     info = QMCKL_INVALID_CONTEXT
     return
  endif

  if (walk_num <= 0) then
     info = QMCKL_INVALID_ARG_2
     return
  endif

  if (elec_num <= 0) then
     info = QMCKL_INVALID_ARG_3
     return
  endif

  if (nucl_num <= 0) then
     info = QMCKL_INVALID_ARG_4
     return
  endif

  if (aord_num <= 0) then
     info = QMCKL_INVALID_ARG_7
     return
  endif

  factor_en = 0.0d0
  
  do nw =1, walk_num
  do a = 1, nucl_num
     do i = 1, elec_num
        x = en_distance_rescaled(nw, i, a)
        power_ser = 0.0d0

        do p = 2, aord_num
          x = x * en_distance_rescaled(nw, i, a)
          power_ser = power_ser + aord_vector(p + 1, type_nucl_vector(a)) * x
        end do

        factor_en(nw) = factor_en(nw) + aord_vector(1, type_nucl_vector(a)) *  &
                                en_distance_rescaled(nw, i, a) /               &
                                (1.0d0 + aord_vector(2, type_nucl_vector(a)) * &
                                en_distance_rescaled(nw, i, a))                &
                                + power_ser

     end do
  end do
  end do

end function qmckl_compute_factor_en_f
    #+end_src

    #+CALL: generate_c_header(table=qmckl_factor_en_args,rettyp=get_value("CRetType"),fname=get_value("Name"))

    #+RESULTS:
    #+begin_src c :tangle (eval h_func) :comments org
    qmckl_exit_code qmckl_compute_factor_en (
	  const qmckl_context context,
	  const int64_t walk_num,
	  const int64_t elec_num,
	  const int64_t nucl_num,
	  const int64_t type_nucl_num,
	  const int64_t* type_nucl_vector,
	  const int64_t aord_num,
	  const double* aord_vector,
	  const double* en_distance_rescaled,
	  double* const factor_en ); 
    #+end_src


    #+CALL: generate_c_interface(table=qmckl_factor_en_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_factor_en &
	(context, &
		 walk_num, &
		 elec_num, &
		 nucl_num, &
		 type_nucl_num, &
		 type_nucl_vector, &
		 aord_num, &
		 aord_vector, &
		 en_distance_rescaled, &
		 factor_en) &
	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 :: walk_num
      integer (c_int64_t) , intent(in)  , value :: elec_num
      integer (c_int64_t) , intent(in)  , value :: nucl_num
      integer (c_int64_t) , intent(in)  , value :: type_nucl_num
      integer (c_int64_t) , intent(in)          :: type_nucl_vector(type_nucl_num)
      integer (c_int64_t) , intent(in)  , value :: aord_num
      real    (c_double ) , intent(in)          :: aord_vector(type_nucl_num,aord_num + 1)
      real    (c_double ) , intent(in)          :: en_distance_rescaled(elec_num,nucl_num,walk_num)
      real    (c_double ) , intent(out)         :: factor_en(walk_num)

      integer(c_int32_t), external :: qmckl_compute_factor_en_f
      info = qmckl_compute_factor_en_f &
	     (context, &
		 walk_num, &
		 elec_num, &
		 nucl_num, &
		 type_nucl_num, &
		 type_nucl_vector, &
		 aord_num, &
		 aord_vector, &
		 en_distance_rescaled, &
		 factor_en)

    end function qmckl_compute_factor_en
    #+end_src

*** Test
    #+begin_src python :results output :exports none :noweb yes
import numpy as np

<<jastrow_data>>

factor_en = 0.0
for a in range(0,nucl_num):
  for i in range(0,elec_num):
    x = en_distance_rescaled[i][a]
    pow_ser = 0.0
    
    for p in range(2,aord_num+1):
      x = x * en_distance_rescaled[i][a]
      pow_ser = pow_ser + aord_vector[(p-1) + 1][type_nucl_vector[a]-1] * x

    factor_en = factor_en + aord_vector[0][type_nucl_vector[a]-1] * en_distance_rescaled[i][a] \
                          / (1.0 + aord_vector[1][type_nucl_vector[a]-1] * en_distance_rescaled[i][a])     \
                          + pow_ser
print("factor_en        :",factor_en)

    #+end_src

    #+RESULTS:
    : factor_en        : -5.865822569188727
    

     #+begin_src c :tangle (eval c_test)
/* Check if Jastrow is properly initialized */
assert(qmckl_jastrow_provided(context));

double factor_en[walk_num];
rc = qmckl_get_jastrow_factor_en(context, factor_en);

// calculate factor_en
assert(fabs(factor_en[0]+5.865822569188727) < 1.e-12);

     #+end_src

** Electron-nucleus component derivative \(f'_{en}\)
  Calculate the electron-electron jastrow component ~factor_en_deriv_e~ derivative
  with respect to the electron coordinates using the ~en_distance_rescaled~ and
 ~en_distance_rescaled_deriv_e~ which are already calculated previously.

  TODO: write equations.
  
*** Get
    #+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_jastrow_factor_en_deriv_e(qmckl_context context, double* const factor_en_deriv_e);
    #+end_src

    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
qmckl_exit_code qmckl_get_jastrow_factor_en_deriv_e(qmckl_context context, double* const factor_en_deriv_e)
{
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return QMCKL_NULL_CONTEXT;
  }

  qmckl_exit_code rc;

  rc = qmckl_provide_factor_en_deriv_e(context);
  if (rc != QMCKL_SUCCESS) return rc;

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

  int64_t sze = ctx->electron.walk_num * 4 * ctx->electron.num;
  memcpy(factor_en_deriv_e, ctx->jastrow.factor_en_deriv_e, 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_factor_en_deriv_e(qmckl_context context);
    #+end_src
    
    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
qmckl_exit_code qmckl_provide_factor_en_deriv_e(qmckl_context context)
{

  qmckl_exit_code rc;

  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return QMCKL_NULL_CONTEXT;
  }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

 /* Check if en rescaled distance is provided */
  rc = qmckl_provide_en_distance_rescaled(context);
  if(rc != QMCKL_SUCCESS) return rc;

 /* Check if en rescaled distance derivatives is provided */
  rc = qmckl_provide_en_distance_rescaled_deriv_e(context);
  if(rc != QMCKL_SUCCESS) return rc;

 /* Compute if necessary */
  if (ctx->date > ctx->jastrow.factor_en_deriv_e_date) {

 /* Allocate array */
    if (ctx->jastrow.factor_en_deriv_e == NULL) {

      qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
      mem_info.size = ctx->electron.walk_num * 4 * ctx->electron.num * sizeof(double);
      double* factor_en_deriv_e = (double*) qmckl_malloc(context, mem_info);

      if (factor_en_deriv_e == NULL) {
        return qmckl_failwith( context,
                               QMCKL_ALLOCATION_FAILED,
                               "qmckl_provide_factor_en_deriv_e",
                               NULL);
      }
      ctx->jastrow.factor_en_deriv_e = factor_en_deriv_e;
    }
    
    qmckl_exit_code rc =
      qmckl_compute_factor_en_deriv_e(context,
                              ctx->electron.walk_num,
                              ctx->electron.num,
                              ctx->nucleus.num,
                              ctx->jastrow.type_nucl_num,
                              ctx->jastrow.type_nucl_vector,
                              ctx->jastrow.aord_num,
                              ctx->jastrow.aord_vector,
                              ctx->electron.en_distance_rescaled,
                              ctx->electron.en_distance_rescaled_deriv_e,
                              ctx->jastrow.factor_en_deriv_e);
    if (rc != QMCKL_SUCCESS) {
      return rc;
    }

    ctx->jastrow.factor_en_deriv_e_date = ctx->date;
  }

  return QMCKL_SUCCESS;
}
    #+end_src

*** Compute
   :PROPERTIES:
   :Name:     qmckl_compute_factor_en_deriv_e
   :CRetType: qmckl_exit_code
   :FRetType: qmckl_exit_code
   :END:

    #+NAME: qmckl_factor_en_deriv_e_args
   | qmckl_context | context                                                       | in  | Global state                          |
   | int64_t       | walk_num                                                      | in  | Number of walkers                     |
   | int64_t       | elec_num                                                      | in  | Number of electrons                   |
   | int64_t       | nucl_num                                                      | in  | Number of nucleii                     |
   | int64_t       | type_nucl_num                                                 | in  | Number of unique nuclei               |
   | int64_t       | type_nucl_vector[type_nucl_num]                               | in  | IDs of unique nucleii                 |
   | int64_t       | aord_num                                                      | in  | Number of coefficients                |
   | double        | aord_vector[aord_num + 1][type_nucl_num]                      | in  | List of coefficients                  |
   | double        | en_distance_rescaled[walk_num][nucl_num][elec_num]            | in  | Electron-nucleus distances            |
   | double        | en_distance_rescaled_deriv_e[walk_num][4][nucl_num][elec_num] | in  | Electron-nucleus distance derivatives |
   | double        | factor_en_deriv_e[walk_num][4][elec_num]                      | out | Electron-nucleus jastrow              |

    #+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_factor_en_deriv_e_f(context, walk_num, elec_num, nucl_num, type_nucl_num,            &
                                           type_nucl_vector, aord_num, aord_vector,                         &
                                           en_distance_rescaled, en_distance_rescaled_deriv_e, factor_en_deriv_e) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer*8             , intent(in)  :: walk_num, elec_num, aord_num, nucl_num, type_nucl_num
  integer*8             , intent(in)  :: type_nucl_vector(type_nucl_num)
  double precision      , intent(in)  :: aord_vector(aord_num, nucl_num)
  double precision      , intent(in)  :: en_distance_rescaled(walk_num, elec_num, nucl_num)
  double precision      , intent(in)  :: en_distance_rescaled_deriv_e(walk_num, 4, elec_num, nucl_num)
  double precision      , intent(out) :: factor_en_deriv_e(elec_num,4,walk_num)

  integer*8 :: i, a, p, ipar, nw, ii
  double precision   :: x, spin_fact, den, invden, invden2, invden3, xinv
  double precision   :: y, lap1, lap2, lap3, third
  double precision, dimension(3) :: power_ser_g
  double precision, dimension(4) :: dx

  info = QMCKL_SUCCESS

  if (context == QMCKL_NULL_CONTEXT) then
     info = QMCKL_INVALID_CONTEXT
     return
  endif

  if (walk_num <= 0) then
     info = QMCKL_INVALID_ARG_2
     return
  endif

  if (elec_num <= 0) then
     info = QMCKL_INVALID_ARG_3
     return
  endif

  if (nucl_num <= 0) then
     info = QMCKL_INVALID_ARG_4
     return
  endif

  if (aord_num <= 0) then
     info = QMCKL_INVALID_ARG_7
     return
  endif

  factor_en_deriv_e = 0.0d0
  third = 1.0d0 / 3.0d0
  
  do nw =1, walk_num
  do a = 1, nucl_num
     do i = 1, elec_num
        x = en_distance_rescaled(nw, i, a)
        if(abs(x) < 1.0d-18) continue
        power_ser_g = 0.0d0
        den = 1.0d0 + aord_vector(2, type_nucl_vector(a)) * x
        invden = 1.0d0 / den
        invden2 = invden * invden
        invden3 = invden2 * invden
        xinv = 1.0d0 / x

        do ii = 1, 4
          dx(ii) = en_distance_rescaled_deriv_e(nw, ii, i, a)
        end do

        lap1 = 0.0d0
        lap2 = 0.0d0
        lap3 = 0.0d0
        do ii = 1, 3
          x = en_distance_rescaled(nw, i, a)
          do p = 2, aord_num
            y = p * aord_vector(p + 1, type_nucl_vector(a)) * x
            power_ser_g(ii) = power_ser_g(ii) + y * dx(ii)
            lap1 = lap1 + (p - 1) * y * xinv * dx(ii) * dx(ii)
            lap2 = lap2 + y
            x = x * en_distance_rescaled(nw, i, a)
          end do

          lap3 = lap3 - 2.0d0 * aord_vector(2, type_nucl_vector(a)) * dx(ii) * dx(ii)

          factor_en_deriv_e(i, ii, nw) = factor_en_deriv_e(i, ii, nw) + aord_vector(1, type_nucl_vector(a)) &
                                  * dx(ii) * invden2                                                        &
                                  + power_ser_g(ii)

        end do

        ii = 4
        lap2 = lap2 * dx(ii) * third
        lap3 = lap3 + den * dx(ii)
        lap3 = lap3 * aord_vector(1, type_nucl_vector(a)) * invden3
        factor_en_deriv_e(i, ii, nw) = factor_en_deriv_e(i, ii, nw) + lap1 + lap2 + lap3

     end do
  end do
  end do

end function qmckl_compute_factor_en_deriv_e_f
    #+end_src

    #+CALL: generate_c_header(table=qmckl_factor_en_deriv_e_args,rettyp=get_value("CRetType"),fname=get_value("Name"))

    #+RESULTS:
    #+begin_src c :tangle (eval h_func) :comments org
    qmckl_exit_code qmckl_compute_factor_en_deriv_e (
	  const qmckl_context context,
	  const int64_t walk_num,
	  const int64_t elec_num,
	  const int64_t nucl_num,
	  const int64_t type_nucl_num,
	  const int64_t* type_nucl_vector,
	  const int64_t aord_num,
	  const double* aord_vector,
	  const double* en_distance_rescaled,
	  const double* en_distance_rescaled_deriv_e,
	  double* const factor_en_deriv_e ); 
    #+end_src


    #+CALL: generate_c_interface(table=qmckl_factor_en_deriv_e_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_factor_en_deriv_e &
	(context, &
		 walk_num, &
		 elec_num, &
		 nucl_num, &
		 type_nucl_num, &
		 type_nucl_vector, &
		 aord_num, &
		 aord_vector, &
		 en_distance_rescaled, &
		 en_distance_rescaled_deriv_e, &
		 factor_en_deriv_e) &
	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 :: walk_num
      integer (c_int64_t) , intent(in)  , value :: elec_num
      integer (c_int64_t) , intent(in)  , value :: nucl_num
      integer (c_int64_t) , intent(in)  , value :: type_nucl_num
      integer (c_int64_t) , intent(in)          :: type_nucl_vector(type_nucl_num)
      integer (c_int64_t) , intent(in)  , value :: aord_num
      real    (c_double ) , intent(in)          :: aord_vector(type_nucl_num,aord_num + 1)
      real    (c_double ) , intent(in)          :: en_distance_rescaled(elec_num,nucl_num,walk_num)
      real    (c_double ) , intent(in)          :: en_distance_rescaled_deriv_e(elec_num,nucl_num,4,walk_num)
      real    (c_double ) , intent(out)         :: factor_en_deriv_e(elec_num,4,walk_num)

      integer(c_int32_t), external :: qmckl_compute_factor_en_deriv_e_f
      info = qmckl_compute_factor_en_deriv_e_f &
	     (context, &
		 walk_num, &
		 elec_num, &
		 nucl_num, &
		 type_nucl_num, &
		 type_nucl_vector, &
		 aord_num, &
		 aord_vector, &
		 en_distance_rescaled, &
		 en_distance_rescaled_deriv_e, &
		 factor_en_deriv_e)

    end function qmckl_compute_factor_en_deriv_e
    #+end_src

*** Test
    #+begin_src python :results output :exports none :noweb yes
import numpy as np

<<jastrow_data>>

kappa = 1.0

elec_coord = np.array(elec_coord)[0]
nucl_coord = np.array(nucl_coord)
elnuc_dist = np.zeros(shape=(elec_num, nucl_num),dtype=float)
for i in range(elec_num):
  for j in range(nucl_num):
    elnuc_dist[i, j] = np.linalg.norm(elec_coord[i] - nucl_coord[:,j])

elnuc_dist_deriv_e = np.zeros(shape=(4, elec_num, nucl_num),dtype=float)
for a in range(nucl_num):
  for i in range(elec_num):
    rij_inv = 1.0 / elnuc_dist[i, a]
    for ii in range(3):
      elnuc_dist_deriv_e[ii, i, a] = (elec_coord[i][ii] - nucl_coord[ii][a]) * rij_inv
    elnuc_dist_deriv_e[3, i, a] = 2.0 * rij_inv

en_distance_rescaled_deriv_e = np.zeros(shape=(4,elec_num,nucl_num),dtype=float)
for a in range(nucl_num):
  for i in range(elec_num):
    f = 1.0 - kappa * en_distance_rescaled[i][a]
    for ii in range(4):
      en_distance_rescaled_deriv_e[ii][i][a] = elnuc_dist_deriv_e[ii][i][a]
    en_distance_rescaled_deriv_e[3][i][a] = en_distance_rescaled_deriv_e[3][i][a] + \
                              (-kappa * en_distance_rescaled_deriv_e[0][i][a] * en_distance_rescaled_deriv_e[0][i][a]) + \
                              (-kappa * en_distance_rescaled_deriv_e[1][i][a] * en_distance_rescaled_deriv_e[1][i][a]) + \
                              (-kappa * en_distance_rescaled_deriv_e[2][i][a] * en_distance_rescaled_deriv_e[2][i][a])
    for ii in range(4):
      en_distance_rescaled_deriv_e[ii][i][a] = en_distance_rescaled_deriv_e[ii][i][a] * f

third = 1.0 / 3.0
factor_en_deriv_e = np.zeros(shape=(4,elec_num),dtype=float)
dx = np.zeros(shape=(4),dtype=float)
pow_ser_g = np.zeros(shape=(3),dtype=float)
for a in range(nucl_num):
   for i in range(elec_num):
      x = en_distance_rescaled[i][a]
      if abs(x) < 1e-18:
        continue
      pow_ser_g = np.zeros(shape=(3),dtype=float)
      den         = 1.0 + aord_vector[1][type_nucl_vector[a]-1] * x
      invden      = 1.0 / den
      invden2     = invden * invden
      invden3     = invden2 * invden
      xinv        = 1.0 / (x + 1.0E-18)

      for ii in range(4):
        dx[ii] = en_distance_rescaled_deriv_e[ii][i][a]

      lap1 = 0.0
      lap2 = 0.0
      lap3 = 0.0
      for ii in range(3):
        x = en_distance_rescaled[i][a]
        if x < 1e-18:
          continue
        for p in range(2,aord_num+1):
          y = p * aord_vector[(p-1) + 1][type_nucl_vector[a]-1] * x
          pow_ser_g[ii] = pow_ser_g[ii] + y * dx[ii]
          lap1 = lap1 + (p - 1) * y * xinv * dx[ii] * dx[ii]
          lap2 = lap2 + y
          x = x * en_distance_rescaled[i][a]
        
        lap3 = lap3 - 2.0 * aord_vector[1][type_nucl_vector[a]-1] * dx[ii] * dx[ii]

        factor_en_deriv_e[ii][i] = factor_en_deriv_e[ii][i] + aord_vector[0][type_nucl_vector[a]-1]  * \
                              dx[ii] * invden2 + pow_ser_g[ii]

      ii = 3
      lap2 = lap2 * dx[ii] * third
      lap3 = lap3 + den * dx[ii]
      lap3 = lap3 * (aord_vector[0][type_nucl_vector[a]-1] * invden3)
      factor_en_deriv_e[ii][i] = factor_en_deriv_e[ii][i] + lap1 + lap2 + lap3

print("factor_en_deriv_e[0][0]:",factor_en_deriv_e[0][0])
print("factor_en_deriv_e[1][0]:",factor_en_deriv_e[1][0])
print("factor_en_deriv_e[2][0]:",factor_en_deriv_e[2][0])
print("factor_en_deriv_e[3][0]:",factor_en_deriv_e[3][0])


    #+end_src

    #+RESULTS:
    : factor_en_deriv_e[0][0]: 0.11609919541763383
    : factor_en_deriv_e[1][0]: -0.23301394780804574
    : factor_en_deriv_e[2][0]: 0.17548337641865783
    : factor_en_deriv_e[3][0]: -0.9667363412285741
    

     #+begin_src c :tangle (eval c_test)
/* Check if Jastrow is properly initialized */
assert(qmckl_jastrow_provided(context));

// calculate factor_en_deriv_e
double factor_en_deriv_e[walk_num][4][elec_num];
rc = qmckl_get_jastrow_factor_en_deriv_e(context, &(factor_en_deriv_e[0][0][0]));

// check factor_en_deriv_e
assert(fabs(factor_en_deriv_e[0][0][0]-0.11609919541763383) < 1.e-12);
assert(fabs(factor_en_deriv_e[0][1][0]+0.23301394780804574) < 1.e-12);
assert(fabs(factor_en_deriv_e[0][2][0]-0.17548337641865783) < 1.e-12);
assert(fabs(factor_en_deriv_e[0][3][0]+0.9667363412285741 ) < 1.e-12);

     #+end_src


* End of files                                                     :noexport:

  #+begin_src c :tangle (eval h_private_type)
#endif
  #+end_src

*** Test
  #+begin_src c :tangle (eval c_test)
    rc = qmckl_context_destroy(context);
    assert (rc == QMCKL_SUCCESS);

    return 0;
}
  #+end_src


# -*- mode: org -*-
# vim: syntax=c