UP | HOME

Jastrow Factor

Table of Contents

1 Context

The following data stored in the context:

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

computed data:

int64_t dim_cord_vect Number of unique C coefficients
uint64_t dim_cord_vect_date Number of unique C coefficients
double asymp_jasb[2] Asymptotic component
uint64_t asymp_jasb_date Asymptotic component
double cord_vect_full[dim_cord_vect][nucl_num] vector of non-zero coefficients
uint64_t cord_vect_full_date Keep track of changes here
int64_t lkpm_combined_index[4][dim_cord_vect] Transform l,k,p, and m into consecutive indices
uint64_t lkpm_combined_index_date Transform l,k,p, and m into consecutive indices
double tmp_c[elec_num][nucl_num][ncord + 1][ncord][walk_num] vector of non-zero coefficients
double dtmp_c[elec_num][4][nucl_num][ncord + 1][ncord][walk_num] vector of non-zero coefficients
double een_rescaled_e[walk_num][elec_num][elec_num][0:cord_num] The electron-electron rescaled distances raised to the powers defined by cord
uint64_t een_rescaled_e_date Keep track of the date of creation
double een_rescaled_n[walk_num][elec_num][nucl_num][0:cord_num] The electron-electron rescaled distances raised to the powers defined by cord
uint64_t een_rescaled_n_date Keep track of the date of creation
double een_rescaled_e_deriv_e[walk_num][elec_num][4][elec_num][0:cord_num] The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons
uint64_t een_rescaled_e_deriv_e_date Keep track of the date of creation
double een_rescaled_n_deriv_e[walk_num][elec_num][4][nucl_num][0:cord_num] The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons
uint64_t een_rescaled_n_deriv_e_date Keep track of the date of creation

For H2O we have the following data: 

import numpy as np

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

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

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

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

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

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

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

kappa     = 1.0
kappa_inv = 1.0/kappa

1.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_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;
  int64_t* type_nucl_vector;
  double * aord_vector;
  double * bord_vector;
  double * cord_vector;
  double * asymp_jasb;
  double * factor_ee;
  double * factor_en;
  double * factor_een;
  double * factor_ee_deriv_e;
  double * factor_en_deriv_e;
  double * factor_een_deriv_e;
  int64_t  dim_cord_vect;
  uint64_t  dim_cord_vect_date;
  double * cord_vect_full;
  uint64_t  cord_vect_full_date;
  int64_t* lkpm_combined_index;
  uint64_t  lkpm_combined_index_date;
  double * tmp_c;
  double * dtmp_c;
  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;
} 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* const) context;
  assert (ctx != NULL);

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

  /* Default values */

  return QMCKL_SUCCESS;
}

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

1.3 Initialization functions

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

qmckl_exit_code  qmckl_set_jastrow_ord_num           (qmckl_context context, const int64_t aord_num, const int64_t bord_num, const int64_t cord_num);
qmckl_exit_code  qmckl_set_jastrow_type_nucl_num     (qmckl_context context, const int64_t type_nucl_num);
qmckl_exit_code  qmckl_set_jastrow_type_nucl_vector  (qmckl_context context, const int64_t* type_nucl_vector, const int64_t nucl_num);
qmckl_exit_code  qmckl_set_jastrow_aord_vector       (qmckl_context context, const double * aord_vector);
qmckl_exit_code  qmckl_set_jastrow_bord_vector       (qmckl_context context, const double * bord_vector);
qmckl_exit_code  qmckl_set_jastrow_cord_vector       (qmckl_context context, const double * cord_vector);
qmckl_exit_code  qmckl_set_jastrow_dependencies      (qmckl_context context);

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.

1.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;
double  rescale_factor_kappa_ee   = 1.0;
double  rescale_factor_kappa_en   = 1.0;
double  nucl_rescale_factor_kappa = 1.0;
double* elec_coord    = &(n2_elec_coord[0][0][0]);

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

/* Provide Electron data */

qmckl_exit_code rc;

assert(!qmckl_electron_provided(context));

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

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

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


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

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

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

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

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

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

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

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

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

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


