qmckl/org/qmckl_jastrow.org

Jastrow Factor

• Introduction
• Context
  • Data structure
  • Access functions
  • Initialization functions
  • Test
• Computation
  • Asymptotic component for \(J_{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
  • Electron-electron rescaled distances for each order
    • Get
    • Compute
    • Test
  • Electron-electron rescaled distances for each order and derivatives
    • Get
    • Compute
    • Test
  • Electron-nucleus rescaled distances for each order
    • Get
    • Compute
    • Test
  • Electron-nucleus rescaled distances for each order and derivatives
    • Get
    • Compute
    • Test
  • Prepare for electron-electron-nucleus Jastrow \(f_{een}\)
    • Get
    • Compute dim_cord_vect
    • Compute cord_vect_full
    • Compute lkpm_combined_index
    • Compute tmp_c
    • Compute tmp_c (OpenACC offload)
    • Compute tmp_c (cuBLAS offload)
    • Compute dtmp_c
    • Compute dtmp_c (OpenACC offload)
    • Compute dtmp_c (cuBLAS offload)
    • Test
  • Electron-electron-nucleus Jastrow \(f_{een}\)
    • Get
    • Compute naive
    • Compute
    • Test
  • Electron-electron-nucleus Jastrow \(f_{een}\) derivative
    • Get
    • Compute Naive
    • Compute
    • Test

Introduction

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

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

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

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

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

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

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

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

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

$f$ and $g$ are scaling function defined as

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

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

Context

The following data stored in the context:

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

computed data:

Variable	Type	In/Out
dim_cord_vect	int64_t	Number of unique C coefficients
dim_cord_vect_date	uint64_t	Number of unique C coefficients
asymp_jasb	double[2]	Asymptotic component
asymp_jasb_date	uint64_t	Asymptotic component
cord_vect_full	double[dim_cord_vect][nucl_num]	vector of non-zero coefficients
cord_vect_full_date	uint64_t	Keep track of changes here
lkpm_combined_index	int64_t[4][dim_cord_vect]	Transform l,k,p, and m into consecutive indices
lkpm_combined_index_date	uint64_t	Transform l,k,p, and m into consecutive indices
tmp_c	double[walk_num][cord_num][cord_num+1][nucl_num][elec_num]	vector of non-zero coefficients
dtmp_c	double[walk_num][elec_num][4][nucl_num][cord_num+1][cord_num]	vector of non-zero coefficients

een_rescaled_n	double[walk_num][cord_num+1][nucl_num][elec_num]	The electron-electron rescaled distances raised to the powers defined by cord
een_rescaled_n_date	uint64_t	Keep track of the date of creation
een_rescaled_e_deriv_e	double[walk_num][cord_num+1][elec_num][4][elec_num]	The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons
een_rescaled_e_deriv_e_date	uint64_t	Keep track of the date of creation
een_rescaled_n_deriv_e	double[walk_num][cord_num+1][nucl_num][4][elec_num]	The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons
een_rescaled_n_deriv_e_date	uint64_t	Keep track of the date of creation

Data structure

typedef enum qmckl_jastrow_offload_type{
  OFFLOAD_NONE,
  OFFLOAD_OPENACC,
  OFFLOAD_CUBLAS
} qmckl_jastrow_offload_type;

typedef struct qmckl_jastrow_struct{
int32_t  uninitialized;
int64_t  aord_num;
int64_t  bord_num;
int64_t  cord_num;
int64_t  type_nucl_num;
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_offload_type offload_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*) context;
assert (ctx != NULL);

ctx->jastrow.uninitialized = (1 << 5) - 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, const int64_t size_max);
qmckl_exit_code  qmckl_set_jastrow_bord_vector       (qmckl_context context, const double * bord_vector, const int64_t size_max);
qmckl_exit_code  qmckl_set_jastrow_cord_vector       (qmckl_context context, const double * cord_vector, const int64_t size_max);
qmckl_exit_code  qmckl_set_jastrow_offload_type      (qmckl_context context, const qmckl_jastrow_offload_type offload_type);

#+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*  nucl_coord    = &(n2_nucl_coord[0][0]);
int64_t size_max;

/* 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, walk_num*3*elec_num);
assert(rc == QMCKL_SUCCESS);

double elec_coord2[walk_num*3*elec_num];

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


/* Provide Nucleus data */

assert(!qmckl_nucleus_provided(context));

rc = qmckl_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, 3*nucl_num);
assert(rc == QMCKL_NOT_PROVIDED);

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

