qmckl/org/qmckl_jastrow.org

Jastrow Factor

• Context
  • Data structure
  • Access functions
  • Initialization functions
  • Test
• Computation
  • Asymptotic component for \(f_{ee}\)
    • Get
    • Compute
    • Test
  • Electron-electron component \(f_{ee}\)
    • Get
    • Compute
    • Test
  • Electron-electron component derivative \(f'_{ee}\)
    • Get
    • Compute
    • Test
  • Electron-nucleus component \(f_{en}\)
    • Get
    • Compute
    • Test
  • Electron-nucleus component derivative \(f'_{en}\)
    • Get
    • Compute
    • Test

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

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
double	factor_ee_date	out	Jastrow factor: electron-electron part
double	factor_en[walk_num]	out	Jastrow factor: electron-nucleus part
double	factor_en_date	out	Jastrow factor: electron-nucleus part
double	factor_een[walk_num]	out	Jastrow factor: electron-electron-nucleus part
double	factor_een_date	out	Jastrow factor: electron-electron-nucleus part
double	factor_ee_deriv_e[4][nelec][walk_num]	out	Derivative of the Jastrow factor: electron-electron-nucleus part
double	factor_ee_deriv_e_date	out	Keep track of the date for the derivative
double	factor_en_deriv_e[4][nelec][walk_num]	out	Derivative of the Jastrow factor: electron-electron-nucleus part
double	factor_en_deriv_e_date	out	Keep track of the date for the en derivative
double	factor_een_deriv_e[4][nelec][walk_num]	out	Derivative of the Jastrow factor: electron-electron-nucleus part
double	factor_een_deriv_e_date	out	Keep track of the date for the een derivative

computed data:

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

For H2O we have the following data:

import numpy as np

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

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

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

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

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

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

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

kappa     = 1.0
kappa_inv = 1.0/kappa

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;
int64_t  asymp_jasb_date;
int64_t  tmp_c_date;
int64_t  dtmp_c_date;
int64_t  factor_ee_date;
int64_t  factor_en_date;
int64_t  factor_een_date;
int64_t  factor_ee_deriv_e_date;
int64_t  factor_en_deriv_e_date;
int64_t  factor_een_deriv_e_date;
int64_t* type_nucl_vector;
double * aord_vector;
double * bord_vector;
double * cord_vector;
double * asymp_jasb;
double * factor_ee;
double * factor_en;
double * factor_een;
double * factor_ee_deriv_e;
double * factor_en_deriv_e;
double * factor_een_deriv_e;
int64_t  dim_cord_vec;
double * coord_vect_full;
double * tmp_c;
double * dtmp_c;
bool     provided;
char *   type;
} qmckl_jastrow_struct;

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

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

#+NAME:post

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

#+NAME:pre2

#+NAME:post2

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.

Test

/* Reference input data */
int64_t walk_num      = n2_walk_num;
int64_t elec_num      = n2_elec_num;
int64_t elec_up_num   = n2_elec_up_num;
int64_t elec_dn_num   = n2_elec_dn_num;
double  rescale_factor_kappa_ee   = 1.0;
double  rescale_factor_kappa_en   = 1.0;
double  nucl_rescale_factor_kappa = 1.0;
double* elec_coord    = &(n2_elec_coord[0][0][0]);

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

/* Provide Electron data */

qmckl_exit_code rc;

assert(!qmckl_electron_provided(context));

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

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

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


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

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

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

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

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

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

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

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

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

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


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


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

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

assert(qmckl_electron_provided(context));

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

double elec_coord2[walk_num*3*elec_num];

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


/* Provide Nucleus data */

assert(!qmckl_nucleus_provided(context));

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


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

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

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


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

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

double nucl_coord2[3*nucl_num];

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

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

assert(!qmckl_nucleus_provided(context));

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

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

double nucl_charge2[nucl_num];

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

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

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

Computation

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

Asymptotic component for \(f_{ee}\)

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

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

Get

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

Compute

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

integer function qmckl_compute_asymp_jasb_f(context, bord_num, bord_vector, rescale_factor_kappa_ee, asymp_jasb) &
 result(info)
use qmckl
implicit none
integer(qmckl_context), intent(in)  :: context
integer*8             , intent(in)  :: bord_num
double precision      , intent(in)  :: bord_vector(bord_num)
double precision      , intent(in)  :: rescale_factor_kappa_ee
double precision      , intent(out) :: asymp_jasb(2)

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

info = QMCKL_SUCCESS

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

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

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

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

end function qmckl_compute_asymp_jasb_f

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

Test

asym_one         :  0.43340325572525706
asymp_jasb[0]    :  0.5323750557252571
asymp_jasb[1]    :  0.31567342786262853

assert(qmckl_electron_provided(context));

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

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

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

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

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

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

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

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

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

Get

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

Compute

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

integer function qmckl_compute_factor_ee_f(context, walk_num, elec_num, up_num, bord_num,            &
                                       bord_vector, ee_distance_rescaled, asymp_jasb, factor_ee) &
 result(info)
use qmckl
implicit none
integer(qmckl_context), intent(in)  :: context
integer*8             , intent(in)  :: walk_num, elec_num, bord_num, up_num
double precision      , intent(in)  :: bord_vector(bord_num)
double precision      , intent(in)  :: ee_distance_rescaled(walk_num, elec_num, elec_num)
double precision      , intent(in)  :: asymp_jasb(2)
double precision      , intent(out) :: factor_ee(walk_num)

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

info = QMCKL_SUCCESS

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

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

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

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

factor_ee = 0.0d0

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

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

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

 end do
end do
end do

end function qmckl_compute_factor_ee_f

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

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

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

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

TODO: Add equation

Get

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

Compute

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

integer function qmckl_compute_factor_ee_deriv_e_f(context, walk_num, elec_num, up_num, bord_num,            &
                                       bord_vector, ee_distance_rescaled, ee_distance_rescaled_deriv_e,  &
	   asymp_jasb, factor_ee_deriv_e) &
 result(info)
use qmckl
implicit none
integer(qmckl_context), intent(in)  :: context
integer*8             , intent(in)  :: walk_num, elec_num, bord_num, up_num
double precision      , intent(in)  :: bord_vector(bord_num)
double precision      , intent(in)  :: ee_distance_rescaled(walk_num, elec_num, elec_num)
double precision      , intent(in)  :: ee_distance_rescaled_deriv_e(walk_num, 4, elec_num, elec_num)
double precision      , intent(in)  :: asymp_jasb(2)
double precision      , intent(out) :: factor_ee_deriv_e(elec_num,4,walk_num)

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

info = QMCKL_SUCCESS

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

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

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

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

factor_ee_deriv_e = 0.0d0
third = 1.0d0 / 3.0d0

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

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

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

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

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

 end do
end do
end do

end function qmckl_compute_factor_ee_deriv_e_f

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

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

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

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

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

Get

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

Compute

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

integer function qmckl_compute_factor_en_f(context, walk_num, elec_num, nucl_num, type_nucl_num,            &
                                       type_nucl_vector, aord_num, aord_vector,                         &
                                       en_distance_rescaled, factor_en) &
 result(info)
use qmckl
implicit none
integer(qmckl_context), intent(in)  :: context
integer*8             , intent(in)  :: walk_num, elec_num, aord_num, nucl_num, type_nucl_num
integer*8             , intent(in)  :: type_nucl_vector(type_nucl_num)
double precision      , intent(in)  :: aord_vector(aord_num, nucl_num)
double precision      , intent(in)  :: en_distance_rescaled(walk_num, elec_num, nucl_num)
double precision      , intent(out) :: factor_en(walk_num)

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

info = QMCKL_SUCCESS

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

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

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

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

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

factor_en = 0.0d0

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

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

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

 end do
end do
end do

end function qmckl_compute_factor_en_f

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

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

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

Calculate the electron-electron jastrow component factor_en_deriv_e derivative with respect to the electron coordinates using the en_distance_rescaled and en_distance_rescaled_deriv_e which are already calculated previously.

TODO: write equations.

Get

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

Compute

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

integer function qmckl_compute_factor_en_deriv_e_f(context, walk_num, elec_num, nucl_num, type_nucl_num,            &
                                       type_nucl_vector, aord_num, aord_vector,                         &
                                       en_distance_rescaled, en_distance_rescaled_deriv_e, factor_en_deriv_e) &
 result(info)
use qmckl
implicit none
integer(qmckl_context), intent(in)  :: context
integer*8             , intent(in)  :: walk_num, elec_num, aord_num, nucl_num, type_nucl_num
integer*8             , intent(in)  :: type_nucl_vector(type_nucl_num)
double precision      , intent(in)  :: aord_vector(aord_num, nucl_num)
double precision      , intent(in)  :: en_distance_rescaled(walk_num, elec_num, nucl_num)
double precision      , intent(in)  :: en_distance_rescaled_deriv_e(walk_num, 4, elec_num, nucl_num)
double precision      , intent(out) :: factor_en_deriv_e(elec_num,4,walk_num)

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

info = QMCKL_SUCCESS

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

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

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

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

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

factor_en_deriv_e = 0.0d0
third = 1.0d0 / 3.0d0

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

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

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

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

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

    end do

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

 end do
end do
end do

end function qmckl_compute_factor_en_deriv_e_f

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

Test

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

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

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