int64_t w;
rc = qmckl_get_electron_walk_num (context, &w);
assert(rc == QMCKL_NOT_PROVIDED);


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

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

assert(qmckl_electron_provided(context));

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

double elec_coord2[walk_num*3*elec_num];

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


/* Provide Nucleus data */

assert(!qmckl_nucleus_provided(context));

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


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

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

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


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

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

double nucl_coord2[3*nucl_num];

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

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

assert(!qmckl_nucleus_provided(context));

rc = qmckl_get_nucleus_coord (context, 'N', nucl_coord2);
assert(rc == QMCKL_SUCCESS);
for (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_get_nucleus_coord (context, 'T', nucl_coord2);
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);
assert(rc == QMCKL_NOT_PROVIDED);

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

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

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

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

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

\[ J_{asymp} = \frac{b_1 \kappa^-1}{1 + b_2 \kappa^-1} \]

2.1.1 Get

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

2.1.2 Compute

qmcklcontext context in Global state
int64t bordnum in Number of electrons
double bordvector[bordnum + 1] in Number of walkers
double rescalefactorkappaee in Electron coordinates
double asympjasb[2] out Electron-electron distances 
integer function qmckl_compute_asymp_jasb_f(context, bord_num, bord_vector, rescale_factor_kappa_ee, asymp_jasb) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer*8             , intent(in)  :: bord_num
  double precision      , intent(in)  :: bord_vector(bord_num + 1)
  double precision      , intent(in)  :: rescale_factor_kappa_ee
  double precision      , intent(out) :: asymp_jasb(2)

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

  info = QMCKL_SUCCESS

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

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

  asym_one = bord_vector(1) * kappa_inv / (1.0d0 + bord_vector(2) * kappa_inv)
  asymp_jasb(:) = (/asym_one, 0.5d0 * asym_one/)

  do i = 1, 2
     x = kappa_inv
     do p = 2, bord_num
        x = x * kappa_inv
        asymp_jasb(i) = asymp_jasb(i) + bord_vector(p + 1) * x 
     end do
  end do

end function qmckl_compute_asymp_jasb_f

qmckl_exit_code qmckl_compute_asymp_jasb (
      const qmckl_context context,
      const int64_t bord_num,
      const double* bord_vector,
      const double rescale_factor_kappa_ee,
      double* const asymp_jasb );

2.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* aord_vector = &(n2_aord_vector[0][0]);
double* bord_vector = &(n2_bord_vector[0]);
double* cord_vector = &(n2_cord_vector[0][0]);

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

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

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

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

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

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

2.2.1 Get

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

2.2.2 Compute

qmcklcontext context in Global state
int64t walknum in Number of walkers
int64t elecnum in Number of electrons
int64t upnum in Number of alpha electrons
int64t bordnum in Number of coefficients
double bordvector[bordnum + 1] in List of coefficients
double eedistancerescaled[walknum][elecnum][elecnum] in Electron-electron distances
double asympjasb[2] in Electron-electron distances
double factoree[walknum] out Electron-electron distances 
integer function qmckl_compute_factor_ee_f(context, walk_num, elec_num, up_num, bord_num,            &
                                           bord_vector, ee_distance_rescaled, asymp_jasb, factor_ee) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer*8             , intent(in)  :: walk_num, elec_num, bord_num, up_num
  double precision      , intent(in)  :: bord_vector(bord_num + 1)
  double precision      , intent(in)  :: ee_distance_rescaled(walk_num, elec_num, elec_num)
  double precision      , intent(in)  :: asymp_jasb(2)
  double precision      , intent(out) :: factor_ee(walk_num)

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

  info = QMCKL_SUCCESS

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

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

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

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

  factor_ee = 0.0d0

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

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

        if(j .LE. up_num .OR. i .GT. up_num) then
          spin_fact = 0.5d0
          ipar = 2
        endif

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

     end do
  end do
  end do

end function qmckl_compute_factor_ee_f

qmckl_exit_code qmckl_compute_factor_ee (
      const qmckl_context context,
      const int64_t walk_num,
      const int64_t elec_num,
      const int64_t up_num,
      const int64_t bord_num,
      const double* bord_vector,
      const double* ee_distance_rescaled,
      const double* asymp_jasb,
      double* const factor_ee );