assert(!qmckl_nucleus_provided(context));

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

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

double nucl_charge2[nucl_num];

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

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

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

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 \(J_{ee}\)

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

\[ J_{\text{ee}}^{\infty} = \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,
                         const int64_t size_max);

Compute

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

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

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

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

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

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

  return QMCKL_SUCCESS;
}

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

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]);
int64_t dim_cord_vect=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,(aord_num+1)*type_nucl_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_jastrow_bord_vector(context, bord_vector,(bord_num+1));
assert(rc == QMCKL_SUCCESS);
rc = qmckl_get_jastrow_dim_cord_vect(context, &dim_cord_vect);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_jastrow_cord_vector(context, cord_vector,dim_cord_vect*type_nucl_num);
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,2);

// 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,
                        const int64_t size_max);

Compute

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

integer function qmckl_compute_factor_ee_f(context, walk_num, elec_num, up_num, bord_num,            &
                                       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(elec_num, elec_num, walk_num)
double precision      , intent(in)  :: asymp_jasb(2)
double precision      , intent(out) :: factor_ee(walk_num)

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

info = QMCKL_SUCCESS

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

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

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

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

factor_ee = 0.0d0

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

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

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

    factor_ee(nw) = factor_ee(nw) + spin_fact * bord_vector(1)  * &
                            ee_distance_rescaled(i,j,nw) / &
                            (1.0d0 + bord_vector(2) *   &
                            ee_distance_rescaled(i,j,nw))  &
                           -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 ) {

  int ipar; // can we use a smaller integer?
  double x, x1, spin_fact, power_ser;

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

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

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

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

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

        if(i < up_num || j >= up_num) {
          spin_fact = 0.5;
          ipar = 1;
        }
        
        factor_ee[nw] = factor_ee[nw] + spin_fact * bord_vector[0]  * \
                                x1 /  \
                                (1.0 + bord_vector[1] *               \
                                x1)   \
                               -asymp_jasb[ipar] + power_ser;

      }
    }
  }

  return QMCKL_SUCCESS;
}

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

// 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,
                                const int64_t size_max);

Compute

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

integer function qmckl_compute_factor_ee_deriv_e_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(elec_num, elec_num,walk_num)
double precision      , intent(in)  :: ee_distance_rescaled_deriv_e(4,elec_num, elec_num,walk_num)   !TODO
double precision      , intent(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(i,j,nw)
    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

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

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

    lap1 = 0.0d0
    lap2 = 0.0d0
    lap3 = 0.0d0
    do ii = 1, 3
      x = ee_distance_rescaled(i, j, nw)
      if(abs(x) < 1.0d-18) cycle
      do p = 2, bord_num
        y = p * 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(i, j, nw)
      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]),walk_num*4*elec_num);

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

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,
                        const int64_t size_max);

Compute

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

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

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

info = QMCKL_SUCCESS

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

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

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

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

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

factor_en = 0.0d0

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

    do p = 2, aord_num
      x = x * en_distance_rescaled(i, a, nw)
      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(i, a, nw) /               &
                            (1.0d0 + aord_vector(2, type_nucl_vector(a)) * &
                            en_distance_rescaled(i, a, nw))                &
                            + 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 ) {

double  x, x1, power_ser;


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

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

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

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

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

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

if (aord_num <= 0) {
 return QMCKL_INVALID_ARG_7;
}

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

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

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


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

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

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

  }
}
}

return QMCKL_SUCCESS;
}

qmckl_exit_code qmckl_compute_factor_en (
const qmckl_context context,
const int64_t walk_num,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t type_nucl_num,
const int64_t* type_nucl_vector,
const int64_t aord_num,
const double* 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,walk_num);

// 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,
                                const int64_t size_max);

Compute

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

