Jastrow Factor

1. Introduction
2. Context
3. Computation

1 Introduction

The Jastrow factor depends on the electronic (\(\mathbf{r}\)) and nuclear (\(\mathbf{R}\)) coordinates. Its defined as \(\exp(J(\mathbf{r},\mathbf{R}))\), where

\[ J(\mathbf{r},\mathbf{R}) = J_{\text{eN}}(\mathbf{r},\mathbf{R}) + J_{\text{ee}}(\mathbf{r}) + J_{\text{eeN}}(\mathbf{r},\mathbf{R}) \]

In the following, we use the notations \(r_{ij} = |\mathbf{r}_i - \mathbf{r}_j|\) and \(R_{i\alpha} = |\mathbf{r}_i - \mathbf{R}_\alpha|\).

\(J_{\text{eN}}\) contains electron-nucleus terms:

\[ J_{\text{eN}}(\mathbf{r},\mathbf{R}) = \sum_{\alpha=1}^{N_\text{nucl}} \sum_{i=1}^{N_\text{elec}} \frac{a_{1,\alpha}\, g_\alpha(R_{i\alpha})}{1+a_{2,\alpha}\, g_\alpha(R_{i\alpha})} + \sum_{p=2}^{N_\text{ord}^a} a_{p+1,\alpha}\, [g_\alpha(R_{i\alpha})]^p - J_{eN}^{\infty \alpha} \]

\(J_{\text{ee}}\) contains electron-electron terms: \[ J_{\text{ee}}(\mathbf{r}) = \sum_{i=1}^{N_\text{elec}} \sum_{j=1}^{i-1} \frac{b_1\, f(r_{ij})}{1+b_2\, f(r_{ij})} + \sum_{p=2}^{N_\text{ord}^b} a_{p+1}\, [f(r_{ij})]^p - J_{ee}^\infty \]

and \(J_{\text{eeN}}\) contains electron-electron-Nucleus terms:

\[ J_{\text{eeN}}(\mathbf{r},\mathbf{R}) = \sum_{\alpha=1}^{N_{\text{nucl}}} \sum_{i=1}^{N_{\text{elec}}} \sum_{j=1}^{i-1} \sum_{p=2}^{N_{\text{ord}}} \sum_{k=0}^{p-1} \sum_{l=0}^{p-k-2\delta_{k,0}} c_{lkp\alpha} \left[ f({r}_{ij}) \right]^k \left[ \left[ g_\alpha({R}_{i\alpha}) \right]^l + \left[ g_\alpha({R}_{j\alpha}) \right]^l \right] \left[ g_\alpha({R}_{i\,\alpha}) \, g_\alpha({R}_{j\alpha}) \right]^{(p-k-l)/2} \]

\(c_{lkp\alpha}\) are non-zero only when \(p-k-l\) is even.

\(f\) and \(g\) are scaling function defined as

\[ f(r) = \frac{1-e^{-\kappa\, r}}{\kappa} \text{ and } g_\alpha(r) = e^{-\kappa_\alpha\, r}. \]

The terms \(J_{\text{ee}}^\infty\) and \(J_{\text{eN}}^\infty\) are shifts to ensure that \(J_{\text{ee}}\) and \(J_{\text{eN}}\) have an asymptotic value of zero.

2 Context

The following data stored in the context:

Variable	Type	In/Out	Description
`uninitialized`	`int32_t`	in	Keeps bits set for uninitialized data
`rescale_factor_ee`	`double`	in	The distance scaling factor
`rescale_factor_en`	`double[type_nucl_num]`	in	The distance scaling factor
`aord_num`	`int64_t`	in	The number of a coeffecients
`bord_num`	`int64_t`	in	The number of b coeffecients
`cord_num`	`int64_t`	in	The number of c coeffecients
`type_nucl_num`	`int64_t`	in	Number of Nuclei types
`type_nucl_vector`	`int64_t[nucl_num]`	in	IDs of types of Nuclei
`a_vector`	`double[aord_num + 1][type_nucl_num]`	in	a polynomial coefficients
`b_vector`	`double[bord_num + 1]`	in	b polynomial coefficients
`c_vector`	`double[cord_num][type_nucl_num]`	in	c polynomial coefficients
`factor_ee`	`double[walker.num]`	out	Jastrow factor: electron-electron part
`factor_ee_date`	`uint64_t`	out	Jastrow factor: electron-electron part
`factor_en`	`double[walker.num]`	out	Jastrow factor: electron-nucleus part
`factor_en_date`	`uint64_t`	out	Jastrow factor: electron-nucleus part
`factor_een`	`double[walker.num]`	out	Jastrow factor: electron-electron-nucleus part
`factor_een_date`	`uint64_t`	out	Jastrow factor: electron-electron-nucleus part
`factor_ee_deriv_e`	`double[4][nelec][walker.num]`	out	Derivative of the Jastrow factor: electron-electron-nucleus part
`factor_ee_deriv_e_date`	`uint64_t`	out	Keep track of the date for the derivative
`factor_en_deriv_e`	`double[4][nelec][walker.num]`	out	Derivative of the Jastrow factor: electron-electron-nucleus part
`factor_en_deriv_e_date`	`uint64_t`	out	Keep track of the date for the en derivative
`factor_een_deriv_e`	`double[4][nelec][walker.num]`	out	Derivative of the Jastrow factor: electron-electron-nucleus part
`factor_een_deriv_e_date`	`uint64_t`	out	Keep track of the date for the een derivative

computed data:

Variable	Type	In/Out
`dim_c_vector`	`int64_t`	Number of unique C coefficients
`dim_c_vector_date`	`uint64_t`	Number of unique C coefficients
`asymp_jasa`	`double[type_nucl_num]`	Asymptotic component
`asymp_jasa_date`	`uint64_t`	Ladt modification of the asymptotic component
`asymp_jasb`	`double[2]`	Asymptotic component (up- or down-spin)
`asymp_jasb_date`	`uint64_t`	Ladt modification of the asymptotic component
`c_vector_full`	`double[dim_c_vector][nucl_num]`	vector of non-zero coefficients
`c_vector_full_date`	`uint64_t`	Keep track of changes here
`lkpm_combined_index`	`int64_t[4][dim_c_vector]`	Transform l,k,p, and m into consecutive indices
`lkpm_combined_index_date`	`uint64_t`	Transform l,k,p, and m into consecutive indices
`tmp_c`	`double[walker.num][cord_num][cord_num+1][nucl_num][elec_num]`	vector of non-zero coefficients
`dtmp_c`	`double[walker.num][elec_num][4][nucl_num][cord_num+1][cord_num]`	vector of non-zero coefficients
`ee_distance_rescaled`	`double[walker.num][num][num]`	Electron-electron rescaled distances
`ee_distance_rescaled_date`	`uint64_t`	Last modification date of the electron-electron distances
`ee_distance_rescaled_deriv_e`	`double[walker.num][4][num][num]`	Electron-electron rescaled distances derivatives
`ee_distance_rescaled_deriv_e_date`	`uint64_t`	Last modification date of the electron-electron distance derivatives
`en_distance_rescaled`	`double[walker.num][nucl_num][num]`	Electron-nucleus distances
`en_distance_rescaled_date`	`uint64_t`	Last modification date of the electron-electron distances
`en_distance_rescaled_deriv_e`	`double[walker.num][4][nucl_num][num]`	Electron-electron rescaled distances derivatives
`en_distance_rescaled_deriv_e_date`	`uint64_t`	Last modification date of the electron-electron distance derivatives
`een_rescaled_n`	`double[walker.num][cord_num+1][nucl_num][elec_num]`	The electron-electron rescaled distances raised to the powers defined by cord
`een_rescaled_n_date`	`uint64_t`	Keep track of the date of creation
`een_rescaled_e_deriv_e`	`double[walker.num][cord_num+1][elec_num][4][elec_num]`	The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons
`een_rescaled_e_deriv_e_date`	`uint64_t`	Keep track of the date of creation
`een_rescaled_n_deriv_e`	`double[walker.num][cord_num+1][nucl_num][4][elec_num]`	The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons
`een_rescaled_n_deriv_e_date`	`uint64_t`	Keep track of the date of creation

2.1 Data structure

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;
  uint64_t  asymp_jasa_date;
  uint64_t  asymp_jasb_date;
  uint64_t  tmp_c_date;
  uint64_t  dtmp_c_date;
  uint64_t  factor_ee_date;
  uint64_t  factor_en_date;
  uint64_t  factor_een_date;
  uint64_t  factor_ee_deriv_e_date;
  uint64_t  factor_en_deriv_e_date;
  uint64_t  factor_een_deriv_e_date;
  double    rescale_factor_ee;
  double*   rescale_factor_en;
  int64_t*  type_nucl_vector;
  double *  a_vector;
  double *  b_vector;
  double *  c_vector;
  double *  asymp_jasa;
  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_c_vector;
  uint64_t  dim_c_vector_date;
  double *  c_vector_full;
  uint64_t  c_vector_full_date;
  int64_t*  lkpm_combined_index;
  uint64_t  lkpm_combined_index_date;
  double *  tmp_c;
  double *  dtmp_c;
  uint64_t  ee_distance_rescaled_date;
  uint64_t  ee_distance_rescaled_deriv_e_date;
  uint64_t  en_distance_rescaled_date;
  uint64_t  en_distance_rescaled_deriv_e_date;
  double*   ee_distance_rescaled;
  double*   ee_distance_rescaled_deriv_e;
  double*   en_distance_rescaled;
  double*   en_distance_rescaled_deriv_e;
  double *  een_rescaled_e;
  double *  een_rescaled_n;
  uint64_t  een_rescaled_e_date;
  uint64_t  een_rescaled_n_date;
  double *  een_rescaled_e_deriv_e;
  double *  een_rescaled_n_deriv_e;
  uint64_t  een_rescaled_e_deriv_e_date;
  uint64_t  een_rescaled_n_deriv_e_date;
  bool      provided;
  char *    type;

  #ifdef HAVE_HPC
  bool     gpu_offload;
  #endif
} qmckl_jastrow_struct;

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.

qmckl_exit_code qmckl_init_jastrow(qmckl_context context);

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*) context;
  assert (ctx != NULL);

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

  /* Default values */
  ctx->jastrow.aord_num = -1;
  ctx->jastrow.bord_num = -1;
  ctx->jastrow.cord_num = -1;
  ctx->jastrow.type_nucl_num = -1;
  ctx->jastrow.dim_c_vector = -1;

  return QMCKL_SUCCESS;
}

2.2 Initialization functions

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

qmckl_exit_code  qmckl_set_jastrow_rescale_factor_ee (qmckl_context context, const double  kappa_ee);
qmckl_exit_code  qmckl_set_jastrow_rescale_factor_en (qmckl_context context, const double* kappa_en, const int64_t size_max);
qmckl_exit_code  qmckl_set_jastrow_aord_num          (qmckl_context context, const int64_t aord_num);
qmckl_exit_code  qmckl_set_jastrow_bord_num          (qmckl_context context, const int64_t bord_num);
qmckl_exit_code  qmckl_set_jastrow_cord_num          (qmckl_context context, 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_a_vector       (qmckl_context context, const double * a_vector, const int64_t size_max);
qmckl_exit_code  qmckl_set_jastrow_b_vector       (qmckl_context context, const double * b_vector, const int64_t size_max);
qmckl_exit_code  qmckl_set_jastrow_c_vector       (qmckl_context context, const double * c_vector, const int64_t size_max);

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.

2.2.0.1 Fortran interface

interface
   integer(qmckl_exit_code) function qmckl_set_jastrow_rescale_factor_ee (context, &
        kappa_ee) bind(C)
     use, intrinsic :: iso_c_binding
     import
     implicit none
     integer (qmckl_context) , intent(in)  , value :: context
     double precision, intent(in), value  :: kappa_ee
   end function qmckl_set_jastrow_rescale_factor_ee

   integer(qmckl_exit_code) function qmckl_set_jastrow_rescale_factor_en (context, &
        kappa_en, size_max) bind(C)
     use, intrinsic :: iso_c_binding
     import
     implicit none
     integer (qmckl_context) , intent(in)  , value :: context
     integer(c_int64_t), intent(in), value  :: size_max
     double precision, intent(in) :: kappa_en(size_max)
   end function qmckl_set_jastrow_rescale_factor_en

   integer(qmckl_exit_code) function qmckl_set_jastrow_aord_num (context, &
        aord_num) bind(C)
     use, intrinsic :: iso_c_binding
     import
     implicit none
     integer (qmckl_context) , intent(in)  , value :: context
     integer(c_int64_t), intent(in), value  :: aord_num
   end function qmckl_set_jastrow_aord_num

   integer(qmckl_exit_code) function qmckl_set_jastrow_bord_num (context, &
        bord_num) bind(C)
     use, intrinsic :: iso_c_binding
     import
     implicit none
     integer (qmckl_context) , intent(in)  , value :: context
     integer(c_int64_t), intent(in), value  :: bord_num
   end function qmckl_set_jastrow_bord_num

   integer(qmckl_exit_code) function qmckl_set_jastrow_cord_num (context, &
        cord_num) bind(C)
     use, intrinsic :: iso_c_binding
     import
     implicit none
     integer (qmckl_context) , intent(in)  , value :: context
     integer(c_int64_t), intent(in), value  :: cord_num
   end function qmckl_set_jastrow_cord_num

   integer(qmckl_exit_code) function qmckl_set_jastrow_type_nucl_num (context, &
        type_nucl_num) bind(C)
     use, intrinsic :: iso_c_binding
     import
     implicit none
     integer (qmckl_context) , intent(in)  , value :: context
     integer(c_int64_t), intent(in), value  :: type_nucl_num
   end function qmckl_set_jastrow_type_nucl_num

   integer(qmckl_exit_code) function qmckl_set_jastrow_type_nucl_vector (context, &
        type_nucl_vector, size_max) bind(C)
     use, intrinsic :: iso_c_binding
     import
     implicit none
     integer (qmckl_context) , intent(in)  , value :: context
     integer(c_int64_t), intent(in), value  :: size_max
     integer(c_int64_t), intent(in) :: type_nucl_vector(size_max)
   end function qmckl_set_jastrow_type_nucl_vector

   integer(qmckl_exit_code) function qmckl_set_jastrow_a_vector(context, &
        a_vector, size_max) bind(C)
     use, intrinsic :: iso_c_binding
     import
     implicit none
     integer (qmckl_context) , intent(in)  , value :: context
     integer(c_int64_t), intent(in), value  :: size_max
     double precision, intent(in) :: a_vector(size_max)
   end function qmckl_set_jastrow_a_vector

   integer(qmckl_exit_code) function qmckl_set_jastrow_b_vector(context, &
        b_vector, size_max) bind(C)
     use, intrinsic :: iso_c_binding
     import
     implicit none
     integer (qmckl_context) , intent(in)  , value :: context
     integer(c_int64_t), intent(in), value  :: size_max
     double precision, intent(in) :: b_vector(size_max)
   end function qmckl_set_jastrow_b_vector

   integer(qmckl_exit_code) function qmckl_set_jastrow_c_vector(context, &
        c_vector, size_max) bind(C)
     use, intrinsic :: iso_c_binding
     import
     implicit none
     integer (qmckl_context) , intent(in)  , value :: context
     integer(c_int64_t), intent(in), value  :: size_max
     double precision, intent(in) :: c_vector(size_max)
   end function qmckl_set_jastrow_c_vector

end interface

2.3 Access functions

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.

bool      qmckl_jastrow_provided           (const qmckl_context context);

2.3.0.1 Fortran interface

interface
   integer(qmckl_exit_code) function qmckl_get_jastrow_rescale_factor_ee (context, &
        kappa_ee) bind(C)
     use, intrinsic :: iso_c_binding
     import
     implicit none
     integer (qmckl_context) , intent(in)  , value :: context
     double precision, intent(out) :: kappa_ee
   end function qmckl_get_jastrow_rescale_factor_ee

   integer(qmckl_exit_code) function qmckl_get_jastrow_rescale_factor_en (context, &
        kappa_en, size_max) bind(C)
     use, intrinsic :: iso_c_binding
     import
     implicit none
     integer (qmckl_context) , intent(in), value :: context
     integer(c_int64_t), intent(in), value       :: size_max
     double precision, intent(out)               :: kappa_en(size_max)
   end function qmckl_get_jastrow_rescale_factor_en

   integer(qmckl_exit_code) function qmckl_get_jastrow_aord_num (context, &
        aord_num) bind(C)
     use, intrinsic :: iso_c_binding
     import
     implicit none
     integer (qmckl_context) , intent(in), value :: context
     integer(c_int64_t), intent(out)             :: aord_num
   end function qmckl_get_jastrow_aord_num

   integer(qmckl_exit_code) function qmckl_get_jastrow_bord_num (context, &
        bord_num) bind(C)
     use, intrinsic :: iso_c_binding
     import
     implicit none
     integer (qmckl_context) , intent(in), value :: context
     integer(c_int64_t), intent(out)             :: bord_num
   end function qmckl_get_jastrow_bord_num

   integer(qmckl_exit_code) function qmckl_get_jastrow_cord_num (context, &
        cord_num) bind(C)
     use, intrinsic :: iso_c_binding
     import
     implicit none
     integer (qmckl_context) , intent(in), value :: context
     integer(c_int64_t), intent(out)             :: cord_num
   end function qmckl_get_jastrow_cord_num

   integer(qmckl_exit_code) function qmckl_get_jastrow_type_nucl_num (context, &
        type_nucl_num) bind(C)
     use, intrinsic :: iso_c_binding
     import
     implicit none
     integer (qmckl_context) , intent(in), value :: context
     integer(c_int64_t), intent(out)              :: type_nucl_num
   end function qmckl_get_jastrow_type_nucl_num

   integer(qmckl_exit_code) function qmckl_get_jastrow_type_nucl_vector (context, &
        type_nucl_vector, size_max) bind(C)
     use, intrinsic :: iso_c_binding
     import
     implicit none
     integer (qmckl_context), intent(in), value :: context
     integer(c_int64_t), intent(in), value      :: size_max
     integer(c_int64_t), intent(out)            :: type_nucl_vector(size_max)
   end function qmckl_get_jastrow_type_nucl_vector

   integer(qmckl_exit_code) function qmckl_get_jastrow_a_vector(context, &
        a_vector, size_max) bind(C)
     use, intrinsic :: iso_c_binding
     import
     implicit none
     integer (qmckl_context) , intent(in), value :: context
     integer(c_int64_t), intent(in), value       :: size_max
     double precision, intent(out)               :: a_vector(size_max)
   end function qmckl_get_jastrow_a_vector

   integer(qmckl_exit_code) function qmckl_get_jastrow_b_vector(context, &
        b_vector, size_max) bind(C)
     use, intrinsic :: iso_c_binding
     import
     implicit none
     integer (qmckl_context) , intent(in)  , value :: context
     integer(c_int64_t), intent(in), value         :: size_max
     double precision, intent(out)                 :: b_vector(size_max)
   end function qmckl_get_jastrow_b_vector

   integer(qmckl_exit_code) function qmckl_get_jastrow_c_vector(context, &
        c_vector, size_max) bind(C)
     use, intrinsic :: iso_c_binding
     import
     implicit none
     integer (qmckl_context) , intent(in)  , value :: context
     integer(c_int64_t), intent(in), value         :: size_max
     double precision, intent(out)                 :: c_vector(size_max)
   end function qmckl_get_jastrow_c_vector

end interface

2.4 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;
int64_t nucl_num      = n2_nucl_num;
double  rescale_factor_ee   = 1.0;
double  rescale_factor_en[2] = { 1.0, 1.0 };
double* elec_coord    = &(n2_elec_coord[0][0][0]);

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

/* Provide Electron data */

qmckl_exit_code rc;

assert(!qmckl_electron_provided(context));

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

assert(qmckl_electron_provided(context));

rc = qmckl_check(context,
                 qmckl_set_electron_coord (context, 'N', walk_num, elec_coord, walk_num*3*elec_num)
                 );
assert(rc == QMCKL_SUCCESS);

double elec_coord2[walk_num*3*elec_num];

rc = qmckl_check(context,
                 qmckl_get_electron_coord (context, 'N', elec_coord2, walk_num*3*elec_num)
                 );
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_check(context,
                 qmckl_set_nucleus_num (context, nucl_num)
                 );
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_nucleus_provided(context));