2.2.3 Test

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

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

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

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

TODO: Add equation

2.3.1 Get

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

2.3.2 Compute

qmcklcontext context in Global state
int64t walknum in Number of walkers
int64t elecnum in Number of electrons
int64t upnum in Number of alpha electrons
int64t bordnum in Number of coefficients
double bordvector[bordnum + 1] in List of coefficients
double eedistancerescaled[walknum][elecnum][elecnum] in Electron-electron distances
double eedistancerescaledderive[walknum][4][elecnum][elecnum] in Electron-electron distances
double asympjasb[2] in Electron-electron distances
double factoreederive[walknum][4][elecnum] out Electron-electron distances 
integer function qmckl_compute_factor_ee_deriv_e_f(context, walk_num, elec_num, up_num, bord_num,            &
                                           bord_vector, ee_distance_rescaled, ee_distance_rescaled_deriv_e,  &
                                           asymp_jasb, factor_ee_deriv_e) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer*8             , intent(in)  :: walk_num, elec_num, bord_num, up_num
  double precision      , intent(in)  :: bord_vector(bord_num + 1)
  double precision      , intent(in)  :: ee_distance_rescaled(walk_num, elec_num, elec_num)
  double precision      , intent(in)  :: ee_distance_rescaled_deriv_e(walk_num, 4, elec_num, elec_num)
  double precision      , intent(in)  :: asymp_jasb(2)
  double precision      , intent(out) :: factor_ee_deriv_e(elec_num,4,walk_num)

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

  info = QMCKL_SUCCESS

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

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

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

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

  factor_ee_deriv_e = 0.0d0
  third = 1.0d0 / 3.0d0

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

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

        if((i .LE. up_num .AND. j .LE. up_num ) .OR.  &
           (i .GT. up_num .AND. j .GT. up_num)) then
          spin_fact = 0.5d0
        endif

        lap1 = 0.0d0
        lap2 = 0.0d0
        lap3 = 0.0d0
        do ii = 1, 3
          x = ee_distance_rescaled(nw, i, j)
          if(abs(x) < 1.0d-18) cycle
          do p = 2, bord_num
            y = p * bord_vector(p + 1) * x
            pow_ser_g(ii) = pow_ser_g(ii) + y * dx(ii)
            lap1 = lap1 + (p - 1) * y * xinv * dx(ii) * dx(ii)
            lap2 = lap2 + y
            x = x * ee_distance_rescaled(nw, i, j)
          end do

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

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

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

     end do
  end do
  end do

end function qmckl_compute_factor_ee_deriv_e_f

qmckl_exit_code qmckl_compute_factor_ee_deriv_e (
      const qmckl_context context,
      const int64_t walk_num,
      const int64_t elec_num,
      const int64_t up_num,
      const int64_t bord_num,
      const double* bord_vector,
      const double* ee_distance_rescaled,
      const double* ee_distance_rescaled_deriv_e,
      const double* asymp_jasb,
      double* const factor_ee_deriv_e );

2.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]));

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

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

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

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

2.4.1 Get

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

2.4.2 Compute

qmcklcontext context in Global state
int64t walknum in Number of walkers
int64t elecnum in Number of electrons
int64t nuclnum in Number of nucleii
int64t typenuclnum in Number of unique nuclei
int64t typenuclvector[nuclnum] in IDs of unique nucleii
int64t aordnum in Number of coefficients
double aordvector[aordnum + 1][typenuclnum] in List of coefficients
double endistancerescaled[walknum][nuclnum][elecnum] in Electron-nucleus distances
double factoren[walknum] 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, aord_vector,                         &
                                           en_distance_rescaled, factor_en) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer*8             , intent(in)  :: walk_num, elec_num, aord_num, nucl_num, type_nucl_num
  integer*8             , intent(in)  :: type_nucl_vector(nucl_num)
  double precision      , intent(in)  :: aord_vector(aord_num + 1, type_nucl_num)
  double precision      , intent(in)  :: en_distance_rescaled(walk_num, nucl_num, elec_num)
  double precision      , intent(out) :: factor_en(walk_num)

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

  info = QMCKL_SUCCESS

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

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

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

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

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

  factor_en = 0.0d0

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

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

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

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

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

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

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