integer function qmckl_compute_factor_en_deriv_e_f( &
 context, walk_num, elec_num, nucl_num, type_nucl_num, &
 type_nucl_vector, aord_num, 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(elec_num, nucl_num, walk_num)
double precision      , intent(in)  :: en_distance_rescaled_deriv_e(4, elec_num, nucl_num, walk_num)
double precision      , intent(out) :: factor_en_deriv_e(elec_num,4,walk_num)

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

info = QMCKL_SUCCESS

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

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

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

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

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

factor_en_deriv_e = 0.0d0
third = 1.0d0 / 3.0d0

do nw =1, walk_num
do a = 1, nucl_num
 do i = 1, elec_num
    x = en_distance_rescaled(i,a,nw)
    if(abs(x) < 1.0d-18) continue
    power_ser_g = 0.0d0
    den = 1.0d0 + 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(ii,i,a,nw)
    end do

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

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

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

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.

Get

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

Compute

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

integer function qmckl_compute_een_rescaled_e_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(elec_num,elec_num,0:cord_num,walk_num)
double precision,dimension(:,:),allocatable :: een_rescaled_e_ij
double precision                    :: x
integer*8                           :: i, j, k, l, nw

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


info = QMCKL_SUCCESS

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

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

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

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

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

k = 0
do j = 1, elec_num
do i = 1, j - 1
  k = k + 1
  een_rescaled_e_ij(k, 2) = dexp(-rescale_factor_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(i, j, l, nw) = x
    een_rescaled_e(j, i, l, nw) = x
  end do
end do
end do

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

end do

end function qmckl_compute_een_rescaled_e_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 );

Test

assert(qmckl_electron_provided(context));


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

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

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

Get

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

Compute

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

integer function qmckl_compute_factor_een_rescaled_e_deriv_e_f( &
 context, walk_num, elec_num, cord_num, rescale_factor_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(elec_num,elec_num,0:cord_num,walk_num)
double precision      , intent(out) :: een_rescaled_e_deriv_e(elec_num,4,elec_num,0:cord_num,walk_num)
double precision,dimension(:,:,:),allocatable  :: elec_dist_deriv_e
double precision                    :: x, rij_inv, kappa_l
integer*8                           :: i, j, k, l, nw, ii

allocate(elec_dist_deriv_e(4,elec_num,elec_num))

info = QMCKL_SUCCESS

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

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

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

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

! Prepare table of exponentiated distances raised to appropriate power
een_rescaled_e_deriv_e     = 0.0d0
do nw = 1, walk_num
do j = 1, elec_num
  do i = 1, elec_num
    rij_inv = 1.0d0 / ee_distance(i, j, nw)
    do ii = 1, 3
      elec_dist_deriv_e(ii, i, j) = (coord_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
      een_rescaled_e_deriv_e(i, 1, j, l, nw) = kappa_l * elec_dist_deriv_e(1, i, j)
      een_rescaled_e_deriv_e(i, 2, j, l, nw) = kappa_l * elec_dist_deriv_e(2, i, j)
      een_rescaled_e_deriv_e(i, 3, j, l, nw) = kappa_l * elec_dist_deriv_e(3, i, j)
      een_rescaled_e_deriv_e(i, 4, j, l, nw) = kappa_l * elec_dist_deriv_e(4, i, j)

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

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

end function qmckl_compute_factor_een_rescaled_e_deriv_e_f

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

Test

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

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

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.

Get

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

Compute

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

integer function qmckl_compute_een_rescaled_n_f( &
 context, walk_num, elec_num, nucl_num, 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(elec_num,nucl_num,0:cord_num,walk_num)
double precision                    :: x
integer*8                           :: i, a, k, l, nw

info = QMCKL_SUCCESS

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

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

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

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

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

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

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

do a = 1, nucl_num
do i = 1, elec_num
  een_rescaled_n(i, a, 1, nw) = dexp(-rescale_factor_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(i, a, l, nw) = een_rescaled_n(i, a, l - 1, nw) * een_rescaled_n(i, a, 1, nw)
  end do
end do
end do
end do

end function qmckl_compute_een_rescaled_n_f

qmckl_exit_code qmckl_compute_een_rescaled_n (
const qmckl_context context,
const int64_t walk_num,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t cord_num,
const double rescale_factor_kappa_en,
const double* en_distance,
double* const een_rescaled_n ) {


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

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

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

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

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

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

for (int nw = 0; nw < walk_num; ++nw) {
for (int a = 0; a < nucl_num; ++a) {
  for (int i = 0; i < elec_num; ++i) {
    // prepare the actual een table
    //een_rescaled_n(:, :, 0, nw) = 1.0d0
    een_rescaled_n[i + a * elec_num + 0 + nw * elec_num*nucl_num*(cord_num+1)] = 1.0;
    //een_rescaled_n(i, a, 1, nw) = dexp(-rescale_factor_kappa_en * en_distance(i, a, nw))
    een_rescaled_n[i + a*elec_num + elec_num*nucl_num + nw*elec_num*nucl_num*(cord_num+1)] = exp(-rescale_factor_kappa_en * \
                                                                          en_distance[i + a*elec_num + nw*elec_num*nucl_num]);
  }
}

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

}

return QMCKL_SUCCESS;
}

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

Test

assert(qmckl_electron_provided(context));

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

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

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:

Get

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

Compute

Variable	Type	In/Out	Description
context	qmckl_context	in	Global state
walk_num	int64_t	in	Number of walkers
elec_num	int64_t	in	Number of electrons
nucl_num	int64_t	in	Number of atoms
cord_num	int64_t	in	Order of polynomials
rescale_factor_kappa_en	double	in	Factor to rescale ee distances
coord_new	double[walk_num][3][elec_num]	in	Electron coordinates
coord	double[3][nucl_num]	in	Nuclear coordinates
en_distance	double[walk_num][elec_num][nucl_num]	in	Electron-nucleus distances
een_rescaled_n	double[walk_num][0:cord_num][nucl_num][elec_num]	in	Electron-nucleus distances
een_rescaled_n_deriv_e	double[walk_num][0:cord_num][nucl_num][4][elec_num]	out	Electron-nucleus rescaled distances

integer function qmckl_compute_factor_een_rescaled_n_deriv_e_f( &
 context, walk_num, elec_num, nucl_num, &
 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(elec_num,nucl_num,0:cord_num,walk_num)
double precision      , intent(out) :: een_rescaled_n_deriv_e(elec_num,4,nucl_num,0:cord_num,walk_num)
double precision,dimension(:,:,:),allocatable :: elnuc_dist_deriv_e
double precision                    :: x, ria_inv, kappa_l
integer*8                           :: i, a, k, l, nw, ii

allocate(elnuc_dist_deriv_e(4, elec_num, nucl_num))

info = QMCKL_SUCCESS

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

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

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

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

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

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

! prepare the actual een table
do a = 1, nucl_num
do i = 1, elec_num
  ria_inv = 1.0d0 / en_distance(i, a, nw)
  do ii = 1, 3
    elnuc_dist_deriv_e(ii, i, a) = (coord_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
    een_rescaled_n_deriv_e(i, 1, a, l, nw) = kappa_l * elnuc_dist_deriv_e(1, i, a)
    een_rescaled_n_deriv_e(i, 2, a, l, nw) = kappa_l * elnuc_dist_deriv_e(2, i, a)
    een_rescaled_n_deriv_e(i, 3, a, l, nw) = kappa_l * elnuc_dist_deriv_e(3, i, a)
    een_rescaled_n_deriv_e(i, 4, a, l, nw) = kappa_l * elnuc_dist_deriv_e(4, i, a)

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

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

end function qmckl_compute_factor_een_rescaled_n_deriv_e_f

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

Test

assert(qmckl_electron_provided(context));

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

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

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.

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);
qmckl_exit_code qmckl_get_jastrow_tmp_c(qmckl_context context, double* const tmp_c);
qmckl_exit_code qmckl_get_jastrow_dtmp_c(qmckl_context context, double* const dtmp_c);

Compute dim_cord_vect

Variable	Type	In/Out	Description
context	qmckl_context	in	Global state
cord_num	int64_t	in	Order of polynomials
dim_cord_vect	int64_t	out	dimension of cord_vect_full 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){

int         lmax;


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

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

*dim_cord_vect = 0;

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

return QMCKL_SUCCESS;
}

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

Compute cord_vect_full

Variable	Type	In/Out	Description
context	qmckl_context	in	Global state
nucl_num	int64_t	in	Number of atoms
dim_cord_vect	int64_t	in	dimension of cord full table
type_nucl_num	int64_t	in	dimension of cord full table
type_nucl_vector	int64_t[nucl_num]	in	dimension of cord full table
cord_vector	double[dim_cord_vect][type_nucl_num]	in	dimension of cord full table
cord_vect_full	double[dim_cord_vect][nucl_num]	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 );

Compute lkpm_combined_index

Variable	Type	In/Out	Description
context	qmckl_context	in	Global state
cord_num	int64_t	in	Order of polynomials
dim_cord_vect	int64_t	in	dimension of cord full table
lkpm_combined_index	int64_t[4][dim_cord_vect]	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 lkpm_combined_index ) {

int kk, lmax, m;

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

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

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

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

return QMCKL_SUCCESS;
}

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

Compute tmp_c

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

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

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

info = QMCKL_SUCCESS

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

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

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

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

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

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

end function qmckl_compute_tmp_c_doc_f

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

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

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

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

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

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

qmckl_exit_code info = QMCKL_SUCCESS;

const char  TransA = 'N';
const char  TransB = 'N';
const double alpha = 1.0;
const double beta  = 0.0;

const int64_t M = elec_num;
const int64_t N = nucl_num*(cord_num + 1);
const int64_t K = elec_num;

const int64_t LDA = elec_num;
const int64_t LDB = elec_num;
const int64_t LDC = elec_num;

const int64_t af = elec_num*elec_num;
const int64_t bf = elec_num*nucl_num*(cord_num+1);
const int64_t cf = bf;

for (int64_t nw=0; nw < walk_num; ++nw) {
for (int64_t i=0; i<cord_num; ++i){
info = qmckl_dgemm(context, TransA, TransB, M, N, K, alpha,     \
                   &(een_rescaled_e[af*(i+nw*(cord_num+1))]), \
                   LDA,                                        \
                   &(een_rescaled_n[bf*nw]), \
                   LDB,                                        \
                   beta,                                       \
                   &(tmp_c[cf*(i+nw*cord_num)]), \
                   LDC);
}
}

return info;
}

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

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

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

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

Compute tmp_c (OpenACC offload)

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

#ifdef HAVE_OPENACC_OFFLOAD
qmckl_exit_code qmckl_compute_tmp_c_acc_offload (
const qmckl_context context,
const int64_t cord_num,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t walk_num,
const double* een_rescaled_e,
const double* een_rescaled_n,
double* const tmp_c ) {

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

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

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

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

// Compute array access strides:
// For tmp_c...
int stride_k_c  = elec_num;
int stride_j_c  = stride_k_c * nucl_num;
int stride_i_c  = stride_j_c * (cord_num+1);
int stride_nw_c = stride_i_c * cord_num;
// For een_rescaled_e...
int stride_m_e  = elec_num;
int stride_i_e  = stride_m_e * elec_num;
int stride_nw_e = stride_i_e * (cord_num+1);
// For een_rescaled_n...
int stride_k_n  = elec_num;
int stride_j_n  = stride_k_n * nucl_num;
int stride_nw_n = stride_j_n * (cord_num+1);

#pragma acc parallel
#pragma acc loop independent gang worker vector collapse(5)
for (int nw=0; nw < walk_num; ++nw) {
for (int i=0; i<cord_num; ++i){

  // Replacement for single DGEMM
  for (int j=0; j<cord_num+1; j++) {
    for (int k=0; k<nucl_num; k++) {
      for (int l=0; l<elec_num; l++) {

        // Single reduction
 tmp_c[l + k*stride_k_c + j*stride_j_c + i*stride_i_c + nw*stride_nw_c] = 0;
        for (int m=0; m<elec_num; m++) {
          tmp_c[l + k*stride_k_c + j*stride_j_c + i*stride_i_c + nw*stride_nw_c] =
          tmp_c[l + k*stride_k_c + j*stride_j_c + i*stride_i_c + nw*stride_nw_c] +
          een_rescaled_e[l + m*stride_m_e + i*stride_i_e + nw*stride_nw_e] *
          een_rescaled_n[m + k*stride_k_n + j*stride_j_n + nw*stride_nw_n];
        }

      }
    }
  }
}
}

return QMCKL_SUCCESS;
}
#endif

#ifdef HAVE_OPENACC_OFFLOAD
qmckl_exit_code qmckl_compute_tmp_c_acc_offload (
const qmckl_context context,
const int64_t cord_num,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t walk_num,
const double* een_rescaled_e,
const double* een_rescaled_n,
double* const tmp_c );
#endif

TODO: FIX dtmp_c dimension in the table.

Compute tmp_c (cuBLAS offload)

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

#ifdef HAVE_CUBLAS_OFFLOAD
qmckl_exit_code qmckl_compute_tmp_c_cublas_offload (
const qmckl_context context,
const int64_t cord_num,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t walk_num,
const double* een_rescaled_e,
const double* een_rescaled_n,
double* const tmp_c ) {

qmckl_exit_code info;
int          i, j, a, l, kk, p, lmax, nw;
char         TransA, TransB;
double       alpha, beta;
int          M, N, K, LDA, LDB, LDC;

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

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

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

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

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

M = elec_num;
N = nucl_num*(cord_num + 1);
K = elec_num;

LDA = sizeof(een_rescaled_e)/sizeof(double);
LDB = sizeof(een_rescaled_n)/sizeof(double);
LDC = sizeof(tmp_c)/sizeof(double);

// TODO Replace with cuBLAS calls
for (int nw=0; nw < walk_num; ++nw) {
for (int i=0; i<cord_num; ++i){
info = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha,     \
             //    &een_rescaled_e[0+0*elec_num+i*elec_num*elec_num+nw*elec_num*elec_num*(cord_num+1)],
                   &een_rescaled_e[             i*elec_num*elec_num+nw*elec_num*elec_num*(cord_num+1)], \
                   LDA,                                        \
             //    &een_rescaled_n[0+0*elec_num+0*elec_num*nucl_num+nw*elec_num*nucl_num*(cord_num+1)],
                   &een_rescaled_n[                                 nw*elec_num*nucl_num*(cord_num+1)], \
                   LDB,                                        \
                   beta,                                       \
             //    &tmp_c[0+0*elec_num+0*elec_num*nucl_num+i*elec_num*nucl_num*(cord_num+1)+nw*elec_num*nucl_num*(cord_num+1)*cord_num],
                   &tmp_c[                                 i*elec_num*nucl_num*(cord_num+1)+nw*elec_num*nucl_num*(cord_num+1)*cord_num],        \
                   LDC);

}
}

return info;
}
#endif

#ifdef HAVE_CUBLAS_OFFLOAD
qmckl_exit_code qmckl_compute_tmp_c_cublas_offload (
const qmckl_context context,
const int64_t cord_num,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t walk_num,
const double* een_rescaled_e,
const double* een_rescaled_n,
double* const tmp_c );
#endif

TODO: FIX dtmp_c dimension in the table.

Compute dtmp_c

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

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

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

info = QMCKL_SUCCESS

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

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

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

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

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

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

end function qmckl_compute_dtmp_c_doc_f

qmckl_exit_code qmckl_compute_dtmp_c_hpc (const qmckl_context context,
                                      const int64_t cord_num,
                                      const int64_t elec_num,
                                      const int64_t nucl_num,
                                      const int64_t walk_num,
                                      const double* een_rescaled_e_deriv_e,
                                      const double* een_rescaled_n,
                                      double* const dtmp_c )
{

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

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

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

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

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

qmckl_exit_code  info = QMCKL_SUCCESS;

const char  TransA = 'N';
const char  TransB = 'N';
const double alpha = 1.0;
const double beta  = 0.0;

const int64_t M = 4*elec_num;
const int64_t N = nucl_num*(cord_num + 1);
const int64_t K = elec_num;

const int64_t LDA = 4*elec_num;
const int64_t LDB = elec_num;
const int64_t LDC = 4*elec_num;

const int64_t af = elec_num*elec_num*4;
const int64_t bf = elec_num*nucl_num*(cord_num+1);
const int64_t cf = elec_num*4*nucl_num*(cord_num+1);

for (int64_t nw=0; nw < walk_num; ++nw) {
for (int64_t i=0; i < cord_num; ++i) {
  info = qmckl_dgemm(context, TransA, TransB, M, N, K, alpha,        \
            &(een_rescaled_e_deriv_e[af*(i+nw*(cord_num+1))]), \
            LDA,                                               \
            &(een_rescaled_n[bf*nw]),   \
            LDB,                                               \
            beta,                                              \
            &(dtmp_c[cf*(i+nw*cord_num)]), \
            LDC);
}
}

return info;
}

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

qmckl_exit_code qmckl_compute_dtmp_c (const qmckl_context context,
                                  const int64_t cord_num,
                                  const int64_t elec_num,
                                  const int64_t nucl_num,
                                  const int64_t walk_num,
                                  const double* een_rescaled_e_deriv_e,
                                  const double* een_rescaled_n,
                                  double* const dtmp_c );

qmckl_exit_code qmckl_compute_dtmp_c_doc (
const qmckl_context context,
const int64_t cord_num,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t walk_num,
const double* een_rescaled_e_deriv_e,
const double* een_rescaled_n,
double* const dtmp_c );

qmckl_exit_code qmckl_compute_dtmp_c_hpc (
const qmckl_context context,
const int64_t cord_num,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t walk_num,
const double* een_rescaled_e_deriv_e,
const double* een_rescaled_n,
double* const dtmp_c );

Compute dtmp_c (OpenACC offload)

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

#ifdef HAVE_OPENACC_OFFLOAD
qmckl_exit_code qmckl_compute_dtmp_c_acc_offload (
const qmckl_context context,
const int64_t cord_num,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t walk_num,
const double* een_rescaled_e_deriv_e,
const double* een_rescaled_n,
double* const dtmp_c ) {

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

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

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

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

// Compute strides...
// For dtmp_c
int stride_l_d  = elec_num;
int stride_k_d  = stride_l_d * 4;
int stride_j_d  = stride_k_d * nucl_num;
int stride_i_d  = stride_j_d * (cord_num+1);
int stride_nw_d = stride_i_d * cord_num;
// For een_rescaled_e_deriv_e
int stride_l_e  = elec_num;
int stride_n_e  = stride_l_e * 4;
int stride_i_e  = stride_n_e * elec_num;
int stride_nw_e = stride_i_e * cord_num;
// For een_rescaled_n
int stride_k_n  = elec_num;
int stride_j_n  = stride_k_n * nucl_num;
int stride_nw_n = stride_j_n * (cord_num+1);


#pragma acc parallel
#pragma loop independent gang worker vector collapse(6)
for (int nw=0; nw < walk_num; nw++) {
for (int i=0; i < cord_num; i++) {

  // Single DGEMM
  for(int j=0; j<cord_num+1; j++) {
    for(int k=0; k<nucl_num; k++) {
      for(int l=0; l<4; l++) {
        for(int m=0; m<elec_num; m++) {

          // Single reduction
dtmp_c[m + l * stride_l_d + k * stride_k_d + j * stride_j_d + i * stride_i_d + nw * stride_nw_d] = 0;
          for(int n=0; n<elec_num; n++){
            dtmp_c[m + l * stride_l_d + k * stride_k_d + j * stride_j_d + i * stride_i_d + nw * stride_nw_d] =
            dtmp_c[m + l * stride_l_d + k * stride_k_d + j * stride_j_d + i * stride_i_d + nw * stride_nw_d] +
            een_rescaled_e_deriv_e[m + l * stride_l_e +  n * stride_n_e + i * stride_i_e + nw * stride_nw_e] *
            een_rescaled_n[n + k * stride_k_n + j * stride_j_n + nw * stride_nw_n];
          }
        }
      }
    }
  }
}
}

return QMCKL_SUCCESS;
}
#endif

#ifdef HAVE_OPENACC_OFFLOAD
qmckl_exit_code qmckl_compute_dtmp_c_acc_offload (
const qmckl_context context,
const int64_t cord_num,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t walk_num,
const double* een_rescaled_e_deriv_e,
const double* een_rescaled_n,
double* const dtmp_c );
#endif

Compute dtmp_c (cuBLAS offload)

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

#ifdef HAVE_CUBLAS_OFFLOAD
qmckl_exit_code qmckl_compute_dtmp_c_cublas_offload (
const qmckl_context context,
const int64_t cord_num,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t walk_num,
const double* een_rescaled_e_deriv_e,
const double* een_rescaled_n,
double* const dtmp_c ) {

qmckl_exit_code  info;
char             TransA, TransB;
double           alpha, beta;
int              M, N, K, LDA, LDB, LDC;

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

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

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

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

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

M = 4*elec_num;
N = nucl_num*(cord_num + 1);
K = elec_num;

LDA = 4*sizeof(een_rescaled_e_deriv_e)/sizeof(double);
LDB = sizeof(een_rescaled_n)/sizeof(double);
LDC = 4*sizeof(dtmp_c)/sizeof(double);

// TODO Replace with cuBLAS calls
for (int nw=0; nw < walk_num; ++nw) {
for (int i=0; nw < cord_num; ++i) {
  info = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha,        \
          //&een_rescaled_e_deriv_e[0+0*elec_num+0*elec_num*4+i*elec_num*4*elec_num+nw*elec_num*4*elec_num*(cord_num+1)],
            &een_rescaled_e_deriv_e[i*elec_num*4*elec_num+nw*elec_num*4*elec_num*(cord_num+1)],                \
            LDA,                                               \
          //&een_rescaled_n[0+0*elec_num+0*elec_num*nucl_num+nw*elec_num*nucl_num*(cord_num+1)],
            &een_rescaled_n[nw*elec_num*nucl_num*(cord_num+1)],   \
            LDB,                                               \
            beta,                                              \
          //&dtmp_c[0+0*elec_num+0*elec_num*4+0*elec_num*4*nucl_num+i*elec_num*4*nucl_num*(cord_num+1)+nw*elec_num*4*nucl_num*(cord_num+1)*cord_num],
            &dtmp_c[i*elec_num*4*nucl_num*(cord_num+1)+nw*elec_num*4*nucl_num*(cord_num+1)*cord_num], \
            LDC);
}
}

return info;
}
#endif