double nucl_coord2[3*nucl_num];

rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2, 3*nucl_num);
assert(rc == QMCKL_NOT_PROVIDED);

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

assert(!qmckl_nucleus_provided(context));

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

rc = qmckl_check(context,
                 qmckl_get_nucleus_coord (context, 'T', nucl_coord2, nucl_num*3)
                 );
assert(rc == QMCKL_SUCCESS);
for (int64_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, nucl_num);
assert(rc == QMCKL_NOT_PROVIDED);

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

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

3 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.

3.1 Asymptotic component for \(J_{ee}\)

Calculate the asymptotic component asymp_jasb to be substracted from the final electron-electron jastrow factor \(J_{\text{ee}}\). The asymptotic component is calculated via the b_vector and the electron-electron rescale factor rescale_factor_ee.

\[ J_{\text{ee}}^{\infty} = \frac{b_1 \kappa^{-1}}{1 + b_2 \kappa^{-1}} \]

3.1.1 Get

qmckl_exit_code
qmckl_get_jastrow_asymp_jasb(qmckl_context context,
                             double* const asymp_jasb,
                             const int64_t size_max);

3.1.1.1 Fortran interface

interface
   integer(qmckl_exit_code) function qmckl_get_jastrow_asymp_jasb(context, &
        asymp_jasb, size_max) bind(C)
     use, intrinsic :: iso_c_binding
     import
     implicit none
     integer (qmckl_context) , intent(in), value :: context
     integer(c_int64_t), intent(in), value       :: size_max
     double precision, intent(out)               :: asymp_jasb(size_max)
   end function qmckl_get_jastrow_asymp_jasb
end interface

3.1.2 Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`bord_num`	`int64_t`	in	Order of the polynomial
`b_vector`	`double[bord_num+1]`	in	Values of b
`rescale_factor_ee`	`double`	in	Electron coordinates
`asymp_jasb`	`double[2]`	out	Asymptotic value

integer function qmckl_compute_jastrow_asymp_jasb_f(context, bord_num, b_vector, rescale_factor_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)  :: b_vector(bord_num + 1)
  double precision      , intent(in)  :: rescale_factor_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_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 = b_vector(1) * kappa_inv / (1.0d0 + b_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) + b_vector(p + 1) * x
     end do
  end do

end function qmckl_compute_jastrow_asymp_jasb_f

qmckl_exit_code qmckl_compute_jastrow_asymp_jasb (
          const qmckl_context context,
          const int64_t bord_num,
          const double* b_vector,
          const double rescale_factor_ee,
          double* const asymp_jasb ) {

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

    if (bord_num < 0) {
      return QMCKL_INVALID_ARG_2;
    }

    const double kappa_inv = 1.0 / rescale_factor_ee;
    const double asym_one = b_vector[0] * kappa_inv / (1.0 + b_vector[1] * kappa_inv);
    asymp_jasb[0] = asym_one;
    asymp_jasb[1] = 0.5 * asym_one;

    for (int i = 0 ; i <= 1; ++i) {
      double x = kappa_inv;
      for (int p = 1; p < bord_num; ++p){
        x *= kappa_inv;
        asymp_jasb[i] = asymp_jasb[i] + b_vector[p + 1] * x;
      }
    }

  return QMCKL_SUCCESS;
}

3.1.3 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* a_vector = &(n2_a_vector[0][0]);
double* b_vector = &(n2_b_vector[0]);
double* c_vector = &(n2_c_vector[0][0]);
int64_t dim_c_vector=0;

assert(!qmckl_jastrow_provided(context));

/* Set the data */
rc = qmckl_check(context,
                 qmckl_set_jastrow_aord_num(context, aord_num)
                 );
rc = qmckl_check(context,
                 qmckl_set_jastrow_bord_num(context, bord_num)
                 );
rc = qmckl_check(context,
                 qmckl_set_jastrow_cord_num(context, cord_num)
                 );
assert(rc == QMCKL_SUCCESS);
rc = qmckl_check(context,
                 qmckl_set_jastrow_type_nucl_num(context, type_nucl_num)
                 );
assert(rc == QMCKL_SUCCESS);
rc = qmckl_check(context,
                 qmckl_set_jastrow_type_nucl_vector(context, type_nucl_vector, nucl_num)
                 );
assert(rc == QMCKL_SUCCESS);
rc = qmckl_check(context,
                 qmckl_set_jastrow_a_vector(context, a_vector,(aord_num+1)*type_nucl_num)
                 );
assert(rc == QMCKL_SUCCESS);
rc = qmckl_check(context,
                 qmckl_set_jastrow_b_vector(context, b_vector,(bord_num+1))
                 );
assert(rc == QMCKL_SUCCESS);
rc = qmckl_check(context,
                 qmckl_get_jastrow_dim_c_vector(context, &dim_c_vector)
                 );
assert(rc == QMCKL_SUCCESS);
rc = qmckl_check(context,
                 qmckl_set_jastrow_c_vector(context, c_vector,dim_c_vector*type_nucl_num)
                 );
assert(rc == QMCKL_SUCCESS);

double k_ee = 0.;
double k_en[2] = { 0., 0. };
rc = qmckl_check(context,
                 qmckl_set_jastrow_rescale_factor_en(context, rescale_factor_en, type_nucl_num)
                 );
assert(rc == QMCKL_SUCCESS);

rc = qmckl_check(context,
                 qmckl_set_jastrow_rescale_factor_ee(context, rescale_factor_ee)
                 );
assert(rc == QMCKL_SUCCESS);

rc = qmckl_check(context,
                 qmckl_get_jastrow_rescale_factor_ee (context, &k_ee)
                 );
assert(rc == QMCKL_SUCCESS);
assert(k_ee == rescale_factor_ee);

rc = qmckl_check(context,
                 qmckl_get_jastrow_rescale_factor_en (context, &(k_en[0]), type_nucl_num)
                 );
assert(rc == QMCKL_SUCCESS);
for (int i=0 ; i<type_nucl_num ; ++i) {
  assert(k_en[i] == rescale_factor_en[i]);
 }

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

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

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

3.2 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

\(\eta\) is the spin factor, \(B\) is the vector of \(b\) parameters, \(C\) is the array of scaled distances.

3.2.1 Get

qmckl_exit_code
qmckl_get_jastrow_factor_ee(qmckl_context context,
                            double* const factor_ee,
                            const int64_t size_max);

3.2.1.1 Fortran interface

interface
   integer(qmckl_exit_code) function qmckl_get_jastrow_factor_ee (context, &
        factor_ee, size_max) bind(C)
     use, intrinsic :: iso_c_binding
     import
     implicit none
     integer (qmckl_context) , intent(in), value :: context
     integer(c_int64_t), intent(in), value       :: size_max
     double precision, intent(out)               :: factor_ee(size_max)
   end function qmckl_get_jastrow_factor_ee
end interface

3.2.2 Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`walk_num`	`int64_t`	in	Number of walkers
`elec_num`	`int64_t`	in	Number of electrons
`up_num`	`int64_t`	in	Number of alpha electrons
`bord_num`	`int64_t`	in	Number of coefficients
`b_vector`	`double[bord_num+1]`	in	List of coefficients
`ee_distance_rescaled`	`double[walk_num][elec_num][elec_num]`	in	Electron-electron distances
`asymp_jasb`	`double[2]`	in	Electron-electron distances
`factor_ee`	`double[walk_num]`	out	Electron-electron distances

integer function qmckl_compute_factor_ee_f(context, walk_num, elec_num, up_num, bord_num,            &
                                           b_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)  :: b_vector(bord_num + 1)
  double precision      , intent(in)  :: ee_distance_rescaled(elec_num, elec_num, walk_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   :: x, power_ser, spin_fact

  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(i,j,nw)
           power_ser = 0.0d0
           spin_fact = 1.0d0
           ipar = 1

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

           if(j <= up_num .or. i > up_num) then
              spin_fact = 0.5d0
              ipar = 2
           endif

           factor_ee(nw) = factor_ee(nw) + spin_fact * b_vector(1)  * &
                ee_distance_rescaled(i,j,nw) / &
                (1.0d0 + b_vector(2) *   &
                ee_distance_rescaled(i,j,nw))  &
                + power_ser - asymp_jasb(ipar) 

        end do
     end do
  end do

end function qmckl_compute_factor_ee_f

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* b_vector,
      const double* ee_distance_rescaled,
      const double* asymp_jasb,
      double* const factor_ee ) {

  int ipar;
  double x, x1, spin_fact, power_ser;

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

  if (walk_num <= 0) {
     return QMCKL_INVALID_ARG_2;
  }

  if (elec_num <= 0) {
     return QMCKL_INVALID_ARG_3;
  }

  if (bord_num < 0) {
     return QMCKL_INVALID_ARG_4;
  }

  for (int nw = 0; nw < walk_num; ++nw) {
    factor_ee[nw] = 0.0; // put init array here.
    for (int i = 0; i < elec_num; ++i ) {
      for (int j = 0; j < i; ++j) {
        //x = ee_distance_rescaled[j * (walk_num * elec_num) + i * (walk_num) + nw];
        x = ee_distance_rescaled[j + i * elec_num + nw*(elec_num * elec_num)];
        x1 = x;
        power_ser = 0.0;
        spin_fact = 1.0;
        ipar = 0; // index of asymp_jasb

        for (int p = 1; p < bord_num; ++p) {
          x = x * x1;
          power_ser = power_ser + b_vector[p + 1] * x;
        }

        if(i < up_num || j >= up_num) {
          spin_fact = 0.5;
          ipar = 1;
        }

        factor_ee[nw] = factor_ee[nw] + spin_fact * b_vector[0]  * 
                                x1 / (1.0 + b_vector[1] * x1)   
                               - asymp_jasb[ipar] + power_ser;

      }
    }
  }

  return QMCKL_SUCCESS;
}

3.2.3 Test

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

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

// calculate factor_ee
printf("%e\n%e\n\n",factor_ee[0],-4.282760865958113 );
assert(fabs(factor_ee[0]+4.282760865958113) < 1.e-12);

3.3 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.

3.3.1 Get

qmckl_exit_code
qmckl_get_jastrow_factor_ee_deriv_e(qmckl_context context,
                                    double* const factor_ee_deriv_e,
                                    const int64_t size_max);

3.3.2 Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`walk_num`	`int64_t`	in	Number of walkers
`elec_num`	`int64_t`	in	Number of electrons
`up_num`	`int64_t`	in	Number of alpha electrons
`bord_num`	`int64_t`	in	Number of coefficients
`b_vector`	`double[bord_num+1]`	in	List of coefficients
`ee_distance_rescaled`	`double[walk_num][elec_num][elec_num]`	in	Electron-electron distances
`ee_distance_rescaled_deriv_e`	`double[walk_num][4][elec_num][elec_num]`	in	Electron-electron distances
`factor_ee_deriv_e`	`double[walk_num][4][elec_num]`	out	Electron-electron distances

integer function qmckl_compute_factor_ee_deriv_e_doc_f( &
     context, walk_num, elec_num, up_num, bord_num, &
     b_vector, ee_distance_rescaled, ee_distance_rescaled_deriv_e,  &
     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)  :: b_vector(bord_num + 1)
  double precision      , intent(in)  :: ee_distance_rescaled(elec_num, elec_num,walk_num)
  double precision      , intent(in)  :: ee_distance_rescaled_deriv_e(4,elec_num, elec_num,walk_num)   !TODO
  double precision      , intent(out) :: factor_ee_deriv_e(elec_num,4,walk_num)

  integer*8 :: i, j, p, 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(i,j,nw)
        if(abs(x) < 1.0d-18) cycle
        pow_ser_g   = 0.0d0
        spin_fact   = 1.0d0
        den         = 1.0d0 + b_vector(2) * x
        invden      = 1.0d0 / den
        invden2     = invden * invden
        invden3     = invden2 * invden
        xinv        = 1.0d0 / (x + 1.0d-18)

        dx(1) = ee_distance_rescaled_deriv_e(1, i, j, nw)
        dx(2) = ee_distance_rescaled_deriv_e(2, i, j, nw)
        dx(3) = ee_distance_rescaled_deriv_e(3, i, j, nw)
        dx(4) = ee_distance_rescaled_deriv_e(4, i, j, nw)

        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(i, j, nw)
          if(abs(x) < 1.0d-18) cycle
          do p = 2, bord_num
            y = p * b_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(i, j, nw)
          end do

          lap3 = lap3 - 2.0d0 * b_vector(2) * dx(ii) * dx(ii)

          factor_ee_deriv_e( j, ii, nw) = factor_ee_deriv_e( j, ii, nw) + spin_fact * b_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 * b_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_doc_f