2.5.1 Get

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

2.5.2 Compute

qmcklcontext context in Global state
int64t walknum in Number of walkers
int64t elecnum in Number of electrons
int64t nuclnum in Number of nucleii
int64t typenuclnum in Number of unique nuclei
int64t typenuclvector[nuclnum] in IDs of unique nucleii
int64t aordnum in Number of coefficients
double aordvector[aordnum + 1][typenuclnum] in List of coefficients
double endistancerescaled[walknum][nuclnum][elecnum] in Electron-nucleus distances
double endistancerescaledderive[walknum][4][nuclnum][elecnum] in Electron-nucleus distance derivatives
double factorenderive[walknum][4][elecnum] 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, aord_vector,                         &
                                           en_distance_rescaled, en_distance_rescaled_deriv_e, factor_en_deriv_e) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer*8             , intent(in)  :: walk_num, elec_num, aord_num, nucl_num, type_nucl_num
  integer*8             , intent(in)  :: type_nucl_vector(nucl_num)
  double precision      , intent(in)  :: aord_vector(aord_num + 1, type_nucl_num)
  double precision      , intent(in)  :: en_distance_rescaled(walk_num, elec_num, nucl_num)
  double precision      , intent(in)  :: en_distance_rescaled_deriv_e(walk_num, 4, elec_num, nucl_num)
  double precision      , intent(out) :: factor_en_deriv_e(elec_num,4,walk_num)

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

  info = QMCKL_SUCCESS

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

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

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

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

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

  factor_en_deriv_e = 0.0d0
  third = 1.0d0 / 3.0d0

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

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

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

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

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

        end do

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

     end do
  end do
  end do

end function qmckl_compute_factor_en_deriv_e_f

qmckl_exit_code qmckl_compute_factor_en_deriv_e (
      const qmckl_context context,
      const int64_t walk_num,
      const int64_t elec_num,
      const int64_t nucl_num,
      const int64_t type_nucl_num,
      const int64_t* type_nucl_vector,
      const int64_t aord_num,
      const double* aord_vector,
      const double* en_distance_rescaled,
      const double* en_distance_rescaled_deriv_e,
      double* const factor_en_deriv_e );

2.5.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]));

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

2.6 Electron-electron 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.

2.6.1 Get

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

2.6.2 Compute

