TREX/qmckl

Fork 0

mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-11-09 07:33:59 +01:00

Anthony Scemama 1f3a08fa30 Avoid FPE in Jastrow GL

2023-08-23 14:32:33 +02:00

393 KiB

Raw Blame History

CHAMP Jastrow Factor

Introduction
Context
Computation

Introduction

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

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

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

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

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

$J_{\text{ee}}$ contains electron-electron terms: \[ J_{\text{ee}}(\mathbf{r}) = \sum_{i=1}^{N_\text{elec}} \sum_{j=1}^{i-1} \frac{\frac{1}{2}(1+\delta^{\uparrow\downarrow}_{ij}) b_1\, f_{\text{ee}}(r_{ij})}{1+b_2\, f_{\text{ee}}(r_{ij})} + \sum_{p=2}^{N_\text{ord}^b} b_{p+1}\, [f_{\text{ee}}(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_\text{ee}({r}_{ij}) \right]^k \left[ \left[ g_\alpha({R}_{i\alpha}) \right]^l + \left[ g_\alpha({R}_{j\alpha}) \right]^l \right] \left[ g_\alpha({R}_{i\,\alpha}) \, g_\alpha({R}_{j\alpha}) \right]^{(p-k-l)/2} \]

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

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

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

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

The eN and eeN parameters are the same of all identical nuclei. The types of nuclei use zero-based indexing.

Context

The following data stored in the context:

Variable	Type	Description
`uninitialized`	`int32_t`	Keeps bits set for uninitialized data
`rescale_factor_ee`	`double`	The distance scaling factor
`rescale_factor_en`	`double[type_nucl_num]`	The distance scaling factor
`aord_num`	`int64_t`	The number of a coeffecients
`bord_num`	`int64_t`	The number of b coeffecients
`cord_num`	`int64_t`	The number of c coeffecients
`type_nucl_num`	`int64_t`	Number of Nuclei types
`type_nucl_vector`	`int64_t[nucl_num]`	IDs of types of Nuclei. These use 0-based indexing as in C.
`a_vector`	`double[aord_num + 1][type_nucl_num]`	a polynomial coefficients
`b_vector`	`double[bord_num + 1]`	b polynomial coefficients
`c_vector`	`double[dim_c_vector][type_nucl_num]`	c polynomial coefficients

Computed data:

Variable	Type	In/Out
`dim_c_vector`	`int64_t`	Number of unique C coefficients
`dim_c_vector_date`	`uint64_t`	Number of unique C coefficients
`asymp_jasa`	`double[type_nucl_num]`	Asymptotic component
`asymp_jasa_date`	`uint64_t`	Ladt modification of the asymptotic component
`asymp_jasb`	`double[2]`	Asymptotic component (up- or down-spin)
`asymp_jasb_date`	`uint64_t`	Ladt modification of the asymptotic component
`c_vector_full`	`double[dim_c_vector][nucl_num]`	vector of non-zero coefficients
`c_vector_full_date`	`uint64_t`	Keep track of changes here
`lkpm_combined_index`	`int64_t[4][dim_c_vector]`	Transform l,k,p, and m into consecutive indices
`lkpm_combined_index_date`	`uint64_t`	Transform l,k,p, and m into consecutive indices
`tmp_c`	`double[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
`ee_distance_rescaled`	`double[walk_num][num][num]`	Electron-electron rescaled distances
`ee_distance_rescaled_date`	`uint64_t`	Last modification date of the electron-electron distances
`ee_distance_rescaled_gl`	`double[walk_num][4][num][num]`	Electron-electron rescaled distances derivatives
`ee_distance_rescaled_gl_date`	`uint64_t`	Last modification date of the electron-electron distance derivatives
`en_distance_rescaled`	`double[walk_num][nucl_num][num]`	Electron-nucleus distances
`en_distance_rescaled_date`	`uint64_t`	Last modification date of the electron-electron distances
`en_distance_rescaled_gl`	`double[walk_num][4][nucl_num][num]`	Electron-electron rescaled distances derivatives
`en_distance_rescaled_gl_date`	`uint64_t`	Last modification date of the electron-electron distance derivatives
`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_gl`	`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_gl_date`	`uint64_t`	Keep track of the date of creation
`een_rescaled_n_gl`	`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_gl_date`	`uint64_t`	Keep track of the date of creation
`factor_ee`	`double[walk_num]`	Jastrow factor: electron-electron part
`factor_ee_date`	`uint64_t`	Jastrow factor: electron-electron part
`factor_en`	`double[walk_num]`	Jastrow factor: electron-nucleus part
`factor_en_date`	`uint64_t`	Jastrow factor: electron-nucleus part
`factor_een`	`double[walk_num]`	Jastrow factor: electron-electron-nucleus part
`factor_een_date`	`uint64_t`	Jastrow factor: electron-electron-nucleus part
`factor_ee_gl`	`double[walk_num][4][elec_num]`	Derivative of the Jastrow factor: electron-electron-nucleus part
`factor_ee_gl_date`	`uint64_t`	Keep track of the date for the derivative
`factor_en_gl`	`double[walk_num][4][elec_num]`	Derivative of the Jastrow factor: electron-electron-nucleus part
`factor_en_gl_date`	`uint64_t`	Keep track of the date for the en derivative
`factor_een_gl`	`double[walk_num][4][elec_num]`	Derivative of the Jastrow factor: electron-electron-nucleus part
`factor_een_gl_date`	`uint64_t`	Keep track of the date for the een derivative
`value`	`double[walk_num]`	Value of the Jastrow factor
`value_date`	`uint64_t`	Keep track of the date
`gl`	`double[walk_num][4][elec_num]`	Gradient and Laplacian of the Jastrow factor
`value_date`	`uint64_t`	Keep track of the date

Data structure

typedef struct qmckl_jastrow_champ_struct{
int64_t * restrict lkpm_combined_index;
int64_t * restrict type_nucl_vector;
double  * restrict asymp_jasa;
double             asymp_jasb[2];
double  * restrict a_vector;
double  * restrict b_vector;
double  * restrict c_vector;
double  * restrict c_vector_full;
double  * restrict dtmp_c;
double  * restrict ee_distance_rescaled;
double  * restrict ee_distance_rescaled_gl;
double  * restrict een_rescaled_e;
double  * restrict een_rescaled_e_gl;
double  * restrict een_rescaled_n;
double  * restrict een_rescaled_n_gl;
double  * restrict en_distance_rescaled;
double  * restrict en_distance_rescaled_gl;
double  * restrict factor_ee;
double  * restrict factor_ee_gl;
double  * restrict factor_een;
double  * restrict factor_een_gl;
double  * restrict factor_en;
double  * restrict factor_en_gl;
double  * restrict rescale_factor_en;
double  * restrict tmp_c;
double  * restrict value;
double  * restrict gl;
int64_t   aord_num;
int64_t   bord_num;
int64_t   cord_num;
int64_t   dim_c_vector;
int64_t   type_nucl_num;
uint64_t  asymp_jasa_date;
uint64_t  asymp_jasb_date;
uint64_t  c_vector_full_date;
uint64_t  dim_c_vector_date;
uint64_t  dtmp_c_date;
uint64_t  ee_distance_rescaled_date;
uint64_t  ee_distance_rescaled_gl_date;
uint64_t  een_rescaled_e_date;
uint64_t  een_rescaled_e_gl_date;
uint64_t  een_rescaled_n_date;
uint64_t  een_rescaled_n_gl_date;
uint64_t  en_distance_rescaled_date;
uint64_t  en_distance_rescaled_gl_date;
uint64_t  factor_ee_date;
uint64_t  factor_ee_gl_date;
uint64_t  factor_een_date;
uint64_t  factor_een_gl_date;
uint64_t  factor_en_date;
uint64_t  factor_en_gl_date;
uint64_t  lkpm_combined_index_date;
uint64_t  tmp_c_date;
uint64_t  value_date;
uint64_t  gl_date;
double    rescale_factor_ee;
int32_t   uninitialized;
bool      provided;

} qmckl_jastrow_champ_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_champ(qmckl_context context);

qmckl_exit_code qmckl_init_jastrow_champ(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_champ.uninitialized = (1 << 10) - 1;

/* Default values */
ctx->jastrow_champ.aord_num = -1;

ctx->jastrow_champ.bord_num = -1;

ctx->jastrow_champ.dim_c_vector = -1;
ctx->jastrow_champ.cord_num = -1;

ctx->jastrow_champ.type_nucl_num = -1;

return QMCKL_SUCCESS;
}

Initialization functions

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

qmckl_exit_code  qmckl_set_jastrow_champ_rescale_factor_ee (qmckl_context context, const double  kappa_ee);
qmckl_exit_code  qmckl_set_jastrow_champ_rescale_factor_en (qmckl_context context, const double* kappa_en, const int64_t size_max);
qmckl_exit_code  qmckl_set_jastrow_champ_aord_num          (qmckl_context context, const int64_t aord_num);
qmckl_exit_code  qmckl_set_jastrow_champ_bord_num          (qmckl_context context, const int64_t bord_num);
qmckl_exit_code  qmckl_set_jastrow_champ_cord_num          (qmckl_context context, const int64_t cord_num);
qmckl_exit_code  qmckl_set_jastrow_champ_type_nucl_num     (qmckl_context context, const int64_t type_nucl_num);
qmckl_exit_code  qmckl_set_jastrow_champ_type_nucl_vector  (qmckl_context context, const int64_t* type_nucl_vector, const int64_t nucl_num);
qmckl_exit_code  qmckl_set_jastrow_champ_a_vector          (qmckl_context context, const double * a_vector, const int64_t size_max);
qmckl_exit_code  qmckl_set_jastrow_champ_b_vector          (qmckl_context context, const double * b_vector, const int64_t size_max);
qmckl_exit_code  qmckl_set_jastrow_champ_c_vector          (qmckl_context context, const double * c_vector, const int64_t size_max);

#+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.

Fortran interface

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

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

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

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

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

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

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

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

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

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

end interface

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_champ_provided           (const qmckl_context context);

Fortran interface

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

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

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

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

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

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

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

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

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

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

end interface

Test

/* Reference input data */
int64_t walk_num      = n2_walk_num;
int64_t elec_num      = n2_elec_num;
int64_t elec_up_num   = n2_elec_up_num;
int64_t elec_dn_num   = n2_elec_dn_num;
int64_t nucl_num      = n2_nucl_num;
double  rescale_factor_ee   = 0.6;
double  rescale_factor_en[2] = { 0.6, 0.6 };
double* elec_coord    = &(n2_elec_coord[0][0][0]);

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

/* Provide Electron data */

qmckl_exit_code rc;

assert(!qmckl_electron_provided(context));

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

assert(qmckl_electron_provided(context));

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

double elec_coord2[walk_num*3*elec_num];

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


/* Provide Nucleus data */

assert(!qmckl_nucleus_provided(context));

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

double nucl_coord2[3*nucl_num];

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

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

assert(!qmckl_nucleus_provided(context));

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

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

double nucl_charge2[nucl_num];

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

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

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

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.

Electron-electron component

Asymptotic component

Calculate the asymptotic component asymp_jasb to be substracted from the electron-electron jastrow factor $J_{\text{ee}}$. Two values are computed. The first one is for parallel spin pairs, and the second one for antiparallel spin pairs.

\[ J_{\text{ee}}^{\infty} = \frac{\frac{1}{2}(1+\delta^{\uparrow \downarrow})\,b_1 \kappa_\text{ee}^{-1}}{1 + b_2\, \kappa_\text{ee}^{-1}} + \sum_{p=2}^{N_\text{ord}^b} b_{p+1}\, \kappa_\text{ee}^{-p} \]

Get

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

Fortran interface

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

Compute

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

integer function qmckl_compute_jastrow_champ_asymp_jasb_doc_f(context, bord_num, b_vector, rescale_factor_ee, asymp_jasb) &
result(info)
use qmckl
implicit none
integer(qmckl_context), intent(in)  :: context
integer*8             , intent(in)  :: bord_num
double precision      , intent(in)  :: b_vector(bord_num + 1)
double precision      , intent(in)  :: rescale_factor_ee
double precision      , intent(out) :: asymp_jasb(2)

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

info = QMCKL_SUCCESS

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

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

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

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

end function qmckl_compute_jastrow_champ_asymp_jasb_doc_f

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

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

if (bord_num < 0) {
return QMCKL_INVALID_ARG_2;
}

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

double f = 0.;
double x = kappa_inv;
for (int k = 2; k <= bord_num; ++k) {
x *= kappa_inv;
f += b_vector[k]*x;
}

asymp_jasb[0] = 0.5 * asym_one + f;
asymp_jasb[1] = asym_one + f;

return QMCKL_SUCCESS;
}

qmckl_exit_code qmckl_compute_jastrow_champ_asymp_jasb (
  const qmckl_context context,
  const int64_t bord_num,
  const double* b_vector,
  const double rescale_factor_ee,
  double* const asymp_jasb )
{
#ifdef HAVE_HPC
 return qmckl_compute_jastrow_champ_asymp_jasb_hpc (context,
  bord_num, b_vector, rescale_factor_ee, asymp_jasb);
#else
 return qmckl_compute_jastrow_champ_asymp_jasb_doc (context,
  bord_num, b_vector, rescale_factor_ee, asymp_jasb);
#endif
}

Test

assert(qmckl_electron_provided(context));

int64_t type_nucl_num = n2_type_nucl_num;
int64_t* type_nucl_vector = &(n2_type_nucl_vector[0]);
int64_t aord_num = n2_aord_num;
int64_t bord_num = n2_bord_num;
int64_t cord_num = n2_cord_num;
double* a_vector = &(n2_a_vector[0][0]);
double* b_vector = &(n2_b_vector[0]);
double* c_vector = &(n2_c_vector[0][0]);
int64_t dim_c_vector=0;

assert(!qmckl_jastrow_champ_provided(context));

/* Set the data */
rc = qmckl_check(context,
           qmckl_set_jastrow_champ_aord_num(context, aord_num)
           );
rc = qmckl_check(context,
           qmckl_set_jastrow_champ_bord_num(context, bord_num)
           );
rc = qmckl_check(context,
           qmckl_set_jastrow_champ_cord_num(context, cord_num)
           );
assert(rc == QMCKL_SUCCESS);
rc = qmckl_check(context,
           qmckl_set_jastrow_champ_type_nucl_num(context, type_nucl_num)
           );
assert(rc == QMCKL_SUCCESS);
rc = qmckl_check(context,
           qmckl_set_jastrow_champ_type_nucl_vector(context, type_nucl_vector, nucl_num)
           );
assert(rc == QMCKL_SUCCESS);
rc = qmckl_check(context,
           qmckl_set_jastrow_champ_a_vector(context, a_vector,(aord_num+1)*type_nucl_num)
           );
assert(rc == QMCKL_SUCCESS);
rc = qmckl_check(context,
           qmckl_set_jastrow_champ_b_vector(context, b_vector,(bord_num+1))
           );
assert(rc == QMCKL_SUCCESS);
rc = qmckl_check(context,
           qmckl_get_jastrow_champ_dim_c_vector(context, &dim_c_vector)
           );
assert(rc == QMCKL_SUCCESS);
rc = qmckl_check(context,
           qmckl_set_jastrow_champ_c_vector(context, c_vector,dim_c_vector*type_nucl_num)
           );
assert(rc == QMCKL_SUCCESS);

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

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

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

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

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

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

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

Electron-electron component

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

\[ f_\text{ee} = \sum_{i,j<i} \left[ \frac{\delta_{ij}^{\uparrow\downarrow} B_0\, C_{ij}}{1 + B_1\, C_{ij}} + \sum_{k=2}^{n_\text{ord}} B_k\, C_{ij}^k - {J_{\text{ee}}^{\infty}}_{ij} \right] \]

$\delta$ is the spin factor, $B$ is the vector of $b$ parameters, $C$ is the array of rescaled distances.

$f_{\text{ee}}$ can be rewritten as:

\[ f_\text{ee} = \frac{1}{2} \left[ \sum_{i,j} \frac{\delta_{ij}^{\uparrow\downarrow} B_0\, C_{ij}}{1 + B_1\, C_{ij}} + \sum_{i,j} \sum_{k=2}^{n_\text{ord}} B_k\, C_{ij}^k \right] - \left[ \frac{n_\uparrow (n_\uparrow-1) + n_\downarrow (n_\downarrow-1)}{2}\, J_{\text{ee}}^{\infty}}_{\uparrow \uparrow} + n_\uparrow\,n_\downarrow\, J_{\text{ee}}^{\infty}}_{\uparrow \downarrow} \right] \]

Get

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

Fortran interface

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

Compute

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

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

integer*8 :: i, j, k, nw
double precision   :: x, xk

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


do nw =1, walk_num

factor_ee(nw) = 0.0d0
do j=1,elec_num
   do i=1,j-1
      x = ee_distance_rescaled(i,j,nw)
      if ( (j <= up_num).or.(i > up_num) ) then
         factor_ee(nw) = factor_ee(nw) + 0.5d0 * b_vector(1) * x / (1.d0 + b_vector(2) * x) - asymp_jasb(1)
      else
         factor_ee(nw) = factor_ee(nw) + b_vector(1) * x / (1.d0 + b_vector(2) * x) - asymp_jasb(2)
      endif

      xk = x
      do k=2,bord_num
         xk = xk * x
         factor_ee(nw) = factor_ee(nw) + b_vector(k+1)* xk
      end do
   end do
end do

end do

end function qmckl_compute_jastrow_champ_factor_ee_doc_f

qmckl_exit_code qmckl_compute_jastrow_champ_factor_ee_hpc (
     const qmckl_context context,
     const int64_t walk_num,
     const int64_t elec_num,
     const int64_t up_num,
     const int64_t bord_num,
     const double* b_vector,
     const double* ee_distance_rescaled,
     const double* asymp_jasb,
     double* const factor_ee ) {

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

 const int64_t dn_num = elec_num - up_num;
 const double fshift = 0.5 * (double) ((dn_num-1)*dn_num + (up_num-1)*up_num) * asymp_jasb[0] +
     (float) (up_num*dn_num) * asymp_jasb[1];
 for (int nw = 0; nw < walk_num; ++nw) {
   factor_ee[nw] = 0.;

   size_t ishift = nw * elec_num * elec_num;
   for (int j = 0; j < up_num; ++j ) {
     const double* xj = &(ee_distance_rescaled[j * elec_num + ishift]);
     for (int i = 0; i < j ; ++i) {
       factor_ee[nw] += 0.5 * b_vector[0]*xj[i] / (1. + b_vector[1]*xj[i]);
     }
   }

   for (int j = up_num ; j < elec_num; ++j ) {
     const double* xj = &(ee_distance_rescaled[j * elec_num + ishift]);
     for (int i = 0; i < up_num; ++i) {
       factor_ee[nw] += b_vector[0]*xj[i] / (1. + b_vector[1]*xj[i]);
     }
     for (int i = up_num ; i < j ; ++i) {
       factor_ee[nw] += 0.5 * b_vector[0]*xj[i] / (1. + b_vector[1]*xj[i]);
     }

   }

   factor_ee[nw] -= fshift;

   for (int j=0; j < elec_num; ++j ) {
     const double* xj = &(ee_distance_rescaled[j * elec_num + ishift]);
     for (int i=0; i < j ; ++i) {
       const double x = xj[i];
       double xk = x;
       for (int k = 2; k <= bord_num; ++k) {
         xk *= x;
         factor_ee[nw] += b_vector[k] * xk;
       }

     }
   }
 }

 return QMCKL_SUCCESS;
}

qmckl_exit_code
qmckl_compute_jastrow_champ_factor_ee (const qmckl_context context,
                        const int64_t walk_num,
                        const int64_t elec_num,
                        const int64_t up_num,
                        const int64_t bord_num,
                        const double* b_vector,
                        const double* ee_distance_rescaled,
                        const double* asymp_jasb,
                        double* const factor_ee )
{

#ifdef HAVE_HPC
 return qmckl_compute_jastrow_champ_factor_ee_hpc(context, walk_num, elec_num,
     up_num, bord_num, b_vector, ee_distance_rescaled, asymp_jasb,
     factor_ee);
#else
 return qmckl_compute_jastrow_champ_factor_ee_doc(context, walk_num, elec_num,
     up_num, bord_num, b_vector, ee_distance_rescaled, asymp_jasb,
     factor_ee);
#endif
}

Test

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

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

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

Derivative

The derivative of factor_ee is computed using the ee_distance_rescaled and the electron-electron rescaled distances derivatives ee_distance_rescaled_gl. There are four components, the gradient which has 3 components in the $x, y, z$ directions and the laplacian as the last component.

\[ \nabla_i f_\text{ee} = \sum_{j\ne i} \left[\frac{\delta_{ij}^{\uparrow\downarrow} B_0\, \nabla_i C_{ij}}{(1 + B_1\, C_{ij})^2} + \sum^{n_\text{ord}}_{k=2} B_k\, k\, C_{ij}^{k-1} \nabla C_{ij} \right] \]

\[ \Delta_i f_\text{ee} = \sum_{j \ne i} \left[ \delta_{ij}^{\uparrow\downarrow} B_0 \left(\frac{ \Delta_i C_{ij}}{(1 + B_1\, C_{ij})^2} -\frac{2\,B_1 \left(\nabla_i C_{ij}\right)^2 }{(1 + B_1\, C_{ij})^3} \right) + \sum^{n_\text{ord}}_{k=2} B_k\, k\, \left((k-1)\, C_{ij}^{k-2} \left(\nabla_i C_{ij}\right)^2 + C_{ij}^{k-1} \Delta_i C_{ij} \right) \right] \]

Get

qmckl_exit_code
qmckl_get_jastrow_champ_factor_ee_gl(qmckl_context context,
                               double* const factor_ee_gl,
                               const int64_t size_max);

Fortran interface

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

Compute

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

integer function qmckl_compute_jastrow_champ_factor_ee_gl_doc_f( &
context, walk_num, elec_num, up_num, bord_num, &
b_vector, ee_distance_rescaled, ee_distance_rescaled_gl,  &
factor_ee_gl) &
result(info)
use qmckl
implicit none
integer(qmckl_context), intent(in)  :: context
integer*8             , intent(in)  :: walk_num, elec_num, bord_num, up_num
double precision      , intent(in)  :: b_vector(bord_num + 1)
double precision      , intent(in)  :: ee_distance_rescaled(elec_num, elec_num,walk_num)
double precision      , intent(in)  :: ee_distance_rescaled_gl(4,elec_num, elec_num,walk_num)   !TODO
double precision      , intent(out) :: factor_ee_gl(elec_num,4,walk_num)

integer*8 :: i, j, k, nw, ii
double precision   :: x, x1, kf
double precision   :: denom, invdenom, invdenom2, f
double precision   :: grad_c2
double precision   :: dx(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 (bord_num < 0) then
info = QMCKL_INVALID_ARG_4
return
endif

do nw =1, walk_num
factor_ee_gl(:,:,nw) = 0.0d0

do j = 1, elec_num
   do i = 1, elec_num
      if (i == j) cycle

      x = ee_distance_rescaled(i,j,nw)

      denom         = 1.0d0 + b_vector(2) * x
      invdenom      = 1.0d0 / denom
      invdenom2     = invdenom * invdenom

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

      grad_c2 = dx(1)*dx(1) + dx(2)*dx(2) + dx(3)*dx(3)

      if((i <= up_num .and. j <= up_num ) .or. (i >  up_num .and. j >  up_num)) then
         f = 0.5d0 * b_vector(1) * invdenom2
      else
         f = b_vector(1) * invdenom2
      end if

      factor_ee_gl(i,1,nw) = factor_ee_gl(i,1,nw) + f * dx(1)
      factor_ee_gl(i,2,nw) = factor_ee_gl(i,2,nw) + f * dx(2)
      factor_ee_gl(i,3,nw) = factor_ee_gl(i,3,nw) + f * dx(3)
      factor_ee_gl(i,4,nw) = factor_ee_gl(i,4,nw) &
           + f * (dx(4) - 2.d0 * b_vector(2) * grad_c2 * invdenom)


      kf = 2.d0
      x1 = x
      x = 1.d0
      do k=2, bord_num
         f = b_vector(k+1) * kf * x
         factor_ee_gl(i,1,nw) = factor_ee_gl(i,1,nw) + f * x1 * dx(1)
         factor_ee_gl(i,2,nw) = factor_ee_gl(i,2,nw) + f * x1 * dx(2)
         factor_ee_gl(i,3,nw) = factor_ee_gl(i,3,nw) + f * x1 * dx(3)
         factor_ee_gl(i,4,nw) = factor_ee_gl(i,4,nw) &
              + f * (x1 * dx(4) + (kf-1.d0) * grad_c2)
         x = x*x1
         kf = kf + 1.d0
      end do

   end do
end do

end do

end function qmckl_compute_jastrow_champ_factor_ee_gl_doc_f

qmckl_exit_code qmckl_compute_jastrow_champ_factor_ee_gl_hpc(
 const qmckl_context context,
 const int64_t walk_num,
 const int64_t elec_num,
 const int64_t up_num,
 const int64_t bord_num,
 const double* b_vector,
 const double* ee_distance_rescaled,
 const double* ee_distance_rescaled_gl,
 double* const factor_ee_gl ) {

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


memset(factor_ee_gl, 0, elec_num*4*walk_num*sizeof(double));

double kf[bord_num+1];
for (int k=0 ; k<=bord_num ; ++k) {
kf[k] = (double) k;
}

for (int nw = 0; nw < walk_num; ++nw) {
for (int j = 0; j < elec_num; ++j) {
 const double* dxj = &ee_distance_rescaled_gl[4*elec_num*(j+nw*elec_num)];
 const double*  xj = &ee_distance_rescaled   [  elec_num*(j+nw*elec_num)];

 double * restrict factor_ee_gl_0 = &(factor_ee_gl[nw*elec_num*4]);
 double * restrict factor_ee_gl_1 = factor_ee_gl_0 + elec_num;
 double * restrict factor_ee_gl_2 = factor_ee_gl_1 + elec_num;
 double * restrict factor_ee_gl_3 = factor_ee_gl_2 + elec_num;

 for (int i = 0; i < elec_num; ++i) {
   if (j == i) continue;

   double x = xj[i];

   const double denom      = 1.0 + b_vector[1]*x;
   const double invdenom   = 1.0  / denom;
   const double invdenom2  = invdenom * invdenom;

   const double* restrict dx = dxj + 4*i;

   const double grad_c2 = dx[0]*dx[0] + dx[1]*dx[1] + dx[2]*dx[2];

   double f =
    (i < up_num && j < up_num ) || (i >= up_num && j >= up_num) ?
    0.5 * b_vector[0] * invdenom2 : b_vector[0] * invdenom2;

   factor_ee_gl_0[i] += f*dx[0];
   factor_ee_gl_1[i] += f*dx[1];
   factor_ee_gl_2[i] += f*dx[2];
   factor_ee_gl_3[i] += f*dx[3];
   factor_ee_gl_3[i] -= f*grad_c2*invdenom*2.0 * b_vector[1];


   double xk[bord_num+1];
   xk[0] = 1.0;
   for (int k=1 ; k<= bord_num ; ++k) {
     xk[k] = xk[k-1]*x;
   }

   for (int k=2 ; k<= bord_num ; ++k) {
     const double f1 = b_vector[k] * kf[k] * xk[k-2];
     const double f2 = f1*xk[1];
     factor_ee_gl_0[i] += f2*dx[0];
     factor_ee_gl_1[i] += f2*dx[1];
     factor_ee_gl_2[i] += f2*dx[2];
     factor_ee_gl_3[i] += f2*dx[3];
     factor_ee_gl_3[i] += f1*kf[k-1]*grad_c2;
   }
 }
}
}

return QMCKL_SUCCESS;
}

qmckl_exit_code qmckl_compute_jastrow_champ_factor_ee_gl_hpc (
 const qmckl_context context,
 const int64_t walk_num,
 const int64_t elec_num,
 const int64_t up_num,
 const int64_t bord_num,
 const double* b_vector,
 const double* ee_distance_rescaled,
 const double* ee_distance_rescaled_gl,
 double* const factor_ee_gl );

qmckl_exit_code qmckl_compute_jastrow_champ_factor_ee_gl_doc (
 const qmckl_context context,
 const int64_t walk_num,
 const int64_t elec_num,
 const int64_t up_num,
 const int64_t bord_num,
 const double* b_vector,
 const double* ee_distance_rescaled,
 const double* ee_distance_rescaled_gl,
 double* const factor_ee_gl );

qmckl_exit_code qmckl_compute_jastrow_champ_factor_ee_gl (
 const qmckl_context context,
 const int64_t walk_num,
 const int64_t elec_num,
 const int64_t up_num,
 const int64_t bord_num,
 const double* b_vector,
 const double* ee_distance_rescaled,
 const double* ee_distance_rescaled_gl,
 double* const factor_ee_gl ) {

 #ifdef HAVE_HPC
 return qmckl_compute_jastrow_champ_factor_ee_gl_hpc(context, walk_num, elec_num, up_num, bord_num, b_vector, ee_distance_rescaled, ee_distance_rescaled_gl, factor_ee_gl );
 #else
 return qmckl_compute_jastrow_champ_factor_ee_gl_doc(context, walk_num, elec_num, up_num, bord_num, b_vector, ee_distance_rescaled, ee_distance_rescaled_gl, factor_ee_gl );
 #endif
}

Test

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

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

// check factor_ee_gl
printf("%f  %f\n", factor_ee_gl[0][0][0], -0.39319353942687446);
assert(fabs(factor_ee_gl[0][0][0]+0.39319353942687446) < 1.e-12);

printf("%f  %f\n", factor_ee_gl[0][1][0], 1.0535615450668214);
assert(fabs(factor_ee_gl[0][1][0]-1.0535615450668214) < 1.e-12);

printf("%f  %f\n", factor_ee_gl[0][2][0],-0.39098406960784515);
assert(fabs(factor_ee_gl[0][2][0]+0.39098406960784515) < 1.e-12);

printf("%f  %f\n", factor_ee_gl[0][3][0],2.8650469630854483);
assert(fabs(factor_ee_gl[0][3][0]-2.8650469630854483) < 1.e-12);

Electron-electron rescaled distances

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

\[ C_{ij} = \frac{ 1 - e^{-\kappa r_{ij}}}{\kappa} \]

where $r_{ij}$ is the matrix of electron-electron distances.

Get

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

Compute

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

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

integer*8 :: k

info = QMCKL_SUCCESS

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

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

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

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

end function qmckl_compute_ee_distance_rescaled_f

Test

assert(qmckl_electron_provided(context));


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

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

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

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

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

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

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

Electron-electron rescaled distance gradients and Laplacian with respect to electron coordinates

The rescaled distances, represented by $C_{ij} = (1 - e^{-\kappa_\text{e} r_{ij}})/\kappa_\text{e}$ are differentiated with respect to the electron coordinates. This information is stored in the tensor ee_distance_rescaled_gl. The initial three sequential elements of this three-dimensional tensor provide the $x$, $y$, and $z$ direction derivatives, while the fourth index corresponds to the Laplacian.

Get

qmckl_exit_code qmckl_get_jastrow_champ_ee_distance_rescaled_gl(qmckl_context context, double* const distance_rescaled_gl);

Compute

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

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

integer*8 :: k

info = QMCKL_SUCCESS

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

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

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

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

end function qmckl_compute_ee_distance_rescaled_gl_f

Test

assert(qmckl_electron_provided(context));


double ee_distance_rescaled_gl[4 * walk_num * elec_num * elec_num];
rc = qmckl_get_jastrow_champ_ee_distance_rescaled_gl(context, ee_distance_rescaled_gl);

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

Electron-nucleus component

Asymptotic component for

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

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

Get

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

Fortran interface

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

Compute

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

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

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

info = QMCKL_SUCCESS

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

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

do i=1,type_nucl_num

kappa_inv = 1.0d0 / rescale_factor_en(i)

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

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

end do

end function qmckl_compute_jastrow_champ_asymp_jasa_f

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

Test

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

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

#+end_src

Electron-nucleus component

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

\[ f_{\alpha}(R_{i\alpha}) = - \sum_{i,j<i} \left[ \frac{ A_0 C_{ij}}{1 - A_1 C_{ij}} + \sum^{N^\alpha_{\text{ord}}}_{k}A_k C_{ij}^k \right] \]

Get

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

Fortran interface

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

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 nuclei
`type_nucl_num`	`int64_t`	in	Number of unique nuclei
`type_nucl_vector`	`int64_t[nucl_num]`	in	IDs of unique nuclei
`aord_num`	`int64_t`	in	Number of coefficients
`a_vector`	`double[aord_num+1][type_nucl_num]`	in	List of coefficients
`en_distance_rescaled`	`double[walk_num][nucl_num][elec_num]`	in	Electron-nucleus distances
`asymp_jasa`	`double[type_nucl_num]`	in	Type of nuclei
`factor_en`	`double[walk_num]`	out	Electron-nucleus jastrow

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

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

info = QMCKL_SUCCESS

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

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

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

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

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

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


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

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

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

   end do
end do
end do

end function qmckl_compute_jastrow_champ_factor_en_doc_f

qmckl_exit_code qmckl_compute_jastrow_champ_factor_en (
     const qmckl_context context,
     const int64_t walk_num,
     const int64_t elec_num,
     const int64_t nucl_num,
     const int64_t type_nucl_num,
     const int64_t* type_nucl_vector,
     const int64_t aord_num,
     const double* a_vector,
     const double* en_distance_rescaled,
     const double* asymp_jasa,
     double* const factor_en );

qmckl_exit_code qmckl_compute_jastrow_champ_factor_en_doc (
     const qmckl_context context,
     const int64_t walk_num,
     const int64_t elec_num,
     const int64_t nucl_num,
     const int64_t type_nucl_num,
     const int64_t* type_nucl_vector,
     const int64_t aord_num,
     const double* a_vector,
     const double* en_distance_rescaled,
     const double* asymp_jasa,
     double* const factor_en );

qmckl_exit_code qmckl_compute_jastrow_champ_factor_en (
     const qmckl_context context,
     const int64_t walk_num,
     const int64_t elec_num,
     const int64_t nucl_num,
     const int64_t type_nucl_num,
     const int64_t* type_nucl_vector,
     const int64_t aord_num,
     const double* a_vector,
     const double* en_distance_rescaled,
     const double* asymp_jasa,
     double* const factor_en )
{
#ifdef HAVE_HPC
return qmckl_compute_jastrow_champ_factor_en_doc (context, walk_num, elec_num, nucl_num, type_nucl_num,
                                  type_nucl_vector, aord_num, a_vector, en_distance_rescaled,
                                  asymp_jasa, factor_en );
#else
return qmckl_compute_jastrow_champ_factor_en_doc (context, walk_num, elec_num, nucl_num, type_nucl_num,
                                  type_nucl_vector, aord_num, a_vector, en_distance_rescaled,
                                  asymp_jasa, factor_en );
#endif
}

Test

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

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

// calculate factor_en
printf("%f %f\n", factor_en[0], 2.781375792083587);
assert(fabs(22.781375792083587 - factor_en[0]) < 1.e-12);

Derivative

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

TODO: write equations.

Get

qmckl_exit_code
qmckl_get_jastrow_champ_factor_en_gl(qmckl_context context,
                               double* const factor_en_gl,
                               const int64_t size_max);

Fortran interface

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

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 nuclei
`type_nucl_num`	`int64_t`	in	Number of unique nuclei
`type_nucl_vector`	`int64_t[nucl_num]`	in	IDs of unique nuclei
`aord_num`	`int64_t`	in	Number of coefficients
`a_vector`	`double[aord_num+1][type_nucl_num]`	in	List of coefficients
`en_distance_rescaled`	`double[walk_num][nucl_num][elec_num]`	in	Electron-nucleus distances
`en_distance_rescaled_gl`	`double[walk_num][4][nucl_num][elec_num]`	in	Electron-nucleus distance derivatives
`factor_en_gl`	`double[walk_num][4][elec_num]`	out	Electron-nucleus jastrow

integer function qmckl_compute_jastrow_champ_factor_en_gl_f( &
context, walk_num, elec_num, nucl_num, type_nucl_num, &
type_nucl_vector, aord_num, a_vector, &
en_distance_rescaled, en_distance_rescaled_gl, factor_en_gl) &
result(info)
use qmckl
implicit none
integer(qmckl_context), intent(in)  :: context
integer*8             , intent(in)  :: walk_num, elec_num, aord_num, nucl_num, type_nucl_num
integer*8             , intent(in)  :: type_nucl_vector(nucl_num)
double precision      , intent(in)  :: a_vector(aord_num + 1, type_nucl_num)
double precision      , intent(in)  :: en_distance_rescaled(elec_num, nucl_num, walk_num)
double precision      , intent(in)  :: en_distance_rescaled_gl(4, elec_num, nucl_num, walk_num)
double precision      , intent(out) :: factor_en_gl(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_gl = 0.0d0
third = 1.0d0 / 3.0d0

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

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

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

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

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

   end do

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

end do
end do
end do

end function qmckl_compute_jastrow_champ_factor_en_gl_f

Test

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

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

// check factor_en_gl
assert(fabs( 0.19656663796630847 - factor_en_gl[0][0][0]) < 1.e-12);
assert(fabs( -0.3945140890522283 - factor_en_gl[0][1][0]) < 1.e-12);
assert(fabs( 0.5082964671286118  - factor_en_gl[0][2][0]) < 1.e-12);
assert(fabs( -1.8409460670666289 - factor_en_gl[0][3][0]) < 1.e-12);

Electron-nucleus rescaled distances

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

\[ C_{i\alpha} = \frac{ 1 - e^{-\kappa_\alpha R_{i\alpha}}}{\kappa_\alpha} \]

where $R_{i\alpha}$ is the matrix of electron-nucleus distances.

Get

qmckl_exit_code qmckl_get_electron_en_distance_rescaled(qmckl_context context, double* distance_rescaled);

Compute

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

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

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

info = QMCKL_SUCCESS

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

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

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

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

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

end function qmckl_compute_en_distance_rescaled_f

Test

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

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

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

// (e,n,w) in Fortran notation
// (1,1,1)
assert(fabs(en_distance_rescaled[0][0][0] - 0.4942158656729477) < 1.e-12);
// (1,2,1)
assert(fabs(en_distance_rescaled[0][1][0] - 1.2464137498005765) < 1.e-12);
// (2,1,1)
assert(fabs(en_distance_rescaled[0][0][1] - 0.5248654474756858) < 1.e-12);
// (1,1,2)
assert(fabs(en_distance_rescaled[0][0][5] - 0.19529459944794733) < 1.e-12);
// (1,2,2)
assert(fabs(en_distance_rescaled[0][1][5] - 1.2091967687767369) < 1.e-12);
// (2,1,2)
assert(fabs(en_distance_rescaled[0][0][6] - 0.4726452953409436) < 1.e-12);

Electron-electron rescaled distance gradients and Laplacian with respect to electron coordinates

The rescaled distances, represented by $C_{i\alpha} = (1 - e^{-\kappa_\alpha R_{i\alpha}})/\kappa$ are differentiated with respect to the electron coordinates. This information is stored in the tensor en_distance_rescaled_gl. The initial three sequential elements of this three-index tensor provide the $x$, $y$, and $z$ direction derivatives, while the fourth index corresponds to the Laplacian.

Get

qmckl_exit_code qmckl_get_electron_en_distance_rescaled_gl(qmckl_context context, double* distance_rescaled_gl);

Compute

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

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

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

info = QMCKL_SUCCESS

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

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

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

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

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

end function qmckl_compute_en_distance_rescaled_gl_f

Test

assert(qmckl_electron_provided(context));

assert(qmckl_nucleus_provided(context));

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

rc = qmckl_check(context,
           qmckl_get_electron_en_distance_rescaled_gl(context, &(en_distance_rescaled_gl[0][0][0][0]))
           );
assert (rc == QMCKL_SUCCESS);

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

Electron-electron-nucleus component

Electron-electron rescaled distances in $J_\text{eeN}$

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[ \exp\left(-\kappa_\text{e}\, r_{ij}\right) \right]^p \]

where $r_{ij}$ is the matrix of electron-electron distances.

Get

qmckl_exit_code
qmckl_get_jastrow_champ_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_ee`	`double`	in	Factor to rescale ee distances
`ee_distance`	`double[walk_num][elec_num][elec_num]`	in	Electron-electron distances for each walker
`een_rescaled_e`	`double[walk_num][0:cord_num][elec_num][elec_num]`	out	Electron-electron rescaled distances for each walker

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

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

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

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


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


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

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

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

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

end do

end function qmckl_compute_een_rescaled_e_doc_f

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

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

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

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

if (cord_num < 0) {
return QMCKL_INVALID_ARG_4;
}

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

memset(een_rescaled_e,0,walk_num*(cord_num+1)*elec_num*elec_num*sizeof(double));

const size_t elec_pairs = (size_t) (elec_num * (elec_num - 1)) / 2;
const size_t len_een_ij = (size_t) elec_pairs * (cord_num + 1);

// number of elements for the een_rescaled_e_ij[N_e*(N_e-1)/2][cord+1]
// probably in C is better [cord+1, Ne*(Ne-1)/2]
// elec_pairs = (elec_num * (elec_num - 1)) / 2;
// len_een_ij = elec_pairs * (cord_num + 1);
const size_t e2 = elec_num*elec_num;

#ifdef HAVE_OPENMP
#pragma omp parallel for
#endif
for (size_t nw = 0; nw < (size_t) walk_num; ++nw) {

double* een_rescaled_e_ij = malloc(len_een_ij*sizeof(double));

memset(&(een_rescaled_e_ij[0]),0,len_een_ij*sizeof(double));
for (size_t kk = 0; kk < elec_pairs ; ++kk) {
 een_rescaled_e_ij[kk]= 1.0;
}

size_t kk = 0;
for (size_t i = 0; i < (size_t) elec_num; ++i) {
#ifdef HAVE_OPENMP
#pragma omp simd
#endif
 for (size_t j = 0; j < i; ++j) {
   een_rescaled_e_ij[j + kk + elec_pairs] = -rescale_factor_ee * ee_distance[j + i*elec_num + nw*e2];
 }
 kk += i;
}

#ifdef HAVE_OPENMP
#pragma omp simd
#endif
for (size_t k = elec_pairs; k < 2*elec_pairs; ++k) {
   een_rescaled_e_ij[k] = exp(een_rescaled_e_ij[k]);
}


for (size_t l = 2; l < (size_t) (cord_num+1); ++l) {
#ifdef HAVE_OPENMP
#pragma omp simd
#endif
 for (size_t k = 0; k < elec_pairs; ++k) {
   // een_rescaled_e_ij(k, l + 1) = een_rescaled_e_ij(k, l + 1 - 1) * een_rescaled_e_ij(k, 2)
   een_rescaled_e_ij[k+l*elec_pairs] = een_rescaled_e_ij[k + (l - 1)*elec_pairs] * \
     een_rescaled_e_ij[k + elec_pairs];
 }
}

double* const een_rescaled_e_ = &(een_rescaled_e[nw*(cord_num+1)*e2]);
// prepare the actual een table
#ifdef HAVE_OPENMP
#pragma omp simd
#endif
for (size_t i = 0; i < e2; ++i){
   een_rescaled_e_[i] = 1.0;
}

for ( size_t l = 1; l < (size_t) (cord_num+1); ++l) {
 double* x = een_rescaled_e_ij + l*elec_pairs;
 double* const een_rescaled_e__ = &(een_rescaled_e_[l*e2]);
 double* een_rescaled_e_i = een_rescaled_e__;
 for (size_t i = 0; i < (size_t) elec_num; ++i) {
   for (size_t j = 0; j < i; ++j) {
     een_rescaled_e_i[j] = *x;
     een_rescaled_e__[i + j*elec_num] = *x;
     x += 1;
   }
   een_rescaled_e_i += elec_num;
 }
}

double* const x0 = &(een_rescaled_e[nw*e2*(cord_num+1)]);
for (size_t l = 0; l < (size_t) (cord_num + 1); ++l) {
 double* x1 = &(x0[l*e2]);
 for (size_t j = 0; j < (size_t) elec_num; ++j) {
   ,*x1 = 0.0;
   x1 += 1+elec_num;
 }
}

free(een_rescaled_e_ij);
}

return QMCKL_SUCCESS;
}

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

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

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

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