qmckl_exit_code qmckl_compute_factor_ee_deriv_e_hpc(
      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* b_vector,
      const double* ee_distance_rescaled,
      const double* ee_distance_rescaled_deriv_e,
      double* const factor_ee_deriv_e ) {

  int64_t ii;
  double  pow_ser_g[3];
  double  dx[4];
  double  x, spin_fact, y;
  double  den, invden, invden2, invden3, xinv;
  double  lap1, lap2, lap3, third;

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

  if (walk_num <= 0) {
    return QMCKL_INVALID_ARG_2;
  } 

  if (elec_num <= 0) {
    return QMCKL_INVALID_ARG_3;
  }

  if (bord_num < 0) {
    return QMCKL_INVALID_ARG_4;
  }


  for (int nw = 0; nw < walk_num; ++nw) {
    for (int ii = 0; ii < 4; ++ii) {
      for (int j = 0; j < elec_num; ++j) {
        factor_ee_deriv_e[j + ii * elec_num + nw * elec_num * 4]  = 0.0;
      }
    }
  }

  third = 1.0 / 3.0;

  for (int nw = 0; nw < walk_num; ++nw) {
    for (int i = 0; i < elec_num; ++i) {
      for (int j = 0; j < elec_num; ++j) {
        x = ee_distance_rescaled[j + i * elec_num + nw * elec_num * elec_num];
        if (fabs(x) < 1.0e-18) continue;
        for (int ii = 0; ii < 3; ++ii){
            pow_ser_g[ii] = 0.0;
        }    
        spin_fact   = 1.0;
        den         = 1.0 + b_vector[1] * x;
        invden      = 1.0 / den;
        invden2     = invden * invden;
        invden3     = invden2 * invden;
        xinv        = 1.0 / (x + 1.0e-18);

        dx[0] = ee_distance_rescaled_deriv_e[0 \
                                           + j * 4 + i * 4 * elec_num \
                                           + nw * 4 * elec_num * elec_num];
        dx[1] = ee_distance_rescaled_deriv_e[1  \
                                           + j * 4 + i * 4 * elec_num \
                                           + nw * 4 * elec_num * elec_num];
        dx[2] = ee_distance_rescaled_deriv_e[2  \
                                           + j * 4 + i * 4 * elec_num \
                                           + nw * 4 * elec_num * elec_num];
        dx[3] = ee_distance_rescaled_deriv_e[3  \
                                           + j * 4 + i * 4 * elec_num \
                                           + nw * 4 * elec_num * elec_num];

        if((i <= (up_num-1) && j <= (up_num-1) ) || (i > (up_num-1) && j > (up_num-1))) {
          spin_fact = 0.5;
        }

        lap1 = 0.0;
        lap2 = 0.0;
        lap3 = 0.0;
        for (int ii = 0; ii < 3; ++ii) {
          x = ee_distance_rescaled[j + i * elec_num + nw * elec_num * elec_num];
          if (fabs(x) < 1.0e-18) continue;
          for (int p = 2; p < bord_num+1; ++p) {
            y = p * b_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 * elec_num + nw * elec_num * elec_num];
          }

          lap3 = lap3 - 2.0 * b_vector[1] * dx[ii] * dx[ii];

          factor_ee_deriv_e[i  + ii * elec_num  + nw * elec_num * 4 ] +=              \
                                       + spin_fact * b_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 * b_vector[0] * invden3);
        factor_ee_deriv_e[i + ii*elec_num + nw * elec_num * 4] += lap1 + lap2 + lap3;

      }
    }
  }

  return QMCKL_SUCCESS;
}

qmckl_exit_code qmckl_compute_factor_ee_deriv_e_hpc (
  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* b_vector,
  const double* ee_distance_rescaled,
  const double* ee_distance_rescaled_deriv_e,
  double* const factor_ee_deriv_e );

qmckl_exit_code qmckl_compute_factor_ee_deriv_e_doc (
  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* b_vector,
  const double* ee_distance_rescaled,
  const double* ee_distance_rescaled_deriv_e,
  double* const factor_ee_deriv_e );

    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* b_vector,
      const double* ee_distance_rescaled,
      const double* ee_distance_rescaled_deriv_e,
      double* const factor_ee_deriv_e ) {

      #ifdef HAVE_HPC
      return qmckl_compute_factor_ee_deriv_e_hpc(context, walk_num, elec_num, up_num, bord_num, b_vector, ee_distance_rescaled, ee_distance_rescaled_deriv_e, factor_ee_deriv_e ); 
      #else
      return qmckl_compute_factor_ee_deriv_e_doc(context, walk_num, elec_num, up_num, bord_num, b_vector, ee_distance_rescaled, ee_distance_rescaled_deriv_e, factor_ee_deriv_e ); 
      #endif
}

3.3.3 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]),walk_num*4*elec_num);

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

3.4 Asymptotic component for \(J_{eN}\)

Calculate the asymptotic component asymp_jasa to be substracted from the final electron-nucleus jastrow factor \(J_{\text{eN}}\). The asymptotic component is calculated via the a_vector and the electron-nucleus rescale factors rescale_factor_en.

\[ J_{\text{en}}^{\infty \alpha} = \frac{a_1 \kappa_\alpha^{-1}}{1 + a_2 \kappa_\alpha^{-1}} \]

3.4.1 Get

qmckl_exit_code
qmckl_get_jastrow_asymp_jasa(qmckl_context context,
                             double* const asymp_jasa,
                             const int64_t size_max);

3.4.1.1 Fortran interface

interface
   integer(qmckl_exit_code) function qmckl_get_jastrow_asymp_jasa(context, &
        asymp_jasa, size_max) bind(C)
     use, intrinsic :: iso_c_binding
     import
     implicit none
     integer (qmckl_context) , intent(in), value :: context
     integer(c_int64_t), intent(in), value       :: size_max
     double precision, intent(out)               :: asymp_jasa(size_max)
   end function qmckl_get_jastrow_asymp_jasa
end interface

3.4.2 Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`aord_num`	`int64_t`	in	Order of the polynomial
`type_nucl_num`	`int64_t`	in	Number of nucleus types
`a_vector`	`double[type_nucl_num][aord_num+1]`	in	Values of a
`rescale_factor_en`	`double[type_nucl_num]`	in	Electron nucleus distances
`asymp_jasa`	`double[type_nucl_num]`	out	Asymptotic value

integer function qmckl_compute_jastrow_asymp_jasa_f(context, aord_num, type_nucl_num, a_vector, &
     rescale_factor_en, asymp_jasa) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer*8             , intent(in)  :: aord_num
  integer*8             , intent(in)  :: type_nucl_num
  double precision      , intent(in)  :: a_vector(aord_num + 1, type_nucl_num)
  double precision      , intent(in)  :: rescale_factor_en(type_nucl_num)
  double precision      , intent(out) :: asymp_jasa(type_nucl_num)

  integer*8 :: i, j, p
  double precision   :: kappa_inv, x, asym_one


  info = QMCKL_SUCCESS

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

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

  do i=1,type_nucl_num

     kappa_inv = 1.0d0 / rescale_factor_en(i)

     asymp_jasa(i) = a_vector(1,i) * kappa_inv / (1.0d0 + a_vector(2,i) * kappa_inv)

     x = kappa_inv
     do p = 2, aord_num
        x = x * kappa_inv
        asymp_jasa(i) = asymp_jasa(i) + a_vector(p+1, i) * x
     end do

  end do

end function qmckl_compute_jastrow_asymp_jasa_f

/*
qmckl_exit_code qmckl_compute_jastrow_asymp_jasa (
          const qmckl_context context,
          const int64_t aord_num,
          const int64_t type_nucl_num,
          const double* a_vector,
          double* const rescale_factor_en,
          double* const asymp_jasa ) {

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

    if (aord_num < 0) {
      return QMCKL_INVALID_ARG_2;
    }

    for (int i = 0 ; i <= type_nucl_num; ++i) {
      const double kappa_inv = 1.0 / rescale_factor_en[i];
      asymp_jasa[i] =  a_vector[aord_num*i] * kappa_inv / (1.0 + a_vector[1 + aord_num*i] * kappa_inv);

      double x = kappa_inv;
      for (int p = 1; p < aord_num; ++p){
        x *= kappa_inv;
        asymp_jasa[i] = asymp_jasa[i] + a_vector[p + 1 + aord_num*i] * x;
      }
    }

  return QMCKL_SUCCESS;
}
*/

qmckl_exit_code qmckl_compute_jastrow_asymp_jasa (
      const qmckl_context context,
      const int64_t aord_num,
      const int64_t type_nucl_num,
      const double* a_vector,
      const double* rescale_factor_en,
      double* const asymp_jasa );

3.4.3 Test

double asymp_jasa[2];
rc = qmckl_get_jastrow_asymp_jasa(context, asymp_jasa, type_nucl_num);

// calculate asymp_jasb
printf("%e %e\n", asymp_jasa[0], -0.548554);
assert(fabs(-0.548554 - asymp_jasa[0]) < 1.e-12);

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

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

\[ f_{en} = \sum_{i,j

3.5.1 Get

qmckl_exit_code
qmckl_get_jastrow_factor_en(qmckl_context context,
                            double* const factor_en,
                            const int64_t size_max);

3.5.1.1 Fortran interface

interface
   integer(qmckl_exit_code) function qmckl_get_jastrow_factor_en (context, &
        factor_en, size_max) bind(C)
     use, intrinsic :: iso_c_binding
     import
     implicit none
     integer (qmckl_context) , intent(in), value :: context
     integer(c_int64_t), intent(in), value       :: size_max
     double precision, intent(out)               :: factor_en(size_max)
   end function qmckl_get_jastrow_factor_en
end interface

3.5.2 Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`walk_num`	`int64_t`	in	Number of walkers
`elec_num`	`int64_t`	in	Number of electrons
`nucl_num`	`int64_t`	in	Number of nucleii
`type_nucl_num`	`int64_t`	in	Number of unique nuclei
`type_nucl_vector`	`int64_t[nucl_num]`	in	IDs of unique nucleii
`aord_num`	`int64_t`	in	Number of coefficients
`a_vector`	`double[aord_num+1][type_nucl_num]`	in	List of coefficients
`en_distance_rescaled`	`double[walk_num][nucl_num][elec_num]`	in	Electron-nucleus distances
`asymp_jasa`	`double[type_nucl_num]`	in	Type of nuclei
`factor_en`	`double[walk_num]`	out	Electron-nucleus jastrow

integer function qmckl_compute_factor_en_f( &
     context, walk_num, elec_num, nucl_num, type_nucl_num, &
     type_nucl_vector, aord_num, a_vector, &
     en_distance_rescaled, asymp_jasa, 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(nucl_num)
  double precision      , intent(in)  :: a_vector(aord_num + 1, type_nucl_num)
  double precision      , intent(in)  :: en_distance_rescaled(elec_num, nucl_num, walk_num)
  double precision      , intent(in)  :: asymp_jasa(type_nucl_num)
  double precision      , intent(out) :: factor_en(walk_num)

  integer*8 :: i, a, p, nw
  double precision   :: x, 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


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

           factor_en(nw) = factor_en(nw) + a_vector(1, type_nucl_vector(a)) * x / &
                (1.0d0 + a_vector(2, type_nucl_vector(a)) * x) - asymp_jasa(type_nucl_vector(a))

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

        end do
     end do
  end do

end function qmckl_compute_factor_en_f

/*
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* a_vector,
      const double* en_distance_rescaled,
      const double* asymp_jasa,
      double* const factor_en ) {

  double  x, x1, power_ser;


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

  if (walk_num <= 0) {
     return QMCKL_INVALID_ARG_2;
  }

  if (elec_num <= 0) {
     return QMCKL_INVALID_ARG_3;
  }

  if (nucl_num <= 0) {
     return QMCKL_INVALID_ARG_4;
  }

  if (type_nucl_num <= 0) {
    return QMCKL_INVALID_ARG_5;
  }

  if (type_nucl_vector == NULL) {
    return QMCKL_INVALID_ARG_6;
  }

  if (aord_num < 0) {
     return QMCKL_INVALID_ARG_7;
  }

  if (a_vector == NULL) {
    return QMCKL_INVALID_ARG_8;
  }

  if (en_distance_rescaled == NULL) {
    return QMCKL_INVALID_ARG_9;
  }

  if (factor_en == NULL) {
    return QMCKL_INVALID_ARG_10;
  }


  for (int nw = 0; nw < walk_num; ++nw ) {
    // init array
    factor_en[nw] = 0.0;
    for (int a = 0; a < nucl_num; ++a ) {
      for (int i = 0; i < elec_num; ++i ) {
        x = en_distance_rescaled[i + a * elec_num + nw * (elec_num * nucl_num)];
        x1 = x;
        power_ser = 0.0;

        for (int p = 2; p < aord_num+1; ++p) {
          x = x * x1;
          power_ser = power_ser + a_vector[p+ (type_nucl_vector[a]-1) * aord_num] * x;
        }

        factor_en[nw] = factor_en[nw] + a_vector[0 + (type_nucl_vector[a]-1)*aord_num] * x1 / \
                        (1.0 + a_vector[1 + (type_nucl_vector[a]-1) * aord_num] * x1) + \
                        power_ser;

      }
      factor_en[nw] = factor_en[nw] + asymp_jasa[type_nucl_vector[a];
    }
  }

  return QMCKL_SUCCESS;
}
*/

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* a_vector,
      const double* en_distance_rescaled,
      const double* asymp_jasa,
      double* const factor_en );

3.5.3 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,walk_num);

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

3.6 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.

3.6.1 Get

qmckl_exit_code
qmckl_get_jastrow_factor_en_deriv_e(qmckl_context context,
                                    double* const factor_en_deriv_e,
                                    const int64_t size_max);

3.6.2 Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`walk_num`	`int64_t`	in	Number of walkers
`elec_num`	`int64_t`	in	Number of electrons
`nucl_num`	`int64_t`	in	Number of nucleii
`type_nucl_num`	`int64_t`	in	Number of unique nuclei
`type_nucl_vector`	`int64_t[nucl_num]`	in	IDs of unique nucleii
`aord_num`	`int64_t`	in	Number of coefficients
`a_vector`	`double[aord_num+1][type_nucl_num]`	in	List of coefficients
`en_distance_rescaled`	`double[walk_num][nucl_num][elec_num]`	in	Electron-nucleus distances
`en_distance_rescaled_deriv_e`	`double[walk_num][4][nucl_num][elec_num]`	in	Electron-nucleus distance derivatives
`factor_en_deriv_e`	`double[walk_num][4][elec_num]`	out	Electron-nucleus jastrow

integer function qmckl_compute_factor_en_deriv_e_f( &
     context, walk_num, elec_num, nucl_num, type_nucl_num, &
     type_nucl_vector, aord_num, a_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(nucl_num)
  double precision      , intent(in)  :: a_vector(aord_num + 1, type_nucl_num)
  double precision      , intent(in)  :: en_distance_rescaled(elec_num, nucl_num, walk_num)
  double precision      , intent(in)  :: en_distance_rescaled_deriv_e(4, elec_num, nucl_num, walk_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, 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(i,a,nw)
        if(abs(x) < 1.0d-18) continue
        power_ser_g = 0.0d0
        den = 1.0d0 + a_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(ii,i,a,nw)
        end do

        lap1 = 0.0d0
        lap2 = 0.0d0
        lap3 = 0.0d0
        do ii = 1, 3
          x = en_distance_rescaled(i, a, nw)
          do p = 2, aord_num
            y = p * a_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(i, a, nw)
          end do

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

          factor_en_deriv_e(i, ii, nw) = factor_en_deriv_e(i, ii, nw) + a_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 * a_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

3.6.3 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]),walk_num*4*elec_num);

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