qmcklcontext context in Global state
int64t walknum in Number of walkers
int64t elecnum in Number of electrons
int64t cordnum in Order of polynomials
double rescalefactorkappaee in Factor to rescale ee distances
double eedistance[walknum][elecnum][elecnum] in Electron-electron distances
double eenrescalede[walknum][elecnum][elecnum][0:cordnum] out Electron-electron rescaled distances 
integer function qmckl_compute_een_rescaled_e_f(context, walk_num, elec_num, cord_num, rescale_factor_kappa_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_kappa_ee
  double precision      , intent(in)  :: ee_distance(elec_num,elec_num,walk_num)
  double precision      , intent(out) :: een_rescaled_e(0:cord_num,elec_num,elec_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_kappa_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(l, i, j, nw) = x
        een_rescaled_e(l, j, i, nw) = x
      end do
    end do
  end do

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

  end do

end function qmckl_compute_een_rescaled_e_f

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_kappa_ee,
      const double* ee_distance,
      double* const een_rescaled_e );

2.6.3 Test

assert(qmckl_electron_provided(context));


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

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

2.7 Electron-electron 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

2.7.1 Get

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

2.7.2 Compute

qmcklcontext context in Global state
int64t walknum in Number of walkers
int64t elecnum in Number of electrons
int64t cordnum in Order of polynomials
double rescalefactorkappaee in Factor to rescale ee distances
double coordnew[walknum][3][elecnum] in Electron coordinates
double eedistance[walknum][elecnum][elecnum] in Electron-electron distances
double eenrescalede[walknum][elecnum][elecnum][0:cordnum] in Electron-electron distances
double eenrescaledederive[walknum][elecnum][4][elecnum][0:cordnum] 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_kappa_ee,  &
     coord_new, 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_kappa_ee
  double precision      , intent(in)  :: coord_new(elec_num,3,walk_num)
  double precision      , intent(in)  :: ee_distance(elec_num,elec_num,walk_num)
  double precision      , intent(in)  :: een_rescaled_e(0:cord_num,elec_num,elec_num,walk_num)
  double precision      , intent(out) :: een_rescaled_e_deriv_e(0:cord_num,elec_num,4,elec_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_new(i, ii, nw) - coord_new(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_kappa_ee
      do j = 1, elec_num
        do i = 1, elec_num
          do ii = 1, 4
            een_rescaled_e_deriv_e(l, i, ii, j, nw) = kappa_l * elec_dist_deriv_e(ii, i, j)
          end do

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

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

end function qmckl_compute_factor_een_rescaled_e_deriv_e_f

qmckl_exit_code qmckl_compute_een_rescaled_e_deriv_e (
      const qmckl_context context,
      const int64_t walk_num,
      const int64_t elec_num,
      const int64_t cord_num,
      const double rescale_factor_kappa_ee,
      const double* coord_new,
      const double* ee_distance,
      const double* een_rescaled_e,
      double* const een_rescaled_e_deriv_e );

2.7.3 Test

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

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

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

2.8.1 Get

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

2.8.2 Compute

qmcklcontext context in Global state
int64t walknum in Number of walkers
int64t elecnum in Number of electrons
int64t nuclnum in Number of atoms
int64t cordnum in Order of polynomials
double rescalefactorkappaen in Factor to rescale ee distances
double endistance[walknum][elecnum][nuclnum] in Electron-nucleus distances
double eenrescaledn[walknum][elecnum][nuclnum][0:cordnum] out Electron-nucleus rescaled distances 
integer function qmckl_compute_een_rescaled_n_f(context, walk_num, elec_num, nucl_num, cord_num, rescale_factor_kappa_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)  :: cord_num
  double precision      , intent(in)  :: rescale_factor_kappa_en
  double precision      , intent(in)  :: en_distance(elec_num,nucl_num,walk_num)
  double precision      , intent(out) :: een_rescaled_n(0:cord_num,nucl_num,elec_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(1, a, i, nw) = dexp(-rescale_factor_kappa_en * 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(l, a, i, nw) = een_rescaled_n(l - 1, a, i, nw) * een_rescaled_n(1, a, i, 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 cord_num,
      const double rescale_factor_kappa_en,
      const double* en_distance,
      double* const een_rescaled_n );

2.8.3 Test

assert(qmckl_electron_provided(context));

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

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

2.9 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:

2.9.1 Get

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

2.9.2 Compute

qmcklcontext context in Global state
int64t walknum in Number of walkers
int64t elecnum in Number of electrons
int64t nuclnum in Number of atoms
int64t cordnum in Order of polynomials
double rescalefactorkappaen in Factor to rescale ee distances
double coordnew[walknum][3][elecnum] in Electron coordinates
double coord[3][nuclnum] in Nuclear coordinates
double endistance[walknum][elecnum][nuclnum] in Electron-nucleus distances
double eenrescaledn[walknum][elecnum][nuclnum][0:cordnum] in Electron-nucleus distances
double eenrescalednderive[walknum][elecnum][4][nuclnum][0:cordnum] out Electron-nucleus rescaled distances 
integer function qmckl_compute_factor_een_rescaled_n_deriv_e_f(context, walk_num, elec_num, nucl_num, &
     cord_num, rescale_factor_kappa_en,                                                               &
     coord_new, coord, 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)  :: cord_num
  double precision      , intent(in)  :: rescale_factor_kappa_en
  double precision      , intent(in)  :: coord_new(elec_num,3,walk_num)
  double precision      , intent(in)  :: coord(nucl_num,3)
  double precision      , intent(in)  :: en_distance(elec_num,nucl_num,walk_num)
  double precision      , intent(in)  :: een_rescaled_n(0:cord_num,nucl_num,elec_num,walk_num)
  double precision      , intent(out) :: een_rescaled_n_deriv_e(0:cord_num,nucl_num,4,elec_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_new(i, ii, nw) - coord(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
    kappa_l = - dble(l) * rescale_factor_kappa_en
    do a = 1, nucl_num
      do i = 1, elec_num
        do ii = 1, 4
          een_rescaled_n_deriv_e(l, a, ii, i, nw) = kappa_l * elnuc_dist_deriv_e(ii, i, a)
        end do

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

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

end function qmckl_compute_factor_een_rescaled_n_deriv_e_f

qmckl_exit_code qmckl_compute_factor_een_rescaled_n_deriv_e (
      const qmckl_context context,
      const int64_t walk_num,
      const int64_t elec_num,
      const int64_t nucl_num,
      const int64_t cord_num,
      const double rescale_factor_kappa_en,
      const double* coord_new,
      const double* coord,
      const double* en_distance,
      const double* een_rescaled_n,
      double* const een_rescaled_n_deriv_e );

2.9.3 Test

assert(qmckl_electron_provided(context));

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

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

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

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

2.10.1 Get

qmckl_exit_code qmckl_get_jastrow_dim_cord_vect(qmckl_context context, int64_t* const dim_cord_vect);
qmckl_exit_code qmckl_get_jastrow_cord_vect_full(qmckl_context context, double* const cord_vect_full);
qmckl_exit_code qmckl_get_jastrow_lkpm_combined_index(qmckl_context context, int64_t* const lkpm_combined_index);

2.10.2 Compute dimcordvect

qmcklcontext context in Global state
int64t cordnum in Order of polynomials
int64t dimcordvect out dimension of cordvectfull table 
integer function qmckl_compute_dim_cord_vect_f(context, cord_num, dim_cord_vect) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer*8             , intent(in)  :: cord_num
  integer*8             , intent(out) :: dim_cord_vect
  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_cord_vect = 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_cord_vect = dim_cord_vect + 1
      end do
    end do
  end do

end function qmckl_compute_dim_cord_vect_f

qmckl_exit_code qmckl_compute_dim_cord_vect (
      const qmckl_context context,
      const int64_t cord_num,
      int64_t* const dim_cord_vect );

2.10.3 Compute cordvectfull

qmcklcontext context in Global state
int64t nuclnum in Number of atoms
int64t dimcordvect in dimension of cord full table
int64t typenuclnum in dimension of cord full table
int64t typenuclvector[nuclnum] in dimension of cord full table
double cordvector[dimcordvect][typenuclnum] in dimension of cord full table
double cordvectfull[dimcordvect][nuclnum] out Full list of coefficients 
integer function qmckl_compute_cord_vect_full_f(context, nucl_num, dim_cord_vect, type_nucl_num,  &
     type_nucl_vector, cord_vector, cord_vect_full) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer*8             , intent(in)  :: nucl_num
  integer*8             , intent(in)  :: dim_cord_vect
  integer*8             , intent(in)  :: type_nucl_num
  integer*8             , intent(in)  :: type_nucl_vector(nucl_num)
  double precision      , intent(in)  :: cord_vector(type_nucl_num, dim_cord_vect)
  double precision      , intent(out) :: cord_vect_full(nucl_num,dim_cord_vect)
  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_cord_vect <= 0) then
     info = QMCKL_INVALID_ARG_5
     return
  endif


  do a = 1, nucl_num
    cord_vect_full(a,1:dim_cord_vect) = cord_vector(type_nucl_vector(a),1:dim_cord_vect)
  end do

end function qmckl_compute_cord_vect_full_f

qmckl_exit_code qmckl_compute_cord_vect_full (
      const qmckl_context context,
      const int64_t nucl_num,
      const int64_t dim_cord_vect,
      const int64_t type_nucl_num,
      const int64_t* type_nucl_vector,
      const double* cord_vector,
      double* const cord_vect_full );

2.10.4 Compute lkpmcombinedindex

qmcklcontext context in Global state
int64t cordnum in Order of polynomials
int64t dimcordvect in dimension of cord full table
int64t lpkmcombinedindex[4][dimcordvect] out Full list of combined indices 
integer function qmckl_compute_lkpm_combined_index_f(context, cord_num, dim_cord_vect,  &
     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_cord_vect
  integer*8             , intent(out) :: lkpm_combined_index(dim_cord_vect, 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_cord_vect <= 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_cord_vect,
      int64_t* const lpkm_combined_index );

2.10.5 Test

//assert(qmckl_electron_provided(context));
//

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

2.11.1 Get

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

2.11.2 Compute

qmcklcontext context in Global state
int64t walknum in Number of walkers
int64t elecnum in Number of electrons
int64t nuclnum in Number of nucleii
int64t cordnum in order of polynomials
int64t dimcordvect in dimension of full coefficient vector
double cordvectfull[dimcordvect][nuclnum] in full coefficient vector
int64t lkpmcombinedindex[4][dimcordvect] in combined indices
double eenrescalede[walknum][elecnum][elecnum][0:cordnum] in Electron-nucleus rescaled
double eenrescaledn[walknum][elecnum][nuclnum][0:cordnum] in Electron-nucleus rescaled factor
double factoreen[walknum] out Electron-nucleus jastrow 
integer function qmckl_compute_factor_een_f(context, walk_num, elec_num, nucl_num, cord_num, dim_cord_vect,            &
                                           cord_vect_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_cord_vect
  integer*8             , intent(in)  :: lkpm_combined_index(dim_cord_vect,4)
  double precision      , intent(in)  :: cord_vect_full(nucl_num, dim_cord_vect)
  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_cord_vect
    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 = cord_vect_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_f

qmckl_exit_code qmckl_compute_factor_een (
      const qmckl_context context,
      const int64_t walk_num,
      const int64_t elec_num,
      const int64_t nucl_num,
      const int64_t cord_num,
      const int64_t dim_cord_vect,
      const double* cord_vect_full,
      const int64_t* lkpm_combined_index,
      const double* een_rescaled_e,
      const double* een_rescaled_n,
      double* const factor_een );

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

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

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

2.12.1 Get

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

2.12.2 Compute

qmcklcontext context in Global state
int64t walknum in Number of walkers
int64t elecnum in Number of electrons
int64t nuclnum in Number of nucleii
int64t cordnum in order of polynomials
int64t dimcordvect in dimension of full coefficient vector
double cordvectfull[dimcordvect][nuclnum] in full coefficient vector
int64t lkpmcombinedindex[4][dimcordvect] in combined indices
double eenrescalede[walknum][elecnum][elecnum][0:cordnum] in Electron-nucleus rescaled
double eenrescaledn[walknum][elecnum][nuclnum][0:cordnum] in Electron-nucleus rescaled factor
double eenrescaledederive[walknum][elecnum][4][elecnum][0:cordnum] in Electron-nucleus rescaled
double eenrescalednderive[walknum][elecnum][4][nuclnum][0:cordnum] in Electron-nucleus rescaled factor
double factoreenderive[walknum][4][elecnum] out Electron-nucleus jastrow 
integer function qmckl_compute_factor_een_deriv_e_f(context, walk_num, elec_num, nucl_num, cord_num, dim_cord_vect,            &
                                           cord_vect_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_cord_vect
  integer*8             , intent(in)  :: lkpm_combined_index(dim_cord_vect, 4)
  double precision      , intent(in)  :: cord_vect_full(nucl_num, dim_cord_vect)
  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_cord_vect
    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 = cord_vect_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_f

qmckl_exit_code qmckl_compute_factor_een_deriv_e (
      const qmckl_context context,
      const int64_t walk_num,
      const int64_t elec_num,
      const int64_t nucl_num,
      const int64_t cord_num,
      const int64_t dim_cord_vect,
      const double* cord_vect_full,
      const int64_t* lkpm_combined_index,
      const double* een_rescaled_e,
      const double* een_rescaled_n,
      const double* een_rescaled_e_deriv_e,
      const double* een_rescaled_n_deriv_e,
      double* const factor_een_deriv_e );

2.12.3 Test

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

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

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

Author: TREX CoE

Created: 2021-10-11 Mon 22:18

Validate