Test

assert(qmckl_electron_provided(context));


double een_rescaled_e[walk_num][(cord_num + 1)][elec_num][elec_num];
rc = qmckl_get_jastrow_champ_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.2211015082992776 ) < 1.e-12);
assert(fabs(een_rescaled_e[0][1][0][3]- 0.2611178387068169 ) < 1.e-12);
assert(fabs(een_rescaled_e[0][1][0][4]- 0.0884012350763747 ) < 1.e-12);
assert(fabs(een_rescaled_e[0][2][1][3]- 0.1016685507354656 ) < 1.e-12);
assert(fabs(een_rescaled_e[0][2][1][4]- 0.0113118073246869 ) < 1.e-12);
assert(fabs(een_rescaled_e[0][2][1][5]- 0.5257156022077619 ) < 1.e-12);

Electron-electron rescaled distances derivatives in $J_\text{eeN}$

een_rescaled_e_gl 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_champ.

\[ \frac{\partial}{\partial x} \left[ {g_\text{e}(r)}\right]^p = -\frac{x}{r} \kappa_\text{e}\, p\,\left[ {g_\text{e}(r)}\right]^p \] \[ \Delta \left[ {g_\text{e}(r)}\right]^p = \frac{2}{r} \kappa_\text{e}\, p\,\left[ {g_\text{e}(r)}\right]^p \right] + \left(\frac{\partial}{\partial x}\left[ {g_\text{e}(r)}\right]^p \right)^2 + \left(\frac{\partial}{\partial y}\left[ {g_\text{e}(r)}\right]^p \right)^2 + \left(\frac{\partial}{\partial z}\left[ {g_\text{e}(r)}\right]^p \right)^2 \]