3.7 Electron-electron rescaled distances

ee_distance_rescaled stores the matrix of the rescaled distances between all pairs of electrons:

\[ C_{ij} = \left( 1 - \exp{-\kappa C_{ij}}\right)/\kappa \]

where \(C_{ij}\) is the matrix of electron-electron distances.

3.7.1 Get

qmckl_exit_code qmckl_get_jastrow_ee_distance_rescaled(qmckl_context context, double* const distance_rescaled);

3.7.2 Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`elec_num`	`int64_t`	in	Number of electrons
`rescale_factor_ee`	`double`	in	Factor to rescale ee distances
`walk_num`	`int64_t`	in	Number of walkers
`coord`	`double[3][walk_num][elec_num]`	in	Electron coordinates
`ee_distance`	`double[walk_num][elec_num][elec_num]`	out	Electron-electron rescaled distances

integer function qmckl_compute_ee_distance_rescaled_f(context, elec_num, rescale_factor_ee, walk_num, &
     coord, ee_distance_rescaled) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer*8             , intent(in)  :: elec_num
  double precision      , intent(in)  :: rescale_factor_ee
  integer*8             , intent(in)  :: walk_num
  double precision      , intent(in)  :: coord(elec_num,walk_num,3)
  double precision      , intent(out) :: ee_distance_rescaled(elec_num,elec_num,walk_num)

  integer*8 :: k

  info = QMCKL_SUCCESS

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

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

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

  do k=1,walk_num
     info = qmckl_distance_rescaled(context, 'T', 'T', elec_num, elec_num, &
          coord(1,k,1), elec_num * walk_num, &
          coord(1,k,1), elec_num * walk_num, &
          ee_distance_rescaled(1,1,k), elec_num, rescale_factor_ee)
     if (info /= QMCKL_SUCCESS) then
        exit
     endif
  end do

end function qmckl_compute_ee_distance_rescaled_f

3.7.3 Test

assert(qmckl_electron_provided(context));


double ee_distance_rescaled[walk_num * elec_num * elec_num];
rc = qmckl_get_jastrow_ee_distance_rescaled(context, ee_distance_rescaled);

// (e1,e2,w)
// (0,0,0) == 0.
assert(ee_distance_rescaled[0] == 0.);

// (1,0,0) == (0,1,0)
assert(ee_distance_rescaled[1] == ee_distance_rescaled[elec_num]);

// value of (1,0,0)
assert(fabs(ee_distance_rescaled[1]-0.5502278003524018) < 1.e-12);

// (0,0,1) == 0.
assert(ee_distance_rescaled[5*elec_num + 5] == 0.);

// (1,0,1) == (0,1,1)
assert(ee_distance_rescaled[5*elec_num+6] == ee_distance_rescaled[6*elec_num+5]);

// value of (1,0,1)
assert(fabs(ee_distance_rescaled[5*elec_num+6]-0.3622098222364193) < 1.e-12);

3.8 Electron-electron rescaled distance gradients and Laplacian with respect to electron coords

The rescaled distances which is given as \(R = (1 - \exp{-\kappa r})/\kappa\) needs to be perturbed with respect to the electorn coordinates. This data is stored in the ee_distance_rescaled_deriv_e tensor. The The first three elements of this three index tensor [4][num][num] gives the derivatives in the x, y, and z directions \(dx, dy, dz\) and the last index gives the Laplacian \(\partial x^2 + \partial y^2 + \partial z^2\).

3.8.1 Get

qmckl_exit_code qmckl_get_jastrow_ee_distance_rescaled_deriv_e(qmckl_context context, double* const distance_rescaled_deriv_e);

3.8.2 Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`elec_num`	`int64_t`	in	Number of electrons
`rescale_factor_ee`	`double`	in	Factor to rescale ee distances
`walk_num`	`int64_t`	in	Number of walkers
`coord`	`double[3][walk_num][elec_num]`	in	Electron coordinates
`ee_distance_deriv_e`	`double[walk_num][4][elec_num][elec_num]`	out	Electron-electron rescaled distance derivatives

integer function qmckl_compute_ee_distance_rescaled_deriv_e_f(context, elec_num, rescale_factor_ee, walk_num, &
     coord, ee_distance_rescaled_deriv_e) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer*8             , intent(in)  :: elec_num
  double precision      , intent(in)  :: rescale_factor_ee
  integer*8             , intent(in)  :: walk_num
  double precision      , intent(in)  :: coord(elec_num,walk_num,3)
  double precision      , intent(out) :: ee_distance_rescaled_deriv_e(4,elec_num,elec_num,walk_num)

  integer*8 :: k

  info = QMCKL_SUCCESS

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

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

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

  do k=1,walk_num
     info = qmckl_distance_rescaled_deriv_e(context, 'T', 'T', elec_num, elec_num, &
          coord(1,k,1), elec_num*walk_num, &
          coord(1,k,1), elec_num*walk_num, &
          ee_distance_rescaled_deriv_e(1,1,1,k), elec_num, rescale_factor_ee)
     if (info /= QMCKL_SUCCESS) then
        exit
     endif
  end do

end function qmckl_compute_ee_distance_rescaled_deriv_e_f

3.8.3 Test

assert(qmckl_electron_provided(context));


double ee_distance_rescaled_deriv_e[4 * walk_num * elec_num * elec_num];
rc = qmckl_get_jastrow_ee_distance_rescaled_deriv_e(context, ee_distance_rescaled_deriv_e);

// TODO: Get exact values
//// (e1,e2,w)
//// (0,0,0) == 0.
//assert(ee_distance[0] == 0.);
//
//// (1,0,0) == (0,1,0)
//assert(ee_distance[1] == ee_distance[elec_num]);
//
//// value of (1,0,0)
//assert(fabs(ee_distance[1]-7.152322512964209) < 1.e-12);
//
//// (0,0,1) == 0.
//assert(ee_distance[elec_num*elec_num] == 0.);
//
//// (1,0,1) == (0,1,1)
//assert(ee_distance[elec_num*elec_num+1] == ee_distance[elec_num*elec_num+elec_num]);
//
//// value of (1,0,1)
//assert(fabs(ee_distance[elec_num*elec_num+1]-6.5517646321055665) < 1.e-12);

3.9 Electron-electron-nucleus rescaled distances for each order

een_rescaled_e stores the table of the rescaled distances between all pairs of electrons and raised to the power \(p\) defined by cord_num:

\[ C_{ij,p} = \left( 1 - \exp{-\kappa C_{ij}} \right)^p \]

where \(C_{ij}\) is the matrix of electron-electron distances.

3.9.1 Get

qmckl_exit_code
qmckl_get_jastrow_een_rescaled_e(qmckl_context context,
                                 double* const distance_rescaled,
                                 const int64_t size_max);

3.9.2 Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`walk_num`	`int64_t`	in	Number of walkers
`elec_num`	`int64_t`	in	Number of electrons
`cord_num`	`int64_t`	in	Order of polynomials
`rescale_factor_ee`	`double`	in	Factor to rescale ee distances
`ee_distance`	`double[walk_num][elec_num][elec_num]`	in	Electron-electron distances
`een_rescaled_e`	`double[walk_num][0:cord_num][elec_num][elec_num]`	out	Electron-electron rescaled distances

integer function qmckl_compute_een_rescaled_e_doc_f( &
     context, walk_num, elec_num, cord_num, rescale_factor_ee,  &
     ee_distance, een_rescaled_e) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer*8             , intent(in)  :: walk_num
  integer*8             , intent(in)  :: elec_num
  integer*8             , intent(in)  :: cord_num
  double precision      , intent(in)  :: rescale_factor_ee
  double precision      , intent(in)  :: ee_distance(elec_num,elec_num,walk_num)
  double precision      , intent(out) :: een_rescaled_e(elec_num,elec_num,0:cord_num,walk_num)
  double precision,dimension(:,:),allocatable :: een_rescaled_e_ij
  double precision                    :: x
  integer*8                           :: i, j, k, l, nw

  allocate(een_rescaled_e_ij(elec_num * (elec_num - 1) / 2, cord_num + 1))

  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 (cord_num < 0) then
     info = QMCKL_INVALID_ARG_4
     return
  endif

  ! Prepare table of exponentiated distances raised to appropriate power
  een_rescaled_e             = 0.0d0
  do nw = 1, walk_num
  een_rescaled_e_ij       = 0.0d0
  een_rescaled_e_ij(:, 1) = 1.0d0


  k = 0
  do j = 1, elec_num
    do i = 1, j - 1
      k = k + 1
      een_rescaled_e_ij(k, 2) = dexp(-rescale_factor_ee * ee_distance(i, j, nw))
    end do
  end do


  do l = 2, cord_num
    do k = 1, elec_num * (elec_num - 1)/2
      een_rescaled_e_ij(k, l + 1) = een_rescaled_e_ij(k, l + 1 - 1) * een_rescaled_e_ij(k, 2)
    end do
  end do

  ! prepare the actual een table
  een_rescaled_e(:, :, 0, nw) = 1.0d0

  do l = 1, cord_num
    k = 0
    do j = 1, elec_num
      do i = 1, j - 1
        k = k + 1
        x = een_rescaled_e_ij(k, l + 1)
        een_rescaled_e(i, j, l, nw) = x
        een_rescaled_e(j, i, l, nw) = x
      end do
    end do
  end do

  do l = 0, cord_num
    do j = 1, elec_num
      een_rescaled_e(j, j, l, nw) = 0.0d0
    end do
  end do

  end do

end function qmckl_compute_een_rescaled_e_doc_f

qmckl_exit_code qmckl_compute_een_rescaled_e_hpc (
          const qmckl_context context,
          const int64_t walk_num,
          const int64_t elec_num,
          const int64_t cord_num,
          const double rescale_factor_ee,
          const double* ee_distance,
          double* const een_rescaled_e ) {

  double   *een_rescaled_e_ij;
  double   x;
  const int64_t   elec_pairs = (elec_num * (elec_num - 1)) / 2;
  const int64_t   len_een_ij = elec_pairs * (cord_num + 1);
  int64_t   k;

  // number of element for the een_rescaled_e_ij[N_e*(N_e-1)/2][cord+1]
  // probably in C is better [cord+1, Ne*(Ne-1)/2]
  //elec_pairs = (elec_num * (elec_num - 1)) / 2;
  //len_een_ij = elec_pairs * (cord_num + 1);
  een_rescaled_e_ij = (double *) malloc (len_een_ij * sizeof(double));

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

  if (walk_num <= 0) {
     return QMCKL_INVALID_ARG_2;
  }

  if (elec_num <= 0) {
     return QMCKL_INVALID_ARG_3;
  }

  if (cord_num < 0) {
     return QMCKL_INVALID_ARG_4;
  }

  // Prepare table of exponentiated distances raised to appropriate power
  // init

  for (int kk = 0; kk < walk_num*(cord_num+1)*elec_num*elec_num; ++kk) {
    een_rescaled_e[kk]= 0.0;
  }

  /*
  for (int nw = 0; nw < walk_num; ++nw) {
    for (int l = 0; l < (cord_num + 1); ++l) {
      for (int i = 0; i < elec_num; ++i) {
        for (int j = 0; j < elec_num; ++j) {
          een_rescaled_e[j + i*elec_num + l*elec_num*elec_num + nw*(cord_num+1)*elec_num*elec_num]= 0.0;
        }
      }
    }
  }
  */

  for (int nw = 0; nw < walk_num; ++nw) {

    for (int kk = 0; kk < len_een_ij; ++kk) {
      // this array initialized at 0 except een_rescaled_e_ij(:, 1) = 1.0d0
      // and the arrangement of indices is [cord_num+1, ne*(ne-1)/2]
      een_rescaled_e_ij[kk]= ( kk < (elec_pairs) ? 1.0 : 0.0 );
    }

    k = 0;
    for (int i = 0; i < elec_num; ++i) {
      for (int j = 0; j < i; ++j) {
        // een_rescaled_e_ij(k, 2) = dexp(-rescale_factor_ee * ee_distance(i, j, nw));
        een_rescaled_e_ij[k + elec_pairs] = exp(-rescale_factor_ee * \
                                    ee_distance[j + i*elec_num + nw*(elec_num*elec_num)]);
        k = k + 1;
      }
    }


    for (int l = 2; l < (cord_num+1); ++l) {
      for (int k = 0; k < elec_pairs; ++k) {
      // een_rescaled_e_ij(k, l + 1) = een_rescaled_e_ij(k, l + 1 - 1) * een_rescaled_e_ij(k, 2)
        een_rescaled_e_ij[k+l*elec_pairs] = een_rescaled_e_ij[k + (l - 1)*elec_pairs] * \
                                                  een_rescaled_e_ij[k + elec_pairs];
      }
    }


  // prepare the actual een table
  for (int i = 0; i < elec_num; ++i){
    for (int j = 0; j < elec_num; ++j) {
      een_rescaled_e[j + i*elec_num + 0 + nw*(cord_num+1)*elec_num*elec_num] = 1.0;
    }
  }

    // Up to here it should work.
  for ( int l = 1; l < (cord_num+1); ++l) {
    k = 0;
    for (int i = 0; i < elec_num; ++i) {
      for (int j = 0; j < i; ++j) {
        x = een_rescaled_e_ij[k + l*elec_pairs];
        een_rescaled_e[j + i*elec_num + l*elec_num*elec_num + nw*elec_num*elec_num*(cord_num+1)] = x;
        een_rescaled_e[i + j*elec_num + l*elec_num*elec_num + nw*elec_num*elec_num*(cord_num+1)] = x;
        k = k + 1;
      }
    }
  }

  for (int l = 0; l < (cord_num + 1); ++l) {
    for (int j = 0; j < elec_num; ++j) {
      een_rescaled_e[j + j*elec_num + l*elec_num*elec_num + nw*elec_num*elec_num*(cord_num+1)] = 0.0;
    }
  }

  }

  free(een_rescaled_e_ij);

  return QMCKL_SUCCESS;
}

qmckl_exit_code qmckl_compute_een_rescaled_e_doc (
      const qmckl_context context,
      const int64_t walk_num,
      const int64_t elec_num,
      const int64_t cord_num,
      const double rescale_factor_ee,
      const double* ee_distance,
      double* const een_rescaled_e );

qmckl_exit_code qmckl_compute_een_rescaled_e_hpc (
      const qmckl_context context,
      const int64_t walk_num,
      const int64_t elec_num,
      const int64_t cord_num,
      const double rescale_factor_ee,
      const double* ee_distance,
      double* const een_rescaled_e );

qmckl_exit_code qmckl_compute_een_rescaled_e (
      const qmckl_context context,
      const int64_t walk_num,
      const int64_t elec_num,
      const int64_t cord_num,
      const double rescale_factor_ee,
      const double* ee_distance,
      double* const een_rescaled_e ) {

#ifdef HAVE_HPC
return qmckl_compute_een_rescaled_e_hpc(context, walk_num, elec_num, cord_num, rescale_factor_ee, ee_distance, een_rescaled_e);
#else
return qmckl_compute_een_rescaled_e_doc(context, walk_num, elec_num, cord_num, rescale_factor_ee, ee_distance, een_rescaled_e);
#endif
}

3.9.3 Test

assert(qmckl_electron_provided(context));


double een_rescaled_e[walk_num][(cord_num + 1)][elec_num][elec_num];
rc = qmckl_get_jastrow_een_rescaled_e(context, &(een_rescaled_e[0][0][0][0]),elec_num*elec_num*(cord_num+1)*walk_num);

// value of (0,2,1)
assert(fabs(een_rescaled_e[0][1][0][2]-0.08084493981483197)   < 1.e-12);
assert(fabs(een_rescaled_e[0][1][0][3]-0.1066745707571846)    < 1.e-12);
assert(fabs(een_rescaled_e[0][1][0][4]-0.01754273169464735)   < 1.e-12);
assert(fabs(een_rescaled_e[0][2][1][3]-0.02214680362033448)   < 1.e-12);
assert(fabs(een_rescaled_e[0][2][1][4]-0.0005700154999202759) < 1.e-12);
assert(fabs(een_rescaled_e[0][2][1][5]-0.3424402276009091)    < 1.e-12);