#ifdef HAVE_CUBLAS_OFFLOAD
qmckl_exit_code qmckl_compute_dtmp_c_cublas_offload (
const qmckl_context context,
const int64_t cord_num,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t walk_num,
const double* een_rescaled_e_deriv_e,
const double* een_rescaled_n,
double* const dtmp_c );
#endif

Test

assert(qmckl_electron_provided(context));

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

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

assert(fabs(tmp_c[0][0][1][0][0] - 2.7083473948352403) < 1e-12);

assert(fabs(dtmp_c[0][1][0][0][0][0] - 0.237440520852232) < 1e-12);

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.

Get

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

Compute naive

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

integer function qmckl_compute_factor_een_naive_f( &
 context, walk_num, elec_num, nucl_num, cord_num,&
 dim_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_naive_f

qmckl_exit_code qmckl_compute_factor_een_naive (
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 );

Compute

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

integer function qmckl_compute_factor_een_f( &
 context, walk_num, elec_num, nucl_num, cord_num,   &
 dim_cord_vect, cord_vect_full, lkpm_combined_index, &
 tmp_c, een_rescaled_n, factor_een) &
 result(info)
use qmckl
implicit none
integer(qmckl_context), intent(in)  :: context
integer*8             , intent(in)  :: walk_num, elec_num, cord_num, nucl_num, dim_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)  :: tmp_c(elec_num, nucl_num,0:cord_num, 0:cord_num-1,  walk_num)
double precision      , intent(in)  :: een_rescaled_n(elec_num, nucl_num, 0:cord_num, walk_num)
double precision      , intent(out) :: factor_een(walk_num)

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