Get

qmckl_exit_code
qmckl_get_jastrow_champ_een_rescaled_e_gl(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_ee`	`double`	in	Factor to rescale ee distances
`coord_ee`	`double[walk_num][3][elec_num]`	in	Electron coordinates
`ee_distance`	`double[walk_num][elec_num][elec_num]`	in	Electron-electron distances
`een_rescaled_e`	`double[walk_num][0:cord_num][elec_num][elec_num]`	in	Electron-electron distances
`een_rescaled_e_gl`	`double[walk_num][0:cord_num][elec_num][4][elec_num]`	out	Electron-electron rescaled distances

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

allocate(elec_dist_gl(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_gl     = 0.0d0
do nw = 1, walk_num
do j = 1, elec_num
 do i = 1, j-1
   rij_inv = 1.0d0 / ee_distance(i, j, nw)
   do ii = 1, 3
     elec_dist_gl(ii, i, j) = (coord_ee(i, ii, nw) - coord_ee(j, ii, nw)) * rij_inv
   end do
   elec_dist_gl(4, i, j) = 2.0d0 * rij_inv
 end do

 elec_dist_gl(:, j, j) = 0.0d0

 do i = j+1, elec_num
   rij_inv = 1.0d0 / ee_distance(i, j, nw)
   do ii = 1, 3
     elec_dist_gl(ii, i, j) = (coord_ee(i, ii, nw) - coord_ee(j, ii, nw)) * rij_inv
   end do
   elec_dist_gl(4, i, j) = 2.0d0 * rij_inv
 end do
end do

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

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

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

end function qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl_f

Test

double een_rescaled_e_gl[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_champ_een_rescaled_e_gl(context,
    &(een_rescaled_e_gl[0][0][0][0][0]),size_max);

// value of (0,0,0,2,1)
assert(fabs(een_rescaled_e_gl[0][1][0][0][2] +  0.09831391870751387   ) < 1.e-12);
assert(fabs(een_rescaled_e_gl[0][1][0][0][3] +  0.017204157459682526  ) < 1.e-12);
assert(fabs(een_rescaled_e_gl[0][1][0][0][4] +  0.013345768421098641  ) < 1.e-12);
assert(fabs(een_rescaled_e_gl[0][2][1][0][3] +  0.03733086358273962   ) < 1.e-12);
assert(fabs(een_rescaled_e_gl[0][2][1][0][4] +  0.004922634822943517  ) < 1.e-12);
assert(fabs(een_rescaled_e_gl[0][2][1][0][5] +  0.5416751547830984    ) < 1.e-12);

Electron-nucleus rescaled distances in $J_\text{eeN}$

een_rescaled_n stores the table of the rescaled distances between electrons and nuclei raised to the power $p$ defined by cord_num:

\[ C_{i\alpha,p} = \left[ \exp\left(-\kappa_\alpha\, R_{i\alpha}\right) \right]^p \]

where $R_{i\alpha}$ is the matrix of electron-nucleus distances.

Get

qmckl_exit_code
qmckl_get_jastrow_champ_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
`type_nucl_num`	`int64_t`	in	Number of atom types
`type_nucl_vector`	`int64_t[nucl_num]`	in	Types of atoms
`cord_num`	`int64_t`	in	Order of polynomials
`rescale_factor_en`	`double[nucl_num]`	in	Factor to rescale ee distances
`en_distance`	`double[walk_num][elec_num][nucl_num]`	in	Electron-nucleus distances
`een_rescaled_n`	`double[walk_num][0:cord_num][nucl_num][elec_num]`	out	Electron-nucleus rescaled distances

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

info = QMCKL_SUCCESS

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

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

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

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

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

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

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

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

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

end do

end function qmckl_compute_een_rescaled_n_f

/*
qmckl_exit_code qmckl_compute_een_rescaled_n (
 const qmckl_context context,
 const int64_t walk_num,
 const int64_t elec_num,
 const int64_t nucl_num,
 const int64_t type_nucl_num,
 int64_t* const type_nucl_vector,
 const int64_t cord_num,
 const double* rescale_factor_en,
 const double* en_distance,
 double* const een_rescaled_n ) {


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

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

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

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

if (cord_num < 0) {
return QMCKL_INVALID_ARG_5;
}

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

for (int nw = 0; nw < walk_num; ++nw) {
for (int a = 0; a < nucl_num; ++a) {
 for (int i = 0; i < elec_num; ++i) {
   een_rescaled_n[i + a*elec_num + nw * elec_num*nucl_num*(cord_num+1)] = 1.0;
   een_rescaled_n[i + a*elec_num + elec_num*nucl_num + nw*elec_num*nucl_num*(cord_num+1)] =
     exp(-rescale_factor_en[type_nucl_vector[a]] * en_distance[a + i*nucl_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;
}
*/

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_champ_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.2603169838750542 )< 1.e-12);
assert(fabs(een_rescaled_n[0][1][0][3]-0.3016180139679065 )< 1.e-12);
assert(fabs(een_rescaled_n[0][1][0][4]-0.10506023826192266)< 1.e-12);
assert(fabs(een_rescaled_n[0][2][1][3]-0.9267719759374164 )< 1.e-12);
assert(fabs(een_rescaled_n[0][2][1][4]-0.11497585238132658)< 1.e-12);
assert(fabs(een_rescaled_n[0][2][1][5]-0.07534033469115217)< 1.e-12);

Electron-nucleus rescaled distances derivatives in $J_\text{eeN}$

een_rescaled_n_gl stores the table of the derivatives of the rescaled distances between all electron-nucleus pairs and raised to the power $p$ defined by cord_num. Here we take its derivatives required for the een jastrow_champ.

\[ \frac{\partial}{\partial x} \left[ {g_\alpha(R_{i\alpha})}\right]^p = -\frac{x}{R_{i\alpha}} \kappa_\alpha\, p\,\left[ {g_\alpha(R_{i\alpha})}\right]^p \] \[ \Delta \left[ {g_\alpha(R_{i\alpha})}\right]^p = \frac{2}{R_{i\alpha}} \kappa_\alpha\, p\,\left[ {g_\alpha(R_{i\alpha})}\right]^p \right] + \left(\frac{\partial}{\partial x}\left[ {g_\alpha(R_{i\alpha})}\right]^p \right)^2 + \left(\frac{\partial}{\partial y}\left[ {g_\alpha(R_{i\alpha})}\right]^p \right)^2 + \left(\frac{\partial}{\partial z}\left[ {g_\alpha(R_{i\alpha})}\right]^p \right)^2 \]

Get

qmckl_exit_code
qmckl_get_jastrow_champ_een_rescaled_n_gl(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
`type_nucl_num`	`int64_t`	in	Number of atom types
`type_nucl_vector`	`int64_t[nucl_num]`	in	Types of atoms
`cord_num`	`int64_t`	in	Order of polynomials
`rescale_factor_en`	`double[nucl_num]`	in	Factor to rescale ee distances
`coord_ee`	`double[walk_num][3][elec_num]`	in	Electron coordinates
`coord_en`	`double[3][nucl_num]`	in	Nuclear coordinates
`en_distance`	`double[walk_num][elec_num][nucl_num]`	in	Electron-nucleus distances
`een_rescaled_n`	`double[walk_num][0:cord_num][nucl_num][elec_num]`	in	Electron-nucleus distances
`een_rescaled_n_gl`	`double[walk_num][0:cord_num][nucl_num][4][elec_num]`	out	Electron-nucleus rescaled distances

integer function qmckl_compute_jastrow_champ_factor_een_rescaled_n_gl_f( &
context, walk_num, elec_num, nucl_num, type_nucl_num, type_nucl_vector, &
cord_num, rescale_factor_en, &
coord_ee, coord_en, en_distance, een_rescaled_n, een_rescaled_n_gl) &
result(info)
use qmckl
implicit none
integer(qmckl_context), intent(in)  :: context
integer*8             , intent(in)  :: walk_num
integer*8             , intent(in)  :: elec_num
integer*8             , intent(in)  :: nucl_num
integer*8             , intent(in)  :: type_nucl_num
integer*8             , intent(in)  :: type_nucl_vector(nucl_num)
integer*8             , intent(in)  :: cord_num
double precision      , intent(in)  :: rescale_factor_en(type_nucl_num)
double precision      , intent(in)  :: coord_ee(elec_num,3,walk_num)
double precision      , intent(in)  :: coord_en(nucl_num,3)
double precision      , intent(in)  :: en_distance(nucl_num,elec_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_gl(elec_num,4,nucl_num,0:cord_num,walk_num)
double precision,dimension(:,:,:),allocatable :: elnuc_dist_gl
double precision                    :: x, ria_inv, kappa_l
integer*8                           :: i, a, k, l, nw, ii

allocate(elnuc_dist_gl(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_gl             = 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(a, i, nw)
 do ii = 1, 3
   elnuc_dist_gl(ii, i, a) = (coord_ee(i, ii, nw) - coord_en(a, ii)) * ria_inv
 end do
 elnuc_dist_gl(4, i, a) = 2.0d0 * ria_inv
end do
end do

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

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

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

end function qmckl_compute_jastrow_champ_factor_een_rescaled_n_gl_f

Test

assert(qmckl_electron_provided(context));

double een_rescaled_n_gl[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_champ_een_rescaled_n_gl(context, &(een_rescaled_n_gl[0][0][0][0][0]),size_max);

// value of (0,2,1)
assert(fabs( -0.11234061209936878  - een_rescaled_n_gl[0][1][0][0][2])  < 1.e-12);
assert(fabs( 0.0004440109367151707 - een_rescaled_n_gl[0][1][0][0][3])  < 1.e-12);
assert(fabs( -0.012868642597346566 - een_rescaled_n_gl[0][1][0][0][4])  < 1.e-12);
assert(fabs( 0.08601122289922644   - een_rescaled_n_gl[0][2][1][0][3])  < 1.e-12);
assert(fabs( -0.058681563677207206 - een_rescaled_n_gl[0][2][1][0][4])  < 1.e-12);
assert(fabs( 0.005359281880312882  - een_rescaled_n_gl[0][2][1][0][5])  < 1.e-12);

Temporary arrays for electron-electron-nucleus Jastrow $f_{een}$

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

Compute dim_c_vector

Computes the dimension of the vector of parameters.

$N_{ord}$	Number of parameters
2	2
3	6
4	13
5	23
6	37
7	55
8	78
9	106
10	140

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`cord_num`	`int64_t`	in	Order of polynomials
`dim_c_vector`	`int64_t`	out	Number of parameters per atom type

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

info = QMCKL_SUCCESS

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

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

dim_c_vector = 0

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

end function qmckl_compute_dim_c_vector_f

Get

qmckl_exit_code qmckl_get_jastrow_champ_tmp_c(qmckl_context context, double* const tmp_c);
qmckl_exit_code qmckl_get_jastrow_champ_dtmp_c(qmckl_context context, double* const dtmp_c);

Compute c_vector_full

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

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

info = QMCKL_SUCCESS

if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT
if (nucl_num <= 0)                 info = QMCKL_INVALID_ARG_2
if (dim_c_vector < 0)              info = QMCKL_INVALID_ARG_3
if (type_nucl_num <= 0)            info = QMCKL_INVALID_ARG_4
if (info /= QMCKL_SUCCESS)         return

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

end function qmckl_compute_c_vector_full_doc_f

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

if (context == QMCKL_NULL_CONTEXT) return QMCKL_INVALID_CONTEXT;
if (nucl_num <= 0)                 return QMCKL_INVALID_ARG_2;
if (dim_c_vector < 0)              return QMCKL_INVALID_ARG_3;
if (type_nucl_num <= 0)            return QMCKL_INVALID_ARG_4;
if (type_nucl_vector == NULL)      return QMCKL_INVALID_ARG_5;
if (c_vector == NULL)              return QMCKL_INVALID_ARG_6;
if (c_vector_full == NULL)         return QMCKL_INVALID_ARG_7;

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

return QMCKL_SUCCESS;
}

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

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

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

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

Compute lkpm_combined_index

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

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

info = QMCKL_SUCCESS

if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT
if (cord_num < 0)                  info = QMCKL_INVALID_ARG_2
if (dim_c_vector < 0)              info = QMCKL_INVALID_ARG_3
if (info /= QMCKL_SUCCESS)         return

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

end function qmckl_compute_lkpm_combined_index_f

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

int kk, lmax, m;

if (context == QMCKL_NULL_CONTEXT) return QMCKL_INVALID_CONTEXT;
if (cord_num < 0) return QMCKL_INVALID_ARG_2;
if (dim_c_vector < 0) return QMCKL_INVALID_ARG_3;

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

return QMCKL_SUCCESS;
}

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 nuclei
`walk_num`	`int64_t`	in	Number of walkers
`een_rescaled_e`	`double[walk_num][0:cord_num][elec_num][elec_num]`	in	Electron-electron rescaled factor
`een_rescaled_n`	`double[walk_num][0:cord_num][nucl_num][elec_num]`	in	Electron-nucleus rescaled factor
`tmp_c`	`double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]`	out	vector of non-zero coefficients

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

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

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

info = QMCKL_SUCCESS

if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT
if (cord_num <  0)                 info = QMCKL_INVALID_ARG_2
if (elec_num <= 0)                 info = QMCKL_INVALID_ARG_3
if (nucl_num <= 0)                 info = QMCKL_INVALID_ARG_4
if (walk_num <= 0)                 info = QMCKL_INVALID_ARG_5
if (info /= QMCKL_SUCCESS)         return


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

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

end function qmckl_compute_tmp_c_doc_f

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

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 nuclei
`walk_num`	`int64_t`	in	Number of walkers
`een_rescaled_e_gl`	`double[walk_num][0:cord_num][elec_num][4][elec_num]`	in	Electron-electron rescaled factor derivatives
`een_rescaled_n`	`double[walk_num][0:cord_num][nucl_num][elec_num]`	in	Electron-nucleus rescaled factor
`dtmp_c`	`double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]`	out	vector of non-zero coefficients

qmckl_exit_code
qmckl_compute_dtmp_c (const qmckl_context context,
                const int64_t cord_num,
                const int64_t elec_num,
                const int64_t nucl_num,
                const int64_t walk_num,
                const double* een_rescaled_e_gl,
                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_gl,
                             een_rescaled_n, dtmp_c );
#else
return qmckl_compute_dtmp_c_doc (context, cord_num, elec_num, nucl_num, walk_num, een_rescaled_e_gl,
                             een_rescaled_n, dtmp_c );
#endif
}

integer function qmckl_compute_dtmp_c_doc_f( &
context, cord_num, elec_num, nucl_num, &
walk_num, een_rescaled_e_gl, 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_gl(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

info = QMCKL_SUCCESS

if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT
if (cord_num <  0)                 info = QMCKL_INVALID_ARG_2
if (elec_num <= 0)                 info = QMCKL_INVALID_ARG_3
if (nucl_num <= 0)                 info = QMCKL_INVALID_ARG_4
if (walk_num <= 0)                 info = QMCKL_INVALID_ARG_5
if (info /= QMCKL_SUCCESS)         return

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

M = 4*elec_num
N = nucl_num*(cord_num + 1)
K = elec_num
LDA = 4*size(een_rescaled_e_gl,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_gl(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

Test

assert(qmckl_electron_provided(context));

double tmp_c[walk_num][cord_num][cord_num+1][nucl_num][elec_num];
rc = qmckl_get_jastrow_champ_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_champ_dtmp_c(context, &(dtmp_c[0][0][0][0][0][0]));

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

printf("%e\n%e\n", dtmp_c[0][1][0][0][0][0],3.278657e-01);
assert(fabs(dtmp_c[0][1][0][0][0][0] - 3.278657e-01 ) < 1e-6);

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_champ_factor_een(qmckl_context context,
                        double* const factor_een,
                        const int64_t size_max);

Fortran interface

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

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 nuclei
`cord_num`	`int64_t`	in	order of polynomials
`dim_c_vector`	`int64_t`	in	dimension of full coefficient vector
`c_vector_full`	`double[dim_c_vector][nucl_num]`	in	full coefficient vector
`lkpm_combined_index`	`int64_t[4][dim_c_vector]`	in	combined indices
`een_rescaled_e`	`double[walk_num][elec_num][elec_num][0:cord_num]`	in	Electron-nucleus rescaled
`een_rescaled_n`	`double[walk_num][elec_num][nucl_num][0:cord_num]`	in	Electron-nucleus rescaled factor
`factor_een`	`double[walk_num]`	out	Electron-nucleus jastrow

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

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

info = QMCKL_SUCCESS

if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT
if (walk_num <= 0)                 info = QMCKL_INVALID_ARG_2
if (elec_num <= 0)                 info = QMCKL_INVALID_ARG_3
if (nucl_num <= 0)                 info = QMCKL_INVALID_ARG_4
if (cord_num <  0)                 info = QMCKL_INVALID_ARG_5
if (info /= QMCKL_SUCCESS)         return

factor_een = 0.0d0

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

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

end function qmckl_compute_jastrow_champ_factor_een_naive_f

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 nuclei
`cord_num`	`int64_t`	in	order of polynomials
`dim_c_vector`	`int64_t`	in	dimension of full coefficient vector
`c_vector_full`	`double[dim_c_vector][nucl_num]`	in	full coefficient vector
`lkpm_combined_index`	`int64_t[4][dim_c_vector]`	in	combined indices
`tmp_c`	`double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]`	in	vector of non-zero coefficients
`een_rescaled_n`	`double[walk_num][0:cord_num][nucl_num][elec_num]`	in	Electron-nucleus rescaled distances
`factor_een`	`double[walk_num]`	out	Electron-nucleus jastrow

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

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

info = QMCKL_SUCCESS

if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT
if (walk_num <= 0)                 info = QMCKL_INVALID_ARG_2
if (elec_num <= 0)                 info = QMCKL_INVALID_ARG_3
if (nucl_num <= 0)                 info = QMCKL_INVALID_ARG_4
if (cord_num <  0)                 info = QMCKL_INVALID_ARG_5
if (info /= QMCKL_SUCCESS)         return

factor_een = 0.0d0

if (cord_num == 0) return

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

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

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

end function qmckl_compute_jastrow_champ_factor_een_doc_f

Test

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

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

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

Electron-electron-nucleus Jastrow $f_{een}$ derivative

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

TODO: write equations.

Get

qmckl_exit_code
qmckl_get_jastrow_champ_factor_een_gl(qmckl_context context,
                                double* const factor_een_gl,
                                const int64_t size_max);

Fortran interface

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

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 nuclei
`cord_num`	`int64_t`	in	order of polynomials
`dim_c_vector`	`int64_t`	in	dimension of full coefficient vector
`c_vector_full`	`double[dim_c_vector][nucl_num]`	in	full coefficient vector
`lkpm_combined_index`	`int64_t[4][dim_c_vector]`	in	combined indices
`een_rescaled_e`	`double[walk_num][elec_num][elec_num][0:cord_num]`	in	Electron-nucleus rescaled
`een_rescaled_n`	`double[walk_num][elec_num][nucl_num][0:cord_num]`	in	Electron-nucleus rescaled factor
`een_rescaled_e_gl`	`double[walk_num][elec_num][4][elec_num][0:cord_num]`	in	Electron-nucleus rescaled
`een_rescaled_n_gl`	`double[walk_num][elec_num][4][nucl_num][0:cord_num]`	in	Electron-nucleus rescaled factor
`factor_een_gl`	`double[walk_num][4][elec_num]`	out	Electron-nucleus jastrow

integer function qmckl_compute_jastrow_champ_factor_een_gl_naive_f( &
context, walk_num, elec_num, nucl_num, cord_num, dim_c_vector, &
c_vector_full, lkpm_combined_index, een_rescaled_e, een_rescaled_n, &
een_rescaled_e_gl, een_rescaled_n_gl, factor_een_gl)&
result(info)
use qmckl
implicit none
integer(qmckl_context), intent(in)  :: context
integer*8             , intent(in)  :: walk_num, elec_num, cord_num, nucl_num, dim_c_vector
integer*8             , intent(in)  :: lkpm_combined_index(dim_c_vector, 4)
double precision      , intent(in)  :: c_vector_full(nucl_num, dim_c_vector)
double precision      , intent(in)  :: een_rescaled_e(0:cord_num, elec_num, elec_num, walk_num)
double precision      , intent(in)  :: een_rescaled_n(0:cord_num, nucl_num, elec_num, walk_num)
double precision      , intent(in)  :: een_rescaled_e_gl(0:cord_num, elec_num, 4, elec_num, walk_num)
double precision      , intent(in)  :: een_rescaled_n_gl(0:cord_num, nucl_num, 4, elec_num, walk_num)
double precision      , intent(out) :: factor_een_gl(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) info = QMCKL_INVALID_CONTEXT
if (walk_num <= 0)                 info = QMCKL_INVALID_ARG_2
if (elec_num <= 0)                 info = QMCKL_INVALID_ARG_3
if (nucl_num <= 0)                 info = QMCKL_INVALID_ARG_4
if (cord_num <  0)                 info = QMCKL_INVALID_ARG_5
if (info /= QMCKL_SUCCESS)         return

factor_een_gl = 0.0d0

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

do a = 1, nucl_num
 cn = c_vector_full(a, n)
 do j = 1, elec_num
   accu = 0.0d0
   accu2 = 0.0d0
   daccu = 0.0d0
   daccu2 = 0.0d0
   do i = 1, elec_num
      accu = accu + een_rescaled_e(k, i, j, nw) * een_rescaled_n(m, a, i, nw)
      accu2 = accu2 + een_rescaled_e(k, i, j, nw) * een_rescaled_n(m + l, a, i, nw)
      daccu(1:4) = daccu(1:4) + een_rescaled_e_gl(k, j, 1:4, i, nw) *   &
           een_rescaled_n(m, a, i, nw)
      daccu2(1:4) = daccu2(1:4) + een_rescaled_e_gl(k, j, 1:4, i, nw) * &
           een_rescaled_n(m + l, a, i, nw)
   end do
   factor_een_gl(j, 1:4, nw) = factor_een_gl(j, 1:4, nw) +   &
        (accu * een_rescaled_n_gl(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_gl(m, a, 1:4, j, nw)) * cn
   
   factor_een_gl(j, 4, nw) = factor_een_gl(j, 4, nw) + 2.0d0 * ( &
        daccu (1) * een_rescaled_n_gl(m + l, a, 1, j, nw) +      &
        daccu (2) * een_rescaled_n_gl(m + l, a, 2, j, nw) +      &
        daccu (3) * een_rescaled_n_gl(m + l, a, 3, j, nw) +      &
        daccu2(1) * een_rescaled_n_gl(m, a, 1, j, nw    ) +      &
        daccu2(2) * een_rescaled_n_gl(m, a, 2, j, nw    ) +      &
        daccu2(3) * een_rescaled_n_gl(m, a, 3, j, nw    ) ) * cn

 end do
end do
end do
end do

end function qmckl_compute_jastrow_champ_factor_een_gl_naive_f

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 nuclei
`cord_num`	`int64_t`	in	order of polynomials
`dim_c_vector`	`int64_t`	in	dimension of full coefficient vector
`c_vector_full`	`double[dim_c_vector][nucl_num]`	in	full coefficient vector
`lkpm_combined_index`	`int64_t[4][dim_c_vector]`	in	combined indices
`tmp_c`	`double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]`	in	Temporary intermediate tensor
`dtmp_c`	`double[walk_num][0:cord_num-1][0:cord_num][nucl_num][4][elec_num]`	in	vector of non-zero coefficients
`een_rescaled_n`	`double[walk_num][0:cord_num][nucl_num][elec_num]`	in	Electron-nucleus rescaled factor
`een_rescaled_n_gl`	`double[walk_num][0:cord_num][nucl_num][4][elec_num]`	in	Derivative of Electron-nucleus rescaled factor
`factor_een_gl`	`double[walk_num][4][elec_num]`	out	Derivative of Electron-nucleus jastrow

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

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

info = QMCKL_SUCCESS

if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT
if (walk_num <= 0)                 info = QMCKL_INVALID_ARG_2
if (elec_num <= 0)                 info = QMCKL_INVALID_ARG_3
if (nucl_num <= 0)                 info = QMCKL_INVALID_ARG_4
if (cord_num <  0)                 info = QMCKL_INVALID_ARG_5
if (info /= QMCKL_SUCCESS)         return

factor_een_gl = 0.0d0

if (cord_num == 0) return

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

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

      do ii = 1, 4
         do j = 1, elec_num
            factor_een_gl(j,ii,nw) = factor_een_gl(j,ii,nw) + (          &
                 tmp_c(j,a,m,k,nw)       * een_rescaled_n_gl(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_gl(j,ii,a,m,nw)     &
                 ) * cn
         end do
      end do

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

end function qmckl_compute_jastrow_champ_factor_een_gl_doc_f

qmckl_exit_code qmckl_compute_jastrow_champ_factor_een_gl_doc (
      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_c_vector,
      const double* c_vector_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_gl,
      double* const factor_een_gl );

qmckl_exit_code qmckl_compute_jastrow_champ_factor_een_gl (
      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_c_vector,
      const double* c_vector_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_gl,
      double* const factor_een_gl );

qmckl_exit_code
qmckl_compute_jastrow_champ_factor_een_gl(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_c_vector,
                                            const double *c_vector_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_gl,
                                            double* const factor_een_gl)
{
#ifdef HAVE_HPC
return qmckl_compute_jastrow_champ_factor_een_gl_hpc(context, walk_num, elec_num, nucl_num,
                                                         cord_num, dim_c_vector, c_vector_full,
                                                         lkpm_combined_index, tmp_c, dtmp_c,
                                                         een_rescaled_n, een_rescaled_n_gl,
                                                         factor_een_gl);
#else
return qmckl_compute_jastrow_champ_factor_een_gl_doc(context, walk_num, elec_num, nucl_num,
                                                         cord_num, dim_c_vector, c_vector_full,
                                                         lkpm_combined_index, tmp_c, dtmp_c,
                                                         een_rescaled_n, een_rescaled_n_gl,
                                                         factor_een_gl);
#endif
}

Test

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

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

printf("%20.15e\n", factor_een_gl[0][0][0]);
assert(fabs(8.967809309100624e-02 - factor_een_gl[0][0][0]) < 1e-12);

printf("%20.15e\n", factor_een_gl[1][0][1]);
assert(fabs(-3.401545636077585e-02 - factor_een_gl[1][0][1]) < 1e-12);

printf("%20.15e\n", factor_een_gl[2][0][2]);
assert(fabs(-2.631321052321952e-01 - factor_een_gl[2][0][2]) < 1e-12);

printf("%20.15e\n", factor_een_gl[3][0][3]);
assert(fabs(-1.016785559040419e+00 - factor_een_gl[3][0][3]) < 1e-12);

Total Jastrow

Value

Value of the total Jastrow factor: $\exp(J)$

Get

qmckl_exit_code
qmckl_get_jastrow_champ_value(qmckl_context context,
                       double* const value,
                       const int64_t size_max);

Fortran interface

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

Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`walk_num`	`int64_t`	in	Number of walkers
`f_ee`	`double[walk_num]`	in	ee component
`f_en`	`double[walk_num]`	in	eN component
`f_een`	`double[walk_num]`	in	eeN component
`value`	`double[walk_num]`	out	Total Jastrow factor

integer function qmckl_compute_jastrow_champ_value_doc_f(context, &
walk_num, f_ee, f_en, f_een, value) &
result(info)
use qmckl
implicit none
integer(qmckl_context), intent(in)  :: context
integer*8             , intent(in)  :: walk_num
double precision      , intent(in)  :: f_ee(walk_num), f_en(walk_num), f_een(walk_num)
double precision      , intent(out) :: value(walk_num)

integer*8 :: i

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

do i = 1, walk_num
value(i) = f_ee(i) + f_en(i) + f_een(i)
end do

do i = 1, walk_num
! Flush to zero to avoid floating-point exception
if (value(i) < -100.d0) then
  value(i) = 0.d0
else
  value(i) = dexp(value(i))
endif
end do

end function qmckl_compute_jastrow_champ_value_doc_f

qmckl_exit_code qmckl_compute_jastrow_champ_value (
      const qmckl_context context,
      const int64_t walk_num,
      const double* f_ee,
      const double* f_en,
      const double* f_een,
      double* const value );

qmckl_exit_code qmckl_compute_jastrow_champ_value_doc (
      const qmckl_context context,
      const int64_t walk_num,
      const double* f_ee,
      const double* f_en,
      const double* f_een,
      double* const value );

qmckl_exit_code qmckl_compute_jastrow_champ_value_hpc (
      const qmckl_context context,
      const int64_t walk_num,
      const double* f_ee,
      const double* f_en,
      const double* f_een,
      double* const value );

qmckl_exit_code qmckl_compute_jastrow_champ_value (
         const qmckl_context context,
         const int64_t walk_num,
         const double* factor_ee,
         const double* factor_en,
         const double* factor_een,
         double* const value)
{

#ifdef HAVE_HPC
 return qmckl_compute_jastrow_champ_value_hpc (
         context, walk_num, factor_ee, factor_en, factor_een, value);
#else
 return qmckl_compute_jastrow_champ_value_doc (
         context, walk_num, factor_ee, factor_en, factor_een, value);
#endif
}

Test

printf("Total Jastrow value\n");
/* Check if Jastrow is properly initialized */
assert(qmckl_jastrow_champ_provided(context));

rc = qmckl_check(context,
           qmckl_get_jastrow_champ_factor_ee(context, &(factor_ee[0]), walk_num)
           );
assert(rc == QMCKL_SUCCESS);

rc = qmckl_check(context,
           qmckl_get_jastrow_champ_factor_en(context, &(factor_en[0]), walk_num)
           );
assert(rc == QMCKL_SUCCESS);

rc = qmckl_check(context,
           qmckl_get_jastrow_champ_factor_een(context, &(factor_een[0]), walk_num)
           );
assert(rc == QMCKL_SUCCESS);

double total_j[walk_num];
rc = qmckl_check(context,
           qmckl_get_jastrow_champ_value(context, &(total_j[0]), walk_num)
           );
assert(rc == QMCKL_SUCCESS);


for (int64_t i=0 ; i< walk_num ; ++i) {
assert (total_j[i] - exp(factor_ee[i] + factor_en[i] + factor_een[i]) < 1.e-12);
}

Derivatives

Gradients and Laplacian of the total Jastrow factor: \[ \nabla \left[ e^{J(\mathbf{r})} \right] = e^{J(\mathbf{r})} \nabla J(\mathbf{r}) \] \[ \Delta \left[ e^{J(\mathbf{r})} \right] = e^{J(\mathbf{r})} \left[ \Delta J(\mathbf{r}) + \nabla J(\mathbf{r}) \cdot \nabla J(\mathbf{r}) \right] \]

Get

qmckl_exit_code
qmckl_get_jastrow_champ_gl(qmckl_context context,
                       double* const gl,
                       const int64_t size_max);

Fortran interface

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

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
`value`	`double[walk_num]`	in	Total Jastrow
`gl_ee`	`double[walk_num][4][elec_num]`	in	ee component
`gl_en`	`double[walk_num][4][elec_num]`	in	eN component
`gl_een`	`double[walk_num][4][elec_num]`	in	eeN component
`gl`	`double[walk_num][4][elec_num]`	out	Total Jastrow factor

integer function qmckl_compute_jastrow_champ_gl_doc_f(context, &
    walk_num, elec_num, value, gl_ee, gl_en, gl_een, gl) &
    result(info)
 use qmckl
 implicit none
 integer(qmckl_context), intent(in)  :: context
 integer*8             , intent(in)  :: walk_num, elec_num
 double precision      , intent(in)  :: value (walk_num)
 double precision      , intent(in)  :: gl_ee (elec_num,4,walk_num)
 double precision      , intent(in)  :: gl_en (elec_num,4,walk_num)
 double precision      , intent(in)  :: gl_een(elec_num,4,walk_num)
 double precision      , intent(out) :: gl    (elec_num,4,walk_num)

 integer*8 :: i, j, k

 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

 do k = 1, walk_num
    do j=1,4
       do i = 1, elec_num
          gl(i,j,k) = gl_ee(i,j,k) + gl_en(i,j,k) + gl_een(i,j,k)
       end do
    end do
    do i = 1, elec_num
        gl(i,4,k) = gl(i,4,k) + &
                    gl(i,1,k) * gl(i,1,k) + &
                    gl(i,2,k) * gl(i,2,k) + &
                    gl(i,3,k) * gl(i,3,k)
    end do
    gl(:,:,k) = gl(:,:,k) * value(k)
 end do


end function qmckl_compute_jastrow_champ_gl_doc_f

qmckl_exit_code qmckl_compute_jastrow_champ_gl (
      const qmckl_context context,
      const int64_t walk_num,
      const int64_t elec_num,
      const double* value,
      const double* gl_ee,
      const double* gl_en,
      const double* gl_een,
      double* const gl );

qmckl_exit_code qmckl_compute_jastrow_champ_gl_doc (
      const qmckl_context context,
      const int64_t walk_num,
      const int64_t elec_num,
      const double* value,
      const double* gl_ee,
      const double* gl_en,
      const double* gl_een,
      double* const gl );

qmckl_exit_code qmckl_compute_jastrow_champ_gl_hpc (
      const qmckl_context context,
      const int64_t walk_num,
      const int64_t elec_num,
      const double* value,
      const double* gl_ee,
      const double* gl_en,
      const double* gl_een,
      double* const gl );

qmckl_exit_code qmckl_compute_jastrow_champ_gl (
     const qmckl_context context,
     const int64_t walk_num,
     const int64_t elec_num,
     const double* value,
     const double* gl_ee,
     const double* gl_en,
     const double* gl_een,
     double* const gl)
{

#ifdef HAVE_HPC
 return qmckl_compute_jastrow_champ_gl_hpc (context,
     walk_num, elec_num, value, gl_ee, gl_en, gl_een, gl);
#else
 return qmckl_compute_jastrow_champ_gl_doc (context,
     walk_num, elec_num, value, gl_ee, gl_en, gl_een, gl);
#endif
}

Test

printf("Total Jastrow derivatives\n");
/* Check if Jastrow is properly initialized */
assert(qmckl_jastrow_champ_provided(context));

rc = qmckl_check(context,
           qmckl_get_jastrow_champ_factor_ee_gl(context, &(factor_ee_gl[0][0][0]), walk_num*elec_num*4)
           );
assert(rc == QMCKL_SUCCESS);

rc = qmckl_check(context,
           qmckl_get_jastrow_champ_factor_en_gl(context, &(factor_en_gl[0][0][0]), walk_num*elec_num*4)
           );
assert(rc == QMCKL_SUCCESS);

rc = qmckl_check(context,
           qmckl_get_jastrow_champ_factor_een_gl(context, &(factor_een_gl[0][0][0]), walk_num*elec_num*4)
           );
assert(rc == QMCKL_SUCCESS);

double total_j_deriv[walk_num][4][elec_num];
rc = qmckl_check(context,
           qmckl_get_jastrow_champ_gl(context, &(total_j_deriv[0][0][0]), walk_num*elec_num*4)
           );
assert(rc == QMCKL_SUCCESS);

rc = qmckl_check(context,
           qmckl_get_jastrow_champ_value(context, &(total_j[0]), walk_num)
           );
assert(rc == QMCKL_SUCCESS);


for (int64_t k=0 ; k< walk_num ; ++k) {
for (int64_t m=0 ; m<4; ++m) {
for (int64_t e=0 ; e<elec_num; ++e) {
if (m < 3) { /* test only gradients */
  assert (total_j_deriv[k][m][e]/total_j[k] - (factor_ee_gl[k][m][e] + factor_en_gl[k][m][e] + factor_een_gl[k][m][e]) < 1.e-12);
}
}
}
}

393 KiB Raw Blame History

CHAMP Jastrow Factor

Introduction

Context

Data structure

Initialization functions

Fortran interface

Access functions

Fortran interface

Test

Computation

Electron-electron component

Asymptotic component

Get

Fortran interface

Compute

Test

Electron-electron component

Get

Fortran interface

Compute

Test

Derivative

Get

Fortran interface

Compute

Test

Electron-electron rescaled distances

Get

Compute

Test

Electron-electron rescaled distance gradients and Laplacian with respect to electron coordinates

Get

Compute

Test

Electron-nucleus component

Asymptotic component for

Get

Fortran interface

Compute

Test

Electron-nucleus component

Get

Fortran interface

Compute

Test

Derivative

Get

Fortran interface

Compute

Test

Electron-nucleus rescaled distances

Get

Compute

Test

Electron-electron rescaled distance gradients and Laplacian with respect to electron coordinates

Get

Compute

Test

Electron-electron-nucleus component

Electron-electron rescaled distances in $J_\text{eeN}$

Get

Compute

Test

Electron-electron rescaled distances derivatives in $J_\text{eeN}$

Get

Compute

Test

Electron-nucleus rescaled distances in $J_\text{eeN}$

Get

Compute

Test

Electron-nucleus rescaled distances derivatives in $J_\text{eeN}$

Get

Compute

Test

Temporary arrays for electron-electron-nucleus Jastrow $f_{een}$

Compute dim_c_vector

Get

Compute c_vector_full

393 KiB

Raw Blame History