3.10 Electron-electron-nucleus rescaled distances for each order and derivatives

een_rescaled_e_deriv_e stores the table of the derivatives of the rescaled distances between all pairs of electrons and raised to the power \(p\) defined by cord_num. Here we take its derivatives required for the een jastrow.

TODO: write formulae

3.10.1 Get

qmckl_exit_code
qmckl_get_jastrow_een_rescaled_e_deriv_e(qmckl_context context,
                                         double* const distance_rescaled,
                                         const int64_t size_max);

3.10.2 Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`walk_num`	`int64_t`	in	Number of walkers
`elec_num`	`int64_t`	in	Number of electrons
`cord_num`	`int64_t`	in	Order of polynomials
`rescale_factor_ee`	`double`	in	Factor to rescale ee distances
`coord_ee`	`double[walk_num][3][elec_num]`	in	Electron coordinates
`ee_distance`	`double[walk_num][elec_num][elec_num]`	in	Electron-electron distances
`een_rescaled_e`	`double[walk_num][0:cord_num][elec_num][elec_num]`	in	Electron-electron distances
`een_rescaled_e_deriv_e`	`double[walk_num][0:cord_num][elec_num][4][elec_num]`	out	Electron-electron rescaled distances

integer function qmckl_compute_factor_een_rescaled_e_deriv_e_f( &
     context, walk_num, elec_num, cord_num, rescale_factor_ee,  &
     coord_ee, ee_distance, een_rescaled_e, een_rescaled_e_deriv_e) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer*8             , intent(in)  :: walk_num
  integer*8             , intent(in)  :: elec_num
  integer*8             , intent(in)  :: cord_num
  double precision      , intent(in)  :: rescale_factor_ee
  double precision      , intent(in)  :: coord_ee(elec_num,3,walk_num)
  double precision      , intent(in)  :: ee_distance(elec_num,elec_num,walk_num)
  double precision      , intent(in)  :: een_rescaled_e(elec_num,elec_num,0:cord_num,walk_num)
  double precision      , intent(out) :: een_rescaled_e_deriv_e(elec_num,4,elec_num,0:cord_num,walk_num)
  double precision,dimension(:,:,:),allocatable  :: elec_dist_deriv_e
  double precision                    :: x, rij_inv, kappa_l
  integer*8                           :: i, j, k, l, nw, ii

  allocate(elec_dist_deriv_e(4,elec_num,elec_num))

  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 (cord_num < 0) then
     info = QMCKL_INVALID_ARG_4
     return
  endif

  ! Prepare table of exponentiated distances raised to appropriate power
  een_rescaled_e_deriv_e     = 0.0d0
  do nw = 1, walk_num
    do j = 1, elec_num
      do i = 1, elec_num
        rij_inv = 1.0d0 / ee_distance(i, j, nw)
        do ii = 1, 3
          elec_dist_deriv_e(ii, i, j) = (coord_ee(i, ii, nw) - coord_ee(j, ii, nw)) * rij_inv
        end do
        elec_dist_deriv_e(4, i, j) = 2.0d0 * rij_inv
      end do
      elec_dist_deriv_e(:, j, j) = 0.0d0
    end do

    ! prepare the actual een table
    do l = 1, cord_num
      kappa_l = - dble(l) * rescale_factor_ee
      do j = 1, elec_num
        do i = 1, elec_num
          een_rescaled_e_deriv_e(i, 1, j, l, nw) = kappa_l * elec_dist_deriv_e(1, i, j)
          een_rescaled_e_deriv_e(i, 2, j, l, nw) = kappa_l * elec_dist_deriv_e(2, i, j)
          een_rescaled_e_deriv_e(i, 3, j, l, nw) = kappa_l * elec_dist_deriv_e(3, i, j)
          een_rescaled_e_deriv_e(i, 4, j, l, nw) = kappa_l * elec_dist_deriv_e(4, i, j)

          een_rescaled_e_deriv_e(i, 4, j, l, nw) = een_rescaled_e_deriv_e(i, 4, j, l, nw)              &
                    + een_rescaled_e_deriv_e(i, 1, j, l, nw) * een_rescaled_e_deriv_e(i, 1, j, l, nw)  &
                    + een_rescaled_e_deriv_e(i, 2, j, l, nw) * een_rescaled_e_deriv_e(i, 2, j, l, nw)  &
                    + een_rescaled_e_deriv_e(i, 3, j, l, nw) * een_rescaled_e_deriv_e(i, 3, j, l, nw)

          een_rescaled_e_deriv_e(i, 1, j, l, nw) = een_rescaled_e_deriv_e(i, 1, j, l, nw) *   &
                                                    een_rescaled_e(i, j, l, nw)
          een_rescaled_e_deriv_e(i, 3, j, l, nw) = een_rescaled_e_deriv_e(i, 2, j, l, nw) *   &
                                                    een_rescaled_e(i, j, l, nw)
          een_rescaled_e_deriv_e(i, 3, j, l, nw) = een_rescaled_e_deriv_e(i, 3, j, l, nw) *   &
                                                    een_rescaled_e(i, j, l, nw)
          een_rescaled_e_deriv_e(i, 4, j, l, nw) = een_rescaled_e_deriv_e(i, 4, j, l, nw) *   &
                                                    een_rescaled_e(i, j, l, nw)
        end do
      end do
    end do
  end do

end function qmckl_compute_factor_een_rescaled_e_deriv_e_f

3.10.3 Test

double een_rescaled_e_deriv_e[walk_num][(cord_num + 1)][elec_num][4][elec_num];
size_max=walk_num*(cord_num + 1)*elec_num*4*elec_num;
rc = qmckl_get_jastrow_een_rescaled_e_deriv_e(context,
          &(een_rescaled_e_deriv_e[0][0][0][0][0]),size_max);

// value of (0,0,0,2,1)
assert(fabs(een_rescaled_e_deriv_e[0][1][0][0][2] + 0.05991352796887283   ) < 1.e-12);
assert(fabs(een_rescaled_e_deriv_e[0][1][0][0][3] + 0.011714035071545248  ) < 1.e-12);
assert(fabs(een_rescaled_e_deriv_e[0][1][0][0][4] + 0.00441398875758468   ) < 1.e-12);
assert(fabs(een_rescaled_e_deriv_e[0][2][1][0][3] + 0.013553180060167595  ) < 1.e-12);
assert(fabs(een_rescaled_e_deriv_e[0][2][1][0][4] + 0.00041342909359870457) < 1.e-12);
assert(fabs(een_rescaled_e_deriv_e[0][2][1][0][5] + 0.5880599146214673    ) < 1.e-12);

3.11 Electron-nucleus rescaled distances

en_distance_rescaled stores the matrix of the rescaled distances between electrons and nuclei.

\[ C_{ij} = \left( 1 - \exp{-\kappa C_{ij}}\right)/\kappa \]

where \(C_{ij}\) is the matrix of electron-nucleus distances.

3.11.1 Get

qmckl_exit_code qmckl_get_electron_en_distance_rescaled(qmckl_context context, double* distance_rescaled);

3.11.2 Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`elec_num`	`int64_t`	in	Number of electrons
`nucl_num`	`int64_t`	in	Number of nuclei
`type_nucl_num`	`int64_t`	in	Number of types of nuclei
`type_nucl_vector`	`int64_t[nucl_num]`	in	Number of types of nuclei
`rescale_factor_en`	`double[type_nucl_num]`	in	The factor for rescaled distances
`walk_num`	`int64_t`	in	Number of walkers
`elec_coord`	`double[3][walk_num][elec_num]`	in	Electron coordinates
`nucl_coord`	`double[3][elec_num]`	in	Nuclear coordinates
`en_distance_rescaled`	`double[walk_num][nucl_num][elec_num]`	out	Electron-nucleus distances

integer function qmckl_compute_en_distance_rescaled_f(context, elec_num, nucl_num, type_nucl_num, &
     type_nucl_vector, rescale_factor_en, walk_num, elec_coord, &
     nucl_coord, en_distance_rescaled) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer*8             , intent(in)  :: elec_num
  integer*8             , intent(in)  :: nucl_num
  integer*8             , intent(in)  :: type_nucl_num
  integer*8             , intent(in)  :: type_nucl_vector(nucl_num)
  double precision      , intent(in)  :: rescale_factor_en(type_nucl_num)
  integer*8             , intent(in)  :: walk_num
  double precision      , intent(in)  :: elec_coord(elec_num,walk_num,3)
  double precision      , intent(in)  :: nucl_coord(nucl_num,3)
  double precision      , intent(out) :: en_distance_rescaled(elec_num,nucl_num,walk_num)

  integer*8 :: i, k
  double precision      :: coord(3)

  info = QMCKL_SUCCESS

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

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

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

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

  do i=1, nucl_num
     coord(1:3) = nucl_coord(i,1:3)
     do k=1,walk_num
        info = qmckl_distance_rescaled(context, 'T', 'T', elec_num, 1_8, &
             elec_coord(1,k,1), elec_num*walk_num, coord, 1_8, &
             en_distance_rescaled(1,i,k), elec_num, rescale_factor_en(type_nucl_vector(i)))
        if (info /= QMCKL_SUCCESS) then
           return
        endif
     end do
  end do

end function qmckl_compute_en_distance_rescaled_f

3.11.3 Test

assert(qmckl_electron_provided(context));
assert(qmckl_nucleus_provided(context));

double en_distance_rescaled[walk_num][nucl_num][elec_num];

rc = qmckl_check(context,
                 qmckl_get_electron_en_distance_rescaled(context, &(en_distance_rescaled[0][0][0]))
                 );
assert (rc == QMCKL_SUCCESS);

// (e,n,w) in Fortran notation
// (1,1,1)
assert(fabs(en_distance_rescaled[0][0][0] - 0.4435709484118112) < 1.e-12);

// (1,2,1)
assert(fabs(en_distance_rescaled[0][1][0] - 0.8993601506374442) < 1.e-12);

// (2,1,1)
assert(fabs(en_distance_rescaled[0][0][1] - 0.46760219699910477) < 1.e-12);

// (1,1,2)
assert(fabs(en_distance_rescaled[0][0][5] - 0.1875631834682101) < 1.e-12);

// (1,2,2)
assert(fabs(en_distance_rescaled[0][1][5] - 0.8840716589810682) < 1.e-12);

// (2,1,2)
assert(fabs(en_distance_rescaled[0][0][6] - 0.42640469987268914) < 1.e-12);

3.12 Electron-electron-nucleus rescaled distance gradients and laplacian with respect to electron coords

The rescaled distances which is given as \(R = (1 - \exp{-\kappa r})/\kappa\) needs to be perturbed with respect to the nuclear coordinates. This data is stored in the en_distance_rescaled_deriv_e tensor. The The first three elements of this three index tensor [4][nucl_num][elec_num] gives the derivatives in the x, y, and z directions \(dx, dy, dz\) and the last index gives the Laplacian \(\partial x^2 + \partial y^2 + \partial z^2\).

3.12.1 Get

qmckl_exit_code qmckl_get_electron_en_distance_rescaled_deriv_e(qmckl_context context, double* distance_rescaled_deriv_e);

3.12.2 Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`elec_num`	`int64_t`	in	Number of electrons
`nucl_num`	`int64_t`	in	Number of nuclei
`type_nucl_num`	`int64_t`	in	Number of nucleus types
`type_nucl_vector`	`int64_t[nucl_num]`	in	Array of nucleus types
`rescale_factor_en`	`double[nucl_num]`	in	The factors for rescaled distances
`walk_num`	`int64_t`	in	Number of walkers
`elec_coord`	`double[3][walk_num][elec_num]`	in	Electron coordinates
`nucl_coord`	`double[3][elec_num]`	in	Nuclear coordinates
`en_distance_rescaled_deriv_e`	`double[walk_num][nucl_num][elec_num][4]`	out	Electron-nucleus distance derivatives

integer function qmckl_compute_en_distance_rescaled_deriv_e_f(context, elec_num, nucl_num, &
     type_nucl_num, type_nucl_vector, rescale_factor_en, walk_num, elec_coord, &
     nucl_coord, en_distance_rescaled_deriv_e) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer*8             , intent(in)  :: elec_num
  integer*8             , intent(in)  :: nucl_num
  integer*8             , intent(in)  :: type_nucl_num
  integer*8             , intent(in)  :: type_nucl_vector(nucl_num)
  double precision      , intent(in)  :: rescale_factor_en(nucl_num)
  integer*8             , intent(in)  :: walk_num
  double precision      , intent(in)  :: elec_coord(elec_num,walk_num,3)
  double precision      , intent(in)  :: nucl_coord(nucl_num,3)
  double precision      , intent(out) :: en_distance_rescaled_deriv_e(4,elec_num,nucl_num,walk_num)

  integer*8 :: i, k
  double precision :: coord(3)

  info = QMCKL_SUCCESS

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

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

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

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

  do i=1, nucl_num
     coord(1:3) = nucl_coord(i,1:3)
     do k=1,walk_num
        info = qmckl_distance_rescaled_deriv_e(context, 'T', 'T', elec_num, 1_8, &
             elec_coord(1,k,1), elec_num*walk_num, coord, 1_8, &
             en_distance_rescaled_deriv_e(1,1,i,k), elec_num, rescale_factor_en(type_nucl_vector(i)))
        if (info /= QMCKL_SUCCESS) then
           return
        endif
     end do
  end do

end function qmckl_compute_en_distance_rescaled_deriv_e_f

3.12.3 Test

assert(qmckl_electron_provided(context));

assert(qmckl_nucleus_provided(context));

double en_distance_rescaled_deriv_e[walk_num][4][nucl_num][elec_num];

rc = qmckl_check(context,
                 qmckl_get_electron_en_distance_rescaled_deriv_e(context, &(en_distance_rescaled_deriv_e[0][0][0][0]))
                 );
assert (rc == QMCKL_SUCCESS);

// TODO: check exact values
//// (e,n,w) in Fortran notation
//// (1,1,1)
//assert(fabs(en_distance_rescaled[0][0][0] - 7.546738741619978) < 1.e-12);
//
//// (1,2,1)
//assert(fabs(en_distance_rescaled[0][1][0] - 8.77102435246984) < 1.e-12);
//
//// (2,1,1)
//assert(fabs(en_distance_rescaled[0][0][1] - 3.698922010513608) < 1.e-12);
//
//// (1,1,2)
//assert(fabs(en_distance_rescaled[1][0][0] - 5.824059436060509) < 1.e-12);
//
//// (1,2,2)
//assert(fabs(en_distance_rescaled[1][1][0] - 7.080482110317645) < 1.e-12);
//
//// (2,1,2)
//assert(fabs(en_distance_rescaled[1][0][1] - 3.1804527583077356) < 1.e-12);

3.13 Electron-electron-nucleus rescaled distances for each order

een_rescaled_n stores the table of the rescaled distances between electrons and nucleii raised to the power \(p\) defined by cord_num:

\[ C_{ia,p} = \left( 1 - \exp{-\kappa C_{ia}} \right)^p \]

where \(C_{ia}\) is the matrix of electron-nucleus distances.

3.13.1 Get

qmckl_exit_code
qmckl_get_jastrow_een_rescaled_n(qmckl_context context,
                                 double* const distance_rescaled,
                                 const int64_t size_max);

3.13.2 Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`walk_num`	`int64_t`	in	Number of walkers
`elec_num`	`int64_t`	in	Number of electrons
`nucl_num`	`int64_t`	in	Number of atoms
`type_nucl_num`	`int64_t`	in	Number of atom types
`type_nucl_vector`	`int64_t[nucl_num]`	in	Types of atoms
`cord_num`	`int64_t`	in	Order of polynomials
`rescale_factor_en`	`double[nucl_num]`	in	Factor to rescale ee distances
`en_distance`	`double[walk_num][elec_num][nucl_num]`	in	Electron-nucleus distances
`een_rescaled_n`	`double[walk_num][0:cord_num][nucl_num][elec_num]`	out	Electron-nucleus rescaled distances