info = QMCKL_SUCCESS

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

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

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

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

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

factor_een = 0.0d0

do nw =1, walk_num
do n = 1, dim_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)
  if(cn == 0.d0) cycle

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

end function qmckl_compute_factor_een_f

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

Test

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

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

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

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.

Get

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

Compute Naive

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

integer function qmckl_compute_factor_een_deriv_e_naive_f( &
 context, walk_num, elec_num, nucl_num, cord_num, dim_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_naive_f

qmckl_exit_code qmckl_compute_factor_een_deriv_e_naive (
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 );

Compute

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

integer function qmckl_compute_factor_een_deriv_e_f( &
 context, walk_num, elec_num, nucl_num, &
 cord_num, dim_cord_vect, cord_vect_full, lkpm_combined_index, &
 tmp_c, dtmp_c, een_rescaled_n, een_rescaled_n_deriv_e, factor_een_deriv_e)&
 result(info)
use qmckl
implicit none
integer(qmckl_context), intent(in)  :: context
integer*8             , intent(in)  :: walk_num, elec_num, cord_num, nucl_num, dim_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)  :: tmp_c(elec_num, nucl_num,0:cord_num, 0:cord_num-1,  walk_num)
double precision      , intent(in)  :: dtmp_c(elec_num, 4, nucl_num,0:cord_num, 0:cord_num-1,  walk_num)
double precision      , intent(in)  :: een_rescaled_n(elec_num, nucl_num, 0:cord_num, walk_num)
double precision      , intent(in)  :: een_rescaled_n_deriv_e(elec_num, 4, nucl_num, 0:cord_num, walk_num)
double precision      , intent(out) :: factor_een_deriv_e(elec_num,4,walk_num)

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

info = QMCKL_SUCCESS

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

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

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

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

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

factor_een_deriv_e = 0.0d0

do nw =1, walk_num
do n = 1, dim_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)
  if(cn == 0.d0) cycle

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

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

end function qmckl_compute_factor_een_deriv_e_f

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* tmp_c,
const double* dtmp_c,
const double* een_rescaled_n,
const double* een_rescaled_n_deriv_e,
double* const factor_een_deriv_e );

Test

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

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

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