integer function qmckl_compute_een_rescaled_n_f( &
     context, walk_num, elec_num, nucl_num, &
     type_nucl_num, type_nucl_vector, cord_num, rescale_factor_en,  &
     en_distance, een_rescaled_n) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer*8             , intent(in)  :: walk_num
  integer*8             , intent(in)  :: elec_num
  integer*8             , intent(in)  :: nucl_num
  integer*8             , intent(in)  :: type_nucl_num
  integer*8             , intent(in)  :: type_nucl_vector(nucl_num)
  integer*8             , intent(in)  :: cord_num
  double precision      , intent(in)  :: rescale_factor_en(type_nucl_num)
  double precision      , intent(in)  :: en_distance(elec_num,nucl_num,walk_num)
  double precision      , intent(out) :: een_rescaled_n(elec_num,nucl_num,0:cord_num,walk_num)
  double precision                    :: x
  integer*8                           :: i, a, k, l, nw

  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 (cord_num < 0) then
     info = QMCKL_INVALID_ARG_5
     return
  endif

  ! Prepare table of exponentiated distances raised to appropriate power
  een_rescaled_n             = 0.0d0
  do nw = 1, walk_num

     ! prepare the actual een table
     een_rescaled_n(:, :, 0, nw) = 1.0d0

     do a = 1, nucl_num
        do i = 1, elec_num
           een_rescaled_n(i, a, 1, nw) = dexp(-rescale_factor_en(type_nucl_vector(a)) * en_distance(i, a, nw))
        end do
     end do

     do l = 2, cord_num
        do a = 1, nucl_num
           do i = 1, elec_num
              een_rescaled_n(i, a, l, nw) = een_rescaled_n(i, a, l - 1, nw) * een_rescaled_n(i, a, 1, nw)
           end do
        end do
     end do

  end do

end function qmckl_compute_een_rescaled_n_f

/*
qmckl_exit_code qmckl_compute_een_rescaled_n (
      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,
      int64_t* const type_nucl_vector,
      const int64_t cord_num,
      const double* rescale_factor_en,
      const double* en_distance,
      double* const een_rescaled_n ) {


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

  if (walk_num <= 0) {
    return QMCKL_INVALID_ARG_2;
  }

  if (elec_num <= 0) {
    return QMCKL_INVALID_ARG_3;
  }

  if (nucl_num <= 0) {
    return QMCKL_INVALID_ARG_4;
  }

  if (cord_num < 0) {
    return QMCKL_INVALID_ARG_5;
  }

  // Prepare table of exponentiated distances raised to appropriate power
  for (int i = 0; i < (walk_num*(cord_num+1)*nucl_num*elec_num); ++i) {
    een_rescaled_n[i] = 1.0;
  }

  for (int nw = 0; nw < walk_num; ++nw) {
    for (int a = 0; a < nucl_num; ++a) {
      for (int i = 0; i < elec_num; ++i) {
        een_rescaled_n[i + a*elec_num + nw * elec_num*nucl_num*(cord_num+1)] = 1.0;
        een_rescaled_n[i + a*elec_num + elec_num*nucl_num + nw*elec_num*nucl_num*(cord_num+1)] =
          exp(-rescale_factor_en[type_nucl_vector[a]] * en_distance[i + a*elec_num + nw*elec_num*nucl_num]);
      }
    }

    for (int l = 2; l < (cord_num+1); ++l){
      for (int a = 0; a < nucl_num; ++a) {
        for (int i = 0; i < elec_num; ++i) {
          een_rescaled_n[i + a*elec_num + l*elec_num*nucl_num + nw*elec_num*nucl_num*(cord_num+1)] =
            een_rescaled_n[i + a*elec_num + (l-1)*elec_num*nucl_num + nw*elec_num*nucl_num*(cord_num+1)] * 
            een_rescaled_n[i + a*elec_num +       elec_num*nucl_num + nw*elec_num*nucl_num*(cord_num+1)];
        }
      }
    }

  }

  return QMCKL_SUCCESS;
}
*/

3.13.3 Test

assert(qmckl_electron_provided(context));

double een_rescaled_n[walk_num][(cord_num + 1)][nucl_num][elec_num];
size_max=walk_num*(cord_num + 1)*nucl_num*elec_num;
rc = qmckl_get_jastrow_een_rescaled_n(context, &(een_rescaled_n[0][0][0][0]),size_max);

// value of (0,2,1)
assert(fabs(een_rescaled_n[0][1][0][2]-0.10612983920006765)  < 1.e-12);
assert(fabs(een_rescaled_n[0][1][0][3]-0.135652809635553)    < 1.e-12);
assert(fabs(een_rescaled_n[0][1][0][4]-0.023391817607642338) < 1.e-12);
assert(fabs(een_rescaled_n[0][2][1][3]-0.880957224822116)    < 1.e-12);
assert(fabs(een_rescaled_n[0][2][1][4]-0.027185942659395074) < 1.e-12);
assert(fabs(een_rescaled_n[0][2][1][5]-0.01343938025140174)  < 1.e-12);

3.14 Electron-nucleus rescaled distances for each order and derivatives

een_rescaled_n_deriv_e stores the table of the rescaled distances between electrons and nucleii raised to the power \(p\) defined by cord_num:

3.14.1 Get

qmckl_exit_code
qmckl_get_jastrow_een_rescaled_n_deriv_e(qmckl_context context,
                                         double* const distance_rescaled,
                                         const int64_t size_max);

3.14.2 Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`walk_num`	`int64_t`	in	Number of walkers
`elec_num`	`int64_t`	in	Number of electrons
`nucl_num`	`int64_t`	in	Number of atoms
`type_nucl_num`	`int64_t`	in	Number of atom types
`type_nucl_vector`	`int64_t[nucl_num]`	in	Types of atoms
`cord_num`	`int64_t`	in	Order of polynomials
`rescale_factor_en`	`double[nucl_num]`	in	Factor to rescale ee distances
`coord_ee`	`double[walk_num][3][elec_num]`	in	Electron coordinates
`coord_en`	`double[3][nucl_num]`	in	Nuclear coordinates
`en_distance`	`double[walk_num][elec_num][nucl_num]`	in	Electron-nucleus distances
`een_rescaled_n`	`double[walk_num][0:cord_num][nucl_num][elec_num]`	in	Electron-nucleus distances
`een_rescaled_n_deriv_e`	`double[walk_num][0:cord_num][nucl_num][4][elec_num]`	out	Electron-nucleus rescaled distances

integer function qmckl_compute_factor_een_rescaled_n_deriv_e_f( &
     context, walk_num, elec_num, nucl_num, type_nucl_num, type_nucl_vector, &
     cord_num, rescale_factor_en, &
     coord_ee, coord_en, en_distance, een_rescaled_n, een_rescaled_n_deriv_e) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer*8             , intent(in)  :: walk_num
  integer*8             , intent(in)  :: elec_num
  integer*8             , intent(in)  :: nucl_num
  integer*8             , intent(in)  :: type_nucl_num
  integer*8             , intent(in)  :: type_nucl_vector(nucl_num)
  integer*8             , intent(in)  :: cord_num
  double precision      , intent(in)  :: rescale_factor_en(type_nucl_num)
  double precision      , intent(in)  :: coord_ee(elec_num,3,walk_num)
  double precision      , intent(in)  :: coord_en(nucl_num,3)
  double precision      , intent(in)  :: en_distance(elec_num,nucl_num,walk_num)
  double precision      , intent(in)  :: een_rescaled_n(elec_num,nucl_num,0:cord_num,walk_num)
  double precision      , intent(out) :: een_rescaled_n_deriv_e(elec_num,4,nucl_num,0:cord_num,walk_num)
  double precision,dimension(:,:,:),allocatable :: elnuc_dist_deriv_e
  double precision                    :: x, ria_inv, kappa_l
  integer*8                           :: i, a, k, l, nw, ii

  allocate(elnuc_dist_deriv_e(4, elec_num, nucl_num))

  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 (cord_num < 0) then
     info = QMCKL_INVALID_ARG_5
     return
  endif

  ! Prepare table of exponentiated distances raised to appropriate power
  een_rescaled_n_deriv_e             = 0.0d0
  do nw = 1, walk_num

  ! prepare the actual een table
  do a = 1, nucl_num
    do i = 1, elec_num
      ria_inv = 1.0d0 / en_distance(i, a, nw)
      do ii = 1, 3
        elnuc_dist_deriv_e(ii, i, a) = (coord_ee(i, ii, nw) - coord_en(a, ii)) * ria_inv
      end do
      elnuc_dist_deriv_e(4, i, a) = 2.0d0 * ria_inv
    end do
  end do

  do l = 0, cord_num
    do a = 1, nucl_num
      kappa_l = - dble(l) * rescale_factor_en(type_nucl_vector(a))
      do i = 1, elec_num
        een_rescaled_n_deriv_e(i, 1, a, l, nw) = kappa_l * elnuc_dist_deriv_e(1, i, a)
        een_rescaled_n_deriv_e(i, 2, a, l, nw) = kappa_l * elnuc_dist_deriv_e(2, i, a)
        een_rescaled_n_deriv_e(i, 3, a, l, nw) = kappa_l * elnuc_dist_deriv_e(3, i, a)
        een_rescaled_n_deriv_e(i, 4, a, l, nw) = kappa_l * elnuc_dist_deriv_e(4, i, a)

        een_rescaled_n_deriv_e(i, 4, a, l, nw) = een_rescaled_n_deriv_e(i, 4, a, l, nw)           &
                + een_rescaled_n_deriv_e(i, 1, a, l, nw) * een_rescaled_n_deriv_e(i, 1, a, l, nw) &
                + een_rescaled_n_deriv_e(i, 2, a, l, nw) * een_rescaled_n_deriv_e(i, 2, a, l, nw) &
                + een_rescaled_n_deriv_e(i, 3, a, l, nw) * een_rescaled_n_deriv_e(i, 3, a, l, nw)

        een_rescaled_n_deriv_e(i, 1, a, l, nw) = een_rescaled_n_deriv_e(i, 1, a, l, nw) * &
                                                  een_rescaled_n(i, a, l, nw)
        een_rescaled_n_deriv_e(i, 2, a, l, nw) = een_rescaled_n_deriv_e(i, 2, a, l, nw) * &
                                                  een_rescaled_n(i, a, l, nw)
        een_rescaled_n_deriv_e(i, 3, a, l, nw) = een_rescaled_n_deriv_e(i, 3, a, l, nw) * &
                                                  een_rescaled_n(i, a, l, nw)
        een_rescaled_n_deriv_e(i, 4, a, l, nw) = een_rescaled_n_deriv_e(i, 4, a, l, nw) * &
                                                  een_rescaled_n(i, a, l, nw)
      end do
    end do
  end do
  end do

end function qmckl_compute_factor_een_rescaled_n_deriv_e_f

3.14.3 Test

assert(qmckl_electron_provided(context));

double een_rescaled_n_deriv_e[walk_num][(cord_num + 1)][nucl_num][4][elec_num];
size_max=walk_num*(cord_num + 1)*nucl_num*4*elec_num;
rc = qmckl_get_jastrow_een_rescaled_n_deriv_e(context, &(een_rescaled_n_deriv_e[0][0][0][0][0]),size_max);

// value of (0,2,1)
assert(fabs(een_rescaled_n_deriv_e[0][1][0][0][2]+0.07633444246999128   )  < 1.e-12);
assert(fabs(een_rescaled_n_deriv_e[0][1][0][0][3]-0.00033282346259738276)  < 1.e-12);
assert(fabs(een_rescaled_n_deriv_e[0][1][0][0][4]+0.004775370547333061  )  < 1.e-12);
assert(fabs(een_rescaled_n_deriv_e[0][2][1][0][3]-0.1362654644223866    )  < 1.e-12);
assert(fabs(een_rescaled_n_deriv_e[0][2][1][0][4]+0.0231253431662794    )  < 1.e-12);
assert(fabs(een_rescaled_n_deriv_e[0][2][1][0][5]-0.001593334817691633  )  < 1.e-12);

3.15 Prepare for electron-electron-nucleus Jastrow \(f_{een}\)

Prepare c_vector_full and lkpm_combined_index tables required for the calculation of the three-body jastrow factor_een and its derivative factor_een_deriv_e.

3.15.1 Get

qmckl_exit_code qmckl_get_jastrow_dim_c_vector(qmckl_context context, int64_t* const dim_c_vector);
qmckl_exit_code qmckl_get_jastrow_c_vector_full(qmckl_context context, double* const c_vector_full);
qmckl_exit_code qmckl_get_jastrow_lkpm_combined_index(qmckl_context context, int64_t* const lkpm_combined_index);
qmckl_exit_code qmckl_get_jastrow_tmp_c(qmckl_context context, double* const tmp_c);
qmckl_exit_code qmckl_get_jastrow_dtmp_c(qmckl_context context, double* const dtmp_c);

3.15.2 Compute dim_c_vector

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`cord_num`	`int64_t`	in	Order of polynomials
`dim_c_vector`	`int64_t`	out	dimension of c_vector_full table

integer function qmckl_compute_dim_c_vector_f( &
     context, cord_num, dim_c_vector) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer*8             , intent(in)  :: cord_num
  integer*8             , intent(out) :: dim_c_vector
  double precision                    :: x
  integer*8                           :: i, a, k, l, p, lmax

  info = QMCKL_SUCCESS

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

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

  dim_c_vector = 0

  do p = 2, cord_num
    do k = p - 1, 0, -1
      if (k .ne. 0) then
        lmax = p - k
      else
        lmax = p - k - 2
      endif
      do l = lmax, 0, -1
        if (iand(p - k - l, 1_8) == 1) cycle
        dim_c_vector = dim_c_vector + 1
      end do
    end do
  end do

end function qmckl_compute_dim_c_vector_f

qmckl_exit_code qmckl_compute_dim_c_vector (
      const qmckl_context context,
      const int64_t cord_num,
      int64_t* const dim_c_vector){

  int         lmax;


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

  if (cord_num < 0) {
    return QMCKL_INVALID_ARG_2;
  }

  *dim_c_vector = 0;

  for (int p=2; p <= cord_num; ++p){
    for (int k=p-1; k >= 0; --k) {
      if (k != 0) {
        lmax = p - k;
      } else {
        lmax = p - k - 2;
      }
      for (int l = lmax; l >= 0; --l) {
        if ( ((p - k - l) & 1)==1) continue;
        *dim_c_vector=*dim_c_vector+1;
      }
    }
  }

  return QMCKL_SUCCESS;
}

3.15.3 Compute c_vector_full

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`nucl_num`	`int64_t`	in	Number of atoms
`dim_c_vector`	`int64_t`	in	dimension of cord full table
`type_nucl_num`	`int64_t`	in	dimension of cord full table
`type_nucl_vector`	`int64_t[nucl_num]`	in	dimension of cord full table
`c_vector`	`double[dim_c_vector][type_nucl_num]`	in	dimension of cord full table
`c_vector_full`	`double[dim_c_vector][nucl_num]`	out	Full list of coefficients

integer function qmckl_compute_c_vector_full_doc_f( &
     context, nucl_num, dim_c_vector, type_nucl_num,  &
     type_nucl_vector, c_vector, c_vector_full) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer*8             , intent(in)  :: nucl_num
  integer*8             , intent(in)  :: dim_c_vector
  integer*8             , intent(in)  :: type_nucl_num
  integer*8             , intent(in)  :: type_nucl_vector(nucl_num)
  double precision      , intent(in)  :: c_vector(type_nucl_num, dim_c_vector)
  double precision      , intent(out) :: c_vector_full(nucl_num,dim_c_vector)
  double precision                    :: x
  integer*8                           :: i, a, k, l, nw

  info = QMCKL_SUCCESS

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

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

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

  if (dim_c_vector < 0) then
     info = QMCKL_INVALID_ARG_5
     return
  endif


  do a = 1, nucl_num
    c_vector_full(a,1:dim_c_vector) = c_vector(type_nucl_vector(a),1:dim_c_vector)
  end do

end function qmckl_compute_c_vector_full_doc_f

qmckl_exit_code qmckl_compute_c_vector_full_hpc (
          const qmckl_context context,
          const int64_t nucl_num,
          const int64_t dim_c_vector,
          const int64_t type_nucl_num,
          const int64_t* type_nucl_vector,
          const double* c_vector,
          double* const c_vector_full ) {

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

  if (nucl_num <= 0) {
     return QMCKL_INVALID_ARG_2;
  }

  if (type_nucl_num <= 0) {
     return QMCKL_INVALID_ARG_4;
  }

  if (dim_c_vector < 0) {
     return QMCKL_INVALID_ARG_5;
  }

  for (int i=0; i < dim_c_vector; ++i) {
    for (int a=0; a < nucl_num; ++a){
      c_vector_full[a + i*nucl_num] = c_vector[(type_nucl_vector[a]-1)+i*type_nucl_num];
    }
  }

  return QMCKL_SUCCESS;
  }

qmckl_exit_code qmckl_compute_c_vector_full_doc (
      const qmckl_context context,
      const int64_t nucl_num,
      const int64_t dim_c_vector,
      const int64_t type_nucl_num,
      const int64_t* type_nucl_vector,
      const double* c_vector,
      double* const c_vector_full );

qmckl_exit_code qmckl_compute_c_vector_full_hpc (
      const qmckl_context context,
      const int64_t nucl_num,
      const int64_t dim_c_vector,
      const int64_t type_nucl_num,
      const int64_t* type_nucl_vector,
      const double* c_vector,
      double* const c_vector_full );

qmckl_exit_code qmckl_compute_c_vector_full (
          const qmckl_context context,
          const int64_t nucl_num,
          const int64_t dim_c_vector,
          const int64_t type_nucl_num,
          const int64_t* type_nucl_vector,
          const double* c_vector,
          double* const c_vector_full ) {

    #ifdef HAVE_HPC
      return qmckl_compute_c_vector_full_hpc(context, nucl_num, dim_c_vector, type_nucl_num, type_nucl_vector, c_vector, c_vector_full);
    #else
      return qmckl_compute_c_vector_full_doc(context, nucl_num, dim_c_vector, type_nucl_num, type_nucl_vector, c_vector, c_vector_full);
    #endif
    }

3.15.4 Compute lkpm_combined_index

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`cord_num`	`int64_t`	in	Order of polynomials
`dim_c_vector`	`int64_t`	in	dimension of cord full table
`lkpm_combined_index`	`int64_t[4][dim_c_vector]`	out	Full list of combined indices

integer function qmckl_compute_lkpm_combined_index_f( &
     context, cord_num, dim_c_vector,  lkpm_combined_index) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer*8             , intent(in)  :: cord_num
  integer*8             , intent(in)  :: dim_c_vector
  integer*8             , intent(out) :: lkpm_combined_index(dim_c_vector, 4)
  double precision                    :: x
  integer*8                           :: i, a, k, l, kk, p, lmax, m

  info = QMCKL_SUCCESS

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

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

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


  kk = 0
  do p = 2, cord_num
    do k = p - 1, 0, -1
      if (k .ne. 0) then
        lmax = p - k
      else
        lmax = p - k - 2
      end if
      do l = lmax, 0, -1
        if (iand(p - k - l, 1_8) .eq. 1) cycle
        m = (p - k - l)/2
        kk = kk + 1
        lkpm_combined_index(kk, 1) = l
        lkpm_combined_index(kk, 2) = k
        lkpm_combined_index(kk, 3) = p
        lkpm_combined_index(kk, 4) = m
      end do
    end do
  end do

end function qmckl_compute_lkpm_combined_index_f

qmckl_exit_code qmckl_compute_lkpm_combined_index (
      const qmckl_context context,
      const int64_t cord_num,
      const int64_t dim_c_vector,
      int64_t* const lkpm_combined_index ) {

  int kk, lmax, m;

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

  if (cord_num < 0) {
     return QMCKL_INVALID_ARG_2;
  }

  if (dim_c_vector < 0) {
     return QMCKL_INVALID_ARG_3;
  }

/*
*/
  kk = 0;
  for (int p = 2; p <= cord_num; ++p) {
    for (int k=(p-1); k >= 0; --k) {
      if (k != 0) {
        lmax = p - k;
      } else {
        lmax = p - k - 2;
      }
      for (int l=lmax; l >= 0; --l) {
        if (((p - k - l) & 1) == 1) continue;
        m = (p - k - l)/2;
        lkpm_combined_index[kk                  ] = l;
        lkpm_combined_index[kk +   dim_c_vector] = k;
        lkpm_combined_index[kk + 2*dim_c_vector] = p;
        lkpm_combined_index[kk + 3*dim_c_vector] = m;
        kk = kk + 1;
      }
    }
  }

  return QMCKL_SUCCESS;
}

3.15.5 Compute tmp_c

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`cord_num`	`int64_t`	in	Order of polynomials
`elec_num`	`int64_t`	in	Number of electrons
`nucl_num`	`int64_t`	in	Number of nucleii
`walk_num`	`int64_t`	in	Number of walkers
`een_rescaled_e`	`double[walk_num][0:cord_num][elec_num][elec_num]`	in	Electron-electron rescaled factor
`een_rescaled_n`	`double[walk_num][0:cord_num][nucl_num][elec_num]`	in	Electron-nucleus rescaled factor
`tmp_c`	`double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]`	out	vector of non-zero coefficients

qmckl_exit_code qmckl_compute_tmp_c (const qmckl_context context,
                                     const int64_t cord_num,
                                     const int64_t elec_num,
                                     const int64_t nucl_num,
                                     const int64_t walk_num,
                                     const double* een_rescaled_e,
                                     const double* een_rescaled_n,
                                     double* const tmp_c )
{
#ifdef HAVE_HPC
  return qmckl_compute_tmp_c_hpc(context, cord_num, elec_num, nucl_num, walk_num, een_rescaled_e, een_rescaled_n, tmp_c);
#else
  return qmckl_compute_tmp_c_doc(context, cord_num, elec_num, nucl_num, walk_num, een_rescaled_e, een_rescaled_n, tmp_c);
#endif
}

integer function qmckl_compute_tmp_c_doc_f( &
     context, cord_num, elec_num, nucl_num, &
     walk_num, een_rescaled_e, een_rescaled_n, tmp_c) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer*8             , intent(in)  :: cord_num
  integer*8             , intent(in)  :: elec_num
  integer*8             , intent(in)  :: nucl_num
  integer*8             , intent(in)  :: walk_num
  double precision      , intent(in)  :: een_rescaled_e(elec_num, elec_num, 0:cord_num, walk_num)
  double precision      , intent(in)  :: een_rescaled_n(elec_num, nucl_num, 0:cord_num, walk_num)
  double precision      , intent(out) :: tmp_c(elec_num, nucl_num,0:cord_num, 0:cord_num-1, walk_num)
  double precision                    :: x
  integer*8                           :: i, j, a, l, kk, p, lmax, nw
  character                           :: TransA, TransB
  double precision                    :: alpha, beta
  integer*8                           :: M, N, K, LDA, LDB, LDC

  TransA = 'N'
  TransB = 'N'
  alpha = 1.0d0
  beta  = 0.0d0

  info = QMCKL_SUCCESS

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

  if (cord_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

  M = elec_num
  N = nucl_num*(cord_num + 1)
  K = elec_num
  LDA = size(een_rescaled_e,1)
  LDB = size(een_rescaled_n,1)
  LDC = size(tmp_c,1)

  do nw=1, walk_num
  do i=0, cord_num-1
  info = qmckl_dgemm(context, TransA, TransB, M, N, K, alpha,     &
                     een_rescaled_e(1,1,i,nw),LDA*1_8,                     &
                     een_rescaled_n(1,1,0,nw),LDB*1_8,                     &
                     beta,                                       &
                     tmp_c(1,1,0,i,nw),LDC)
  end do
  end do

end function qmckl_compute_tmp_c_doc_f

qmckl_exit_code qmckl_compute_tmp_c_doc (
          const qmckl_context context,
          const int64_t cord_num,
          const int64_t elec_num,
          const int64_t nucl_num,
          const int64_t walk_num,
          const double* een_rescaled_e,
          const double* een_rescaled_n,
          double* const tmp_c );

3.15.6 Compute dtmp_c

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`cord_num`	`int64_t`	in	Order of polynomials
`elec_num`	`int64_t`	in	Number of electrons
`nucl_num`	`int64_t`	in	Number of nucleii
`walk_num`	`int64_t`	in	Number of walkers
`een_rescaled_e_deriv_e`	`double[walk_num][0:cord_num][elec_num][4][elec_num]`	in	Electron-electron rescaled factor derivatives
`een_rescaled_n`	`double[walk_num][0:cord_num][nucl_num][elec_num]`	in	Electron-nucleus rescaled factor
`dtmp_c`	`double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]`	out	vector of non-zero coefficients

qmckl_exit_code
qmckl_compute_dtmp_c (const qmckl_context context,
                      const int64_t cord_num,
                      const int64_t elec_num,
                      const int64_t nucl_num,
                      const int64_t walk_num,
                      const double* een_rescaled_e_deriv_e,
                      const double* een_rescaled_n,
                      double* const dtmp_c )
{
#ifdef HAVE_HPC
  return qmckl_compute_dtmp_c_hpc (context, cord_num, elec_num, nucl_num, walk_num, een_rescaled_e_deriv_e,
                                   een_rescaled_n, dtmp_c );
#else
  return qmckl_compute_dtmp_c_doc (context, cord_num, elec_num, nucl_num, walk_num, een_rescaled_e_deriv_e,
                                   een_rescaled_n, dtmp_c );
#endif
}

integer function qmckl_compute_dtmp_c_doc_f( &
     context, cord_num, elec_num, nucl_num, &
     walk_num, een_rescaled_e_deriv_e, een_rescaled_n, dtmp_c) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer*8             , intent(in)  :: cord_num
  integer*8             , intent(in)  :: elec_num
  integer*8             , intent(in)  :: nucl_num
  integer*8             , intent(in)  :: walk_num
  double precision      , intent(in)  :: een_rescaled_e_deriv_e(elec_num, 4, elec_num, 0:cord_num, walk_num)
  double precision      , intent(in)  :: een_rescaled_n(elec_num, nucl_num, 0:cord_num, walk_num)
  double precision      , intent(out) :: dtmp_c(elec_num, 4, nucl_num,0:cord_num, 0:cord_num-1,  walk_num)
  double precision                    :: x
  integer*8                           :: i, j, a, l, kk, p, lmax, nw, ii
  character                           :: TransA, TransB
  double precision                    :: alpha, beta
  integer*8                           :: M, N, K, LDA, LDB, LDC

  TransA = 'N'
  TransB = 'N'
  alpha = 1.0d0
  beta  = 0.0d0

  info = QMCKL_SUCCESS

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

  if (cord_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

  M = 4*elec_num
  N = nucl_num*(cord_num + 1)
  K = elec_num
  LDA = 4*size(een_rescaled_e_deriv_e,1)
  LDB = size(een_rescaled_n,1)
  LDC = 4*size(dtmp_c,1)

  do nw=1, walk_num
     do i=0, cord_num-1
        info = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha,  &
             een_rescaled_e_deriv_e(1,1,1,i,nw),LDA*1_8,            &
             een_rescaled_n(1,1,0,nw),LDB*1_8,                      &
             beta,                                                  &
             dtmp_c(1,1,1,0,i,nw),LDC)
     end do
  end do

end function qmckl_compute_dtmp_c_doc_f

3.15.7 Test

assert(qmckl_electron_provided(context));

double tmp_c[walk_num][cord_num][cord_num+1][nucl_num][elec_num];
rc = qmckl_get_jastrow_tmp_c(context, &(tmp_c[0][0][0][0][0]));

double dtmp_c[walk_num][cord_num][cord_num+1][nucl_num][4][elec_num];
rc = qmckl_get_jastrow_dtmp_c(context, &(dtmp_c[0][0][0][0][0][0]));

printf("%e\n%e\n", tmp_c[0][0][1][0][0], 2.7083473948352403);
assert(fabs(tmp_c[0][0][1][0][0] - 2.7083473948352403) < 1e-12);

printf("%e\n%e\n", dtmp_c[0][1][0][0][0][0],0.237440520852232);
assert(fabs(dtmp_c[0][1][0][0][0][0] - 0.237440520852232) < 1e-12);

3.16 Electron-electron-nucleus Jastrow \(f_{een}\)

Calculate the electron-electron-nuclear three-body jastrow component factor_een using the above prepared tables.

TODO: write equations.

3.16.1 Get

qmckl_exit_code
qmckl_get_jastrow_factor_een(qmckl_context context,
                             double* const factor_een,
                             const int64_t size_max);

3.16.1.1 Fortran interface

interface
   integer(qmckl_exit_code) function qmckl_get_jastrow_factor_een (context, &
        factor_een, size_max) bind(C)
     use, intrinsic :: iso_c_binding
     import
     implicit none
     integer (qmckl_context) , intent(in), value :: context
     integer(c_int64_t), intent(in), value       :: size_max
     double precision, intent(out)               :: factor_een(size_max)
   end function qmckl_get_jastrow_factor_een
end interface

3.16.2 Compute naive

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`walk_num`	`int64_t`	in	Number of walkers
`elec_num`	`int64_t`	in	Number of electrons
`nucl_num`	`int64_t`	in	Number of nucleii
`cord_num`	`int64_t`	in	order of polynomials
`dim_c_vector`	`int64_t`	in	dimension of full coefficient vector
`c_vector_full`	`double[dim_c_vector][nucl_num]`	in	full coefficient vector
`lkpm_combined_index`	`int64_t[4][dim_c_vector]`	in	combined indices
`een_rescaled_e`	`double[walk_num][elec_num][elec_num][0:cord_num]`	in	Electron-nucleus rescaled
`een_rescaled_n`	`double[walk_num][elec_num][nucl_num][0:cord_num]`	in	Electron-nucleus rescaled factor
`factor_een`	`double[walk_num]`	out	Electron-nucleus jastrow

integer function qmckl_compute_factor_een_naive_f( &
     context, walk_num, elec_num, nucl_num, cord_num,&
     dim_c_vector, c_vector_full, lkpm_combined_index, &
     een_rescaled_e, een_rescaled_n, factor_een) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer*8             , intent(in)  :: walk_num, elec_num, cord_num, nucl_num, dim_c_vector
  integer*8             , intent(in)  :: lkpm_combined_index(dim_c_vector,4)
  double precision      , intent(in)  :: c_vector_full(nucl_num, dim_c_vector)
  double precision      , intent(in)  :: een_rescaled_e(0:cord_num, elec_num, elec_num, walk_num)
  double precision      , intent(in)  :: een_rescaled_n(0:cord_num, nucl_num, elec_num, walk_num)
  double precision      , intent(out) :: factor_een(walk_num)

  integer*8 :: i, a, j, l, k, p, m, n, nw
  double precision :: accu, accu2, cn

  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 (cord_num < 0) then
     info = QMCKL_INVALID_ARG_5
     return
  endif

  factor_een = 0.0d0

  do nw =1, walk_num
  do n = 1, dim_c_vector
    l = lkpm_combined_index(n, 1)
    k = lkpm_combined_index(n, 2)
    p = lkpm_combined_index(n, 3)
    m = lkpm_combined_index(n, 4)

    do a = 1, nucl_num
      accu2 = 0.0d0
      cn = c_vector_full(a, n)
      do j = 1, elec_num
        accu = 0.0d0
        do i = 1, elec_num
          accu = accu + een_rescaled_e(k,i,j,nw) *       &
                        een_rescaled_n(m,a,i,nw)
          !if(nw .eq. 1) then
          !  print *,l,k,p,m,j,i,een_rescaled_e(k,i,j,nw), een_rescaled_n(m,a,i,nw), accu
          !endif
        end do
        accu2 = accu2 + accu * een_rescaled_n(m + l,a,j,nw)
        !print *, l,m,nw,accu, accu2, een_rescaled_n(m + l, a, j, nw), cn, factor_een(nw)
      end do
      factor_een(nw) = factor_een(nw) + accu2 * cn
    end do
  end do
  end do

end function qmckl_compute_factor_een_naive_f

3.16.3 Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`walk_num`	`int64_t`	in	Number of walkers
`elec_num`	`int64_t`	in	Number of electrons
`nucl_num`	`int64_t`	in	Number of nucleii
`cord_num`	`int64_t`	in	order of polynomials
`dim_c_vector`	`int64_t`	in	dimension of full coefficient vector
`c_vector_full`	`double[dim_c_vector][nucl_num]`	in	full coefficient vector
`lkpm_combined_index`	`int64_t[4][dim_c_vector]`	in	combined indices
`tmp_c`	`double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]`	vector of non-zero coefficients
`een_rescaled_n`	`double[walk_num][0:cord_num][nucl_num][elec_num]`	in	Electron-nucleus rescaled factor
`factor_een`	`double[walk_num]`	out	Electron-nucleus jastrow

integer function qmckl_compute_factor_een_f( &
     context, walk_num, elec_num, nucl_num, cord_num,   &
     dim_c_vector, c_vector_full, lkpm_combined_index, &
     tmp_c, een_rescaled_n, factor_een) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer*8             , intent(in)  :: walk_num, elec_num, cord_num, nucl_num, dim_c_vector
  integer*8             , intent(in)  :: lkpm_combined_index(dim_c_vector,4)
  double precision      , intent(in)  :: c_vector_full(nucl_num, dim_c_vector)
  double precision      , intent(in)  :: tmp_c(elec_num, nucl_num,0:cord_num, 0:cord_num-1,  walk_num)
  double precision      , intent(in)  :: een_rescaled_n(elec_num, nucl_num, 0:cord_num, walk_num)
  double precision      , intent(out) :: factor_een(walk_num)

  integer*8 :: i, a, j, l, k, p, m, n, nw
  double precision :: accu, accu2, cn

  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 (cord_num < 0) then
     info = QMCKL_INVALID_ARG_5
     return
  endif

  factor_een = 0.0d0

  do nw =1, walk_num
  do n = 1, dim_c_vector
    l = lkpm_combined_index(n, 1)
    k = lkpm_combined_index(n, 2)
    p = lkpm_combined_index(n, 3)
    m = lkpm_combined_index(n, 4)

    do a = 1, nucl_num
      cn = c_vector_full(a, n)
      if(cn == 0.d0) cycle

      accu = 0.0d0
      do j = 1, elec_num
        accu = accu + een_rescaled_n(j,a,m,nw) * tmp_c(j,a,m+l,k,nw)
      end do
      factor_een(nw) = factor_een(nw) + accu * cn
    end do
  end do
  end do

end function qmckl_compute_factor_een_f

3.16.4 Test

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

double factor_een[walk_num];
rc = qmckl_get_jastrow_factor_een(context, &(factor_een[0]),walk_num);

assert(fabs(factor_een[0] + 0.37407972141304213) < 1e-12);

3.17 Electron-electron-nucleus Jastrow \(f_{een}\) derivative

Calculate the electron-electron-nuclear three-body jastrow component factor_een_deriv_e using the above prepared tables.

TODO: write equations.

3.17.1 Get

qmckl_exit_code
qmckl_get_jastrow_factor_een_deriv_e(qmckl_context context,
                                     double* const factor_een_deriv_e,
                                     const int64_t size_max);

3.17.2 Compute Naive

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`walk_num`	`int64_t`	in	Number of walkers
`elec_num`	`int64_t`	in	Number of electrons
`nucl_num`	`int64_t`	in	Number of nucleii
`cord_num`	`int64_t`	in	order of polynomials
`dim_c_vector`	`int64_t`	in	dimension of full coefficient vector
`c_vector_full`	`double[dim_c_vector][nucl_num]`	in	full coefficient vector
`lkpm_combined_index`	`int64_t[4][dim_c_vector]`	in	combined indices
`een_rescaled_e`	`double[walk_num][elec_num][elec_num][0:cord_num]`	in	Electron-nucleus rescaled
`een_rescaled_n`	`double[walk_num][elec_num][nucl_num][0:cord_num]`	in	Electron-nucleus rescaled factor
`een_rescaled_e_deriv_e`	`double[walk_num][elec_num][4][elec_num][0:cord_num]`	in	Electron-nucleus rescaled
`een_rescaled_n_deriv_e`	`double[walk_num][elec_num][4][nucl_num][0:cord_num]`	in	Electron-nucleus rescaled factor
`factor_een_deriv_e`	`double[walk_num][4][elec_num]`	out	Electron-nucleus jastrow

integer function qmckl_compute_factor_een_deriv_e_naive_f( &
     context, walk_num, elec_num, nucl_num, cord_num, dim_c_vector, &
     c_vector_full, lkpm_combined_index, een_rescaled_e, een_rescaled_n, &
     een_rescaled_e_deriv_e, een_rescaled_n_deriv_e, factor_een_deriv_e)&
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer*8             , intent(in)  :: walk_num, elec_num, cord_num, nucl_num, dim_c_vector
  integer*8             , intent(in)  :: lkpm_combined_index(dim_c_vector, 4)
  double precision      , intent(in)  :: c_vector_full(nucl_num, dim_c_vector)
  double precision      , intent(in)  :: een_rescaled_e(0:cord_num, elec_num, elec_num, walk_num)
  double precision      , intent(in)  :: een_rescaled_n(0:cord_num, nucl_num, elec_num, walk_num)
  double precision      , intent(in)  :: een_rescaled_e_deriv_e(0:cord_num, elec_num, 4, elec_num, walk_num)
  double precision      , intent(in)  :: een_rescaled_n_deriv_e(0:cord_num, nucl_num, 4, elec_num, walk_num)
  double precision      , intent(out) :: factor_een_deriv_e(elec_num, 4, walk_num)

  integer*8 :: i, a, j, l, k, p, m, n, nw
  double precision :: accu, accu2, cn
  double precision :: daccu(1:4), daccu2(1:4)

  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 (cord_num < 0) then
     info = QMCKL_INVALID_ARG_5
     return
  endif

  factor_een_deriv_e = 0.0d0

  do nw =1, walk_num
  do n = 1, dim_c_vector
    l = lkpm_combined_index(n, 1)
    k = lkpm_combined_index(n, 2)
    p = lkpm_combined_index(n, 3)
    m = lkpm_combined_index(n, 4)

    do a = 1, nucl_num
      cn = c_vector_full(a, n)
      do j = 1, elec_num
        accu = 0.0d0
        accu2 = 0.0d0
        daccu = 0.0d0
        daccu2 = 0.0d0
        do i = 1, elec_num
          accu = accu + een_rescaled_e(k, i, j, nw) *         &
                        een_rescaled_n(m, a, i, nw)
          accu2 = accu2 + een_rescaled_e(k, i, j, nw) *       &
                          een_rescaled_n(m + l, a, i, nw)
          daccu(1:4) = daccu(1:4) + een_rescaled_e_deriv_e(k, j, 1:4, i, nw) *         &
                                    een_rescaled_n(m, a, i, nw)
          daccu2(1:4) = daccu2(1:4) + een_rescaled_e_deriv_e(k, j, 1:4, i, nw) *       &
                                      een_rescaled_n(m + l, a, i, nw)
        end do
        factor_een_deriv_e(j, 1:4, nw) = factor_een_deriv_e(j, 1:4, nw) +              &
               (accu * een_rescaled_n_deriv_e(m + l, a, 1:4, j, nw)                    &
                + daccu(1:4) * een_rescaled_n(m + l, a, j, nw)                         &
                + daccu2(1:4) * een_rescaled_n(m, a, j, nw)                            &
                + accu2 * een_rescaled_n_deriv_e(m, a, 1:4, j, nw)) * cn

        factor_een_deriv_e(j, 4, nw) = factor_een_deriv_e(j, 4, nw) + 2.0d0 * (        &
            daccu (1) * een_rescaled_n_deriv_e(m + l, a, 1, j, nw) +                    &
            daccu (2) * een_rescaled_n_deriv_e(m + l, a, 2, j, nw) +                    &
            daccu (3) * een_rescaled_n_deriv_e(m + l, a, 3, j, nw) +                    &
            daccu2(1) * een_rescaled_n_deriv_e(m, a, 1, j, nw    ) +                    &
            daccu2(2) * een_rescaled_n_deriv_e(m, a, 2, j, nw    ) +                    &
            daccu2(3) * een_rescaled_n_deriv_e(m, a, 3, j, nw    ) ) * cn

      end do
    end do
  end do
  end do

end function qmckl_compute_factor_een_deriv_e_naive_f

3.17.3 Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`walk_num`	`int64_t`	in	Number of walkers
`elec_num`	`int64_t`	in	Number of electrons
`nucl_num`	`int64_t`	in	Number of nucleii
`cord_num`	`int64_t`	in	order of polynomials
`dim_c_vector`	`int64_t`	in	dimension of full coefficient vector
`c_vector_full`	`double[dim_c_vector][nucl_num]`	in	full coefficient vector
`lkpm_combined_index`	`int64_t[4][dim_c_vector]`	in	combined indices
`tmp_c`	`double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]`	in	Temporary intermediate tensor
`dtmp_c`	`double[walk_num][0:cord_num-1][0:cord_num][nucl_num][4][elec_num]`	in	vector of non-zero coefficients
`een_rescaled_n`	`double[walk_num][0:cord_num][nucl_num][elec_num]`	in	Electron-nucleus rescaled factor
`een_rescaled_n_deriv_e`	`double[walk_num][0:cord_num][nucl_num][4][elec_num]`	in	Derivative of Electron-nucleus rescaled factor
`factor_een_deriv_e`	`double[walk_num][4][elec_num]`	out	Derivative of Electron-nucleus jastrow

integer function qmckl_compute_factor_een_deriv_e_f( &
     context, walk_num, elec_num, nucl_num, &
     cord_num, dim_c_vector, c_vector_full, lkpm_combined_index, &
     tmp_c, dtmp_c, een_rescaled_n, een_rescaled_n_deriv_e, factor_een_deriv_e)&
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer*8             , intent(in)  :: walk_num, elec_num, cord_num, nucl_num, dim_c_vector
  integer*8             , intent(in)  :: lkpm_combined_index(dim_c_vector,4)
  double precision      , intent(in)  :: c_vector_full(nucl_num, dim_c_vector)
  double precision      , intent(in)  :: tmp_c(elec_num, nucl_num,0:cord_num, 0:cord_num-1,  walk_num)
  double precision      , intent(in)  :: dtmp_c(elec_num, 4, nucl_num,0:cord_num, 0:cord_num-1,  walk_num)
  double precision      , intent(in)  :: een_rescaled_n(elec_num, nucl_num, 0:cord_num, walk_num)
  double precision      , intent(in)  :: een_rescaled_n_deriv_e(elec_num, 4, nucl_num, 0:cord_num, walk_num)
  double precision      , intent(out) :: factor_een_deriv_e(elec_num,4,walk_num)

  integer*8 :: i, a, j, l, k, p, m, n, nw, ii
  double precision :: accu, accu2, cn

  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 (cord_num < 0) then
     info = QMCKL_INVALID_ARG_5
     return
  endif

  factor_een_deriv_e = 0.0d0

  do nw =1, walk_num
  do n = 1, dim_c_vector
    l = lkpm_combined_index(n, 1)
    k = lkpm_combined_index(n, 2)
    p = lkpm_combined_index(n, 3)
    m = lkpm_combined_index(n, 4)

    do a = 1, nucl_num
      cn = c_vector_full(a, n)
      if(cn == 0.d0) cycle

      do ii = 1, 4
        do j = 1, elec_num
          factor_een_deriv_e(j,ii,nw) = factor_een_deriv_e(j,ii,nw) +                           (&
                                        tmp_c(j,a,m,k,nw)       * een_rescaled_n_deriv_e(j,ii,a,m+l,nw) + &
                                        (dtmp_c(j,ii,a,m,k,nw))   * een_rescaled_n(j,a,m+l,nw)            + &
                                        (dtmp_c(j,ii,a,m+l,k,nw)) * een_rescaled_n(j,a,m  ,nw)              + &
                                        tmp_c(j,a,m+l,k,nw)     * een_rescaled_n_deriv_e(j,ii,a,m,nw)     &
                                        ) * cn
        end do
      end do

      cn = cn + cn
      do j = 1, elec_num
        factor_een_deriv_e(j,4,nw) = factor_een_deriv_e(j,4,nw) +                              (&
                                      (dtmp_c(j,1,a,m  ,k,nw)) * een_rescaled_n_deriv_e(j,1,a,m+l,nw)  + &
                                      (dtmp_c(j,2,a,m  ,k,nw)) * een_rescaled_n_deriv_e(j,2,a,m+l,nw)  + &
                                      (dtmp_c(j,3,a,m  ,k,nw)) * een_rescaled_n_deriv_e(j,3,a,m+l,nw)  + &
                                      (dtmp_c(j,1,a,m+l,k,nw)) * een_rescaled_n_deriv_e(j,1,a,m  ,nw)  + &
                                      (dtmp_c(j,2,a,m+l,k,nw)) * een_rescaled_n_deriv_e(j,2,a,m  ,nw)  + &
                                      (dtmp_c(j,3,a,m+l,k,nw)) * een_rescaled_n_deriv_e(j,3,a,m  ,nw)    &
                                      ) * cn
      end do
    end do
  end do
  end do

end function qmckl_compute_factor_een_deriv_e_f

3.17.4 Test

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

double factor_een_deriv_e[4][walk_num][elec_num];
rc = qmckl_get_jastrow_factor_een_deriv_e(context, &(factor_een_deriv_e[0][0][0]),4*walk_num*elec_num);

assert(fabs(factor_een_deriv_e[0][0][0] + 0.0005481671107226865) < 1e-12);

Jastrow Factor

Table of Contents

1 Introduction

2 Context

2.1 Data structure

2.2 Initialization functions

2.2.0.1 Fortran interface

2.3 Access functions

2.3.0.1 Fortran interface

2.4 Test

3 Computation

3.1 Asymptotic component for \(J_{ee}\)

3.1.1 Get

3.1.1.1 Fortran interface

3.1.2 Compute

3.1.3 Test

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

3.2.1 Get

3.2.1.1 Fortran interface

3.2.2 Compute

3.2.3 Test

3.3 Electron-electron component derivative \(f'_{ee}\)

3.3.1 Get

3.3.2 Compute

3.3.3 Test

3.4 Asymptotic component for \(J_{eN}\)

3.4.1 Get

3.4.1.1 Fortran interface

3.4.2 Compute

3.4.3 Test

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

3.5.1 Get

3.5.1.1 Fortran interface

3.5.2 Compute

3.5.3 Test

3.6 Electron-nucleus component derivative \(f'_{en}\)

3.6.1 Get

3.6.2 Compute

3.6.3 Test

3.7 Electron-electron rescaled distances

3.7.1 Get

3.7.2 Compute

3.7.3 Test

3.8 Electron-electron rescaled distance gradients and Laplacian with respect to electron coords

3.8.1 Get

3.8.2 Compute

3.8.3 Test

3.9 Electron-electron-nucleus rescaled distances for each order

3.9.1 Get

3.9.2 Compute

3.9.3 Test

3.10 Electron-electron-nucleus rescaled distances for each order and derivatives

3.10.1 Get

3.10.2 Compute

3.10.3 Test

3.11 Electron-nucleus rescaled distances

3.11.1 Get

3.11.2 Compute

3.11.3 Test

3.12 Electron-electron-nucleus rescaled distance gradients and laplacian with respect to electron coords

3.12.1 Get

3.12.2 Compute

3.12.3 Test

3.13 Electron-electron-nucleus rescaled distances for each order

3.13.1 Get

3.13.2 Compute

3.13.3 Test

3.14 Electron-nucleus rescaled distances for each order and derivatives

3.14.1 Get

3.14.2 Compute

3.14.3 Test

3.15 Prepare for electron-electron-nucleus Jastrow \(f_{een}\)

3.15.1 Get

3.15.2 Compute dimcvector

3.15.3 Compute cvectorfull

3.15.4 Compute lkpmcombinedindex

3.15.5 Compute tmpc

3.15.6 Compute dtmpc

3.15.7 Test

3.16 Electron-electron-nucleus Jastrow \(f_{een}\)

3.15.2 Compute dim_c_vector

3.15.3 Compute c_vector_full

3.15.4 Compute lkpm_combined_index

3.15.5 Compute tmp_c

3.15.6 Compute dtmp_c