CHAMP Jastrow Factor

1. Introduction
- 1.1. Reformulation of the three-body part
2.
3. Context
4. Computation

1. Introduction

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

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

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

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

\[ J_{\text{eN}}(\mathbf{r},\mathbf{R}) = \sum_{\alpha=1}^{N_\text{nucl}} \sum_{i=1}^{N_\text{elec}} \frac{a_{1\,\alpha}\, 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{e}({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. Warning: The types of nuclei use zero-based indexing.

The derivatives are computed with respect to the electron \(i\) for \[ r_{ij} = |r_i - r_j| \]

1.1. Reformulation of the three-body part

To accelerate the computation of the three-body part, the Jastrow factor is re-expressed as follows, with \(m=(p-k)/2 -l/2\):

\begin{eqnarray*} J_{kpl} & = & \sum_{\alpha=1}^{N_{\text{nucl}}} \sum_{i=1}^{N_{\text{elec}}} \sum_{j=1}^{i-1} c_{lkp\alpha} \left[ g_\text{e}({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]^{m} \\ & = & \frac{1}{2} \sum_{\alpha=1}^{N_{\text{nucl}}} \sum_{i=1}^{N_{\text{elec}}} \sum_{j=1}^{N_{\text{elec}}} c_{lkp\alpha} \left[ g_\text{e}({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]^{m} \\ & = & \frac{1}{2} \sum_{\alpha=1}^{N_{\text{nucl}}} \sum_{i=1}^{N_{\text{elec}}} \sum_{j=1}^{N_{\text{elec}}} c_{lkp\alpha} \left[ g_\text{e}({r}_{ij}) \right]^k \left[ \left[ g_\alpha({R}_{i\alpha}) \right]^{l+m} \left[ g_\alpha({R}_{j\alpha}) \right]^{m} + \left[ g_\alpha({R}_{i\alpha}) \right]^{l} \left[ g_\alpha({R}_{j\alpha}) \right]^{l+m} \right] \\ & = & \sum_{\alpha=1}^{N_{\text{nucl}}} c_{lkp\alpha} \sum_{i=1}^{N_{\text{elec}}} \sum_{j=1}^{N_{\text{elec}}} \left[ g_\alpha({R}_{i\alpha}) \right]^{l+m} \left[ g_\text{e}({r}_{ij}) \right]^k \left[ g_\alpha({R}_{j\alpha}) \right]^{m} \\ & = & \sum_{\alpha=1}^{N_{\text{nucl}}} c_{lkp\alpha} \sum_{i=1}^{N_{\text{elec}}} \left[ g_\alpha({R}_{i\alpha}) \right]^{l+m} P_{i\alpha}^{km}, \text{ with } P_{i\alpha}^{km} = \sum_{j=1}^{N_{\text{elec}}} \left[ g_\text{e}({r}_{ij}) \right]^k \left[ g_\alpha({R}_{j\alpha}) \right]^{m}. \\ J & = & \sum_{p=2}^{N_{\text{ord}}} \sum_{k=0}^{p-1} \sum_{l=0}^{p-k-2\delta_{k,0}} \sum_{\alpha=1}^{N_{\text{nucl}}} c_{lkp\alpha} \sum_{i=1}^{N_{\text{elec}}} \left[ g_\alpha({R}_{i\alpha}) \right]^{(p-k+l)/2} P_{i\alpha}^{k, (p-k-l)/2} \end{eqnarray*}

The computation of \(P\) scales as \(\mathcal{O}(N_\text{elec}^2 N_\text{nucl}n^2)\), and the computation of \(J\) scales as \(\mathcal{O}(N_\text{elec}N_\text{nucl}n^2)\).

2. 3. 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
`spin_independent`	`int32_t`	If 1, use same parameters for parallel and anti-parallel spins. Otherwise, 0.

Computed data:

Variable	Type	In/Out
`dim_c_vector`	`int64_t`	Number of unique C coefficients
`dim_c_vector_date`	`uint64_t`	Last modification of the number of unique C coefficients
`asymp_jasa`	`double[type_nucl_num]`	Asymptotic component
`asymp_jasa_date`	`uint64_t`	Last modification of the asymptotic component
`asymp_jasa_pderiv`	`double[type_nucl_num][aord_num+1]`	Derivatives w.r.t. the A parameters of the Asymptotic component
`asymp_jasa_pderiv_date`	`uint64_t`	Last modification of the asymptotic component derivative
`asymp_jasb`	`double[2]`	Asymptotic component
`asymp_jasb_date`	`uint64_t`	Last modification of derivative of the asymptotic jastrow
`asymp_jasb_pderiv`	`double[2][bord_num+1]`	Derivative of the asymptotic jastrow w.r.t the jastrow parameters (up- or down-spin)
`asymp_jasb_pderiv_date`	`double[2][bord_num+1]`	Last modification of psrameter derivative of the asymptotic jastrow
`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][num][num][4]`	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][nucl_num][num][4]`	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_ee_pderiv`	`double[bord_num+1]`	Parameter derivatives Jastrow factor: electron-electron part
`factor_ee_pderiv_date`	`uint64_t`	Parameter derivatives 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_en_pderiv`	`double[type_nucl_num][aord_num+1]`	Parameter derivatives Jastrow factor: electron-nucleus part
`factor_en_pderiv_date`	`uint64_t`	Parameter derivativesJastrow 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_een_pderiv`	`double[type_nucl_num][dim_c_vector]`	Parameter derivatives Jastrow factor: electron-electron-nucleus part
`factor_een_pderiv_date`	`uint64_t`	Parameter derivatives Jastrow factor: electron-electron-nucleus part
`factor_ee_gl`	`double[walk_num][4][elec_num]`	Parameter derivatives Jastrow factor: electron-electron-nucleus part
`factor_ee_gl_date`	`uint64_t`	Keep track of the date for the derivative
`factor_ee_gl_pderiv`	`double[bord_num+1][elec_num][4]`	Parameter derivative of the gradient of the jastrow factor	Derivative of the gradient and Laplacian with respect to the Jastrow parameters
`factor_ee_gl_pderiv_date`	`uint64_t`	Parameter derivative of the gradient of the jastrow factor	Derivative of the gradient and Laplacian with respect to the Jastrow parameters
`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_en_gl_pderiv`	`double[type_nucl_num][aord_num+1][elec_num][4]`	Derivative of the Jastrow factor: electron-electron-nucleus part
`factor_en_gl_pderiv_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
`factor_een_grad`	`double[walk_num][3][elec_num]`	Gradient of the Jastrow factor: electron-electron-nucleus part
`factor_een_grad_date`	`uint64_t`	Keep track of the date for the een derivative
`factor_een_gl_pderiv`	`double[type_nucl_num][dim_c_vector][elec_num][4]`	Parameter Derivative of the Jastrow factor: electron-electron-nucleus part
`factor_een_gl_pderiv_date`	`uint64_t`	Keep track of the date for the en 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
`grad`	`double[walk_num][3][elec_num]`	Gradient of the Jastrow factor
`value_date`	`uint64_t`	Keep track of the date

3.1. Data structure

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.

3.2. Initialization functions

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

qmckl_exit_code  qmckl_set_jastrow_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 size_max);
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);
qmckl_exit_code  qmckl_set_jastrow_champ_spin_independent  (qmckl_context context, const int32_t spin_independent);

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.

3.2.0.1. 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
     real(c_double),   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
     real(c_double),   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
     real(c_double),   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
     real(c_double),   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
     real(c_double),   intent(in) :: c_vector(size_max)
   end function qmckl_set_jastrow_champ_c_vector

   integer(qmckl_exit_code) function qmckl_set_jastrow_champ_spin_independent(context, &
        spin_independent) bind(C)
     use, intrinsic :: iso_c_binding
     import
     implicit none
     integer(qmckl_context) , intent(in)  , value :: context
     integer(c_int32_t),    intent(in), value  :: spin_independent
   end function qmckl_set_jastrow_champ_spin_independent

end interface

3.3. Access functions

Along with these core functions, calculation of the jastrow factor requires the following additional information to be set:

When all the data for the AOs have been provided, the following function returns true.

bool      qmckl_jastrow_champ_provided           (const qmckl_context context);

3.3.0.1. 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
     real(c_double),   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
     real(c_double),   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
     real(c_double),   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
     real(c_double),   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
     real(c_double),   intent(out)                 :: c_vector(size_max)
   end function qmckl_get_jastrow_champ_c_vector

   integer(qmckl_exit_code) function qmckl_get_jastrow_champ_spin_independent(context, &
        spin_independent) bind(C)
     use, intrinsic :: iso_c_binding
     import
     implicit none
     integer(qmckl_context) , intent(in)  , value :: context
     integer(c_int32_t),   intent(out)                :: spin_independent
   end function qmckl_get_jastrow_champ_spin_independent

end interface

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

4.1. Electron-electron component

4.1.1. Asymptotic component

Calculate the asymptotic component asymp_jasb to be subtracted 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. If the spin_independent variable is set to 1, then \(\delta^{\uparrow \downarrow}\) is always equal to one.

\[ 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} \]

4.1.1.1. 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
     real(c_double),   intent(out)               :: asymp_jasb(size_max)
   end function qmckl_get_jastrow_champ_asymp_jasb
end interface

4.1.1.2. Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`bord_num`	`int64_t`	in	Order of the polynomial
`b_vector`	`double[bord_num+1]`	in	Values of b
`rescale_factor_ee`	`double`	in	Electron coordinates
`spin_independent`	`int32_t`	in	If 1, same parameters for parallel and anti-parallel pairs
`asymp_jasb`	`double[2]`	out	Asymptotic value

function qmckl_compute_jastrow_champ_asymp_jasb_doc(context, &
     bord_num, b_vector, rescale_factor_ee, spin_independent, asymp_jasb) &
     bind(C) result(info)
  use qmckl
  implicit none

  integer (qmckl_context) , intent(in)  , value :: context
  integer (c_int64_t) , intent(in)  , value :: bord_num
  real    (c_double ) , intent(in)          :: b_vector(bord_num+1)
  real    (c_double ) , intent(in)  , value :: rescale_factor_ee
  integer (c_int32_t) , intent(in)  , value :: spin_independent
  real    (c_double ) , intent(out)         :: asymp_jasb(2)
  integer(qmckl_exit_code)                  :: info

  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)
  if (spin_independent == 1) then
     asymp_jasb(:) = (/asym_one, asym_one/)
  else
     asymp_jasb(:) = (/0.5d0*asym_one, asym_one/)
  end if

  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

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

4.1.2.1. Get

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

4.1.2.2. Compute

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

function qmckl_compute_ee_distance_rescaled_doc(context, &
     elec_num, rescale_factor_ee, walk_num, &
     coord, ee_distance_rescaled) &
     bind(C) result(info)
  use qmckl
  implicit none

  integer(qmckl_context), intent(in), value :: context
  integer (c_int64_t) , intent(in)  , value :: elec_num
  real    (c_double ) , intent(in)  , value :: rescale_factor_ee
  integer (c_int64_t) , intent(in)  , value :: walk_num
  real    (c_double ) , intent(in)          :: coord(elec_num,walk_num,3)
  real    (c_double ) , intent(out)         :: ee_distance_rescaled(elec_num,elec_num,walk_num)
  integer(qmckl_exit_code)                  :: info

  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_doc

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

4.1.3.1. Get

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

4.1.3.2. Compute

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

function qmckl_compute_ee_distance_rescaled_gl_doc(context,  &
     elec_num, rescale_factor_ee, walk_num, coord, ee_distance_rescaled_gl) &
     bind(C) result(info)
  use qmckl
  implicit none

  integer(qmckl_context), intent(in), value :: context
  integer (c_int64_t) , intent(in)  , value :: elec_num
  real    (c_double ) , intent(in)  , value :: rescale_factor_ee
  integer (c_int64_t) , intent(in)  , value :: walk_num
  real    (c_double ) , intent(in)          :: coord(elec_num,walk_num,3)
  real    (c_double ) , intent(out)         :: ee_distance_rescaled_gl(4,elec_num,elec_num,walk_num)
  integer(qmckl_exit_code)                  :: info

  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_doc

4.1.4. 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. If the spin_independent variable is set to 1, then \(\delta^{\uparrow \downarrow}\) is always equal to one.

\[ f_\text{ee} = \sum_{i,j

\(\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] \]

4.1.4.1. 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
     real(c_double),   intent(out)               :: factor_ee(size_max)
   end function qmckl_get_jastrow_champ_factor_ee
end interface

4.1.4.2. Compute

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

function qmckl_compute_jastrow_champ_factor_ee_doc(context, &
     walk_num, elec_num, up_num, bord_num, b_vector, &
     ee_distance_rescaled, asymp_jasb, spin_independent, factor_ee) &
     bind(C) result(info)
  use qmckl
  implicit none

  integer (qmckl_context), intent(in), value :: context
  integer (c_int64_t)    , intent(in), value :: walk_num
  integer (c_int64_t)    , intent(in), value :: elec_num
  integer (c_int64_t)    , intent(in), value :: up_num
  integer (c_int64_t)    , intent(in), value :: bord_num
  real    (c_double )    , intent(in)        :: b_vector(bord_num+1)
  real    (c_double )    , intent(in)        :: ee_distance_rescaled(elec_num,elec_num,walk_num)
  real    (c_double )    , intent(in)        :: asymp_jasb(2)
  integer (c_int32_t)    , intent(in), value :: spin_independent
  real    (c_double )    , intent(out)       :: factor_ee(walk_num)
  integer(qmckl_exit_code)                   :: info

  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 (spin_independent == 1) then
              factor_ee(nw) = factor_ee(nw) + b_vector(1) * x / (1.d0 + b_vector(2) * x) - asymp_jasb(2)
           else
              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
           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

4.1.5. 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] \]

4.1.5.1. 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
     real(c_double),   intent(out)               :: factor_ee_gl(size_max)
   end function qmckl_get_jastrow_champ_factor_ee_gl
end interface

4.1.5.2. Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`walk_num`	`int64_t`	in	Number of walkers
`elec_num`	`int64_t`	in	Number of electrons
`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][elec_num][elec_num][4]`	in	Electron-electron distances
`spin_independent`	`int32_t`	in	If 1, same parameters for parallel and antiparallel spins
`factor_ee_gl`	`double[walk_num][4][elec_num]`	out	Electron-electron distances

function 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,  &
     spin_independent, factor_ee_gl) &
     bind(C) result(info)
  use qmckl
  implicit none

  integer (qmckl_context), intent(in), value :: context
  integer (c_int64_t) , intent(in)  , value :: walk_num
  integer (c_int64_t) , intent(in)  , value :: elec_num
  integer (c_int64_t) , intent(in)  , value :: up_num
  integer (c_int64_t) , intent(in)  , value :: bord_num
  real    (c_double ) , intent(in)          :: b_vector(bord_num+1)
  real    (c_double ) , intent(in)          :: ee_distance_rescaled(elec_num,elec_num,walk_num)
  real    (c_double ) , intent(in)          :: ee_distance_rescaled_gl(4,elec_num,elec_num,walk_num)
  integer (c_int32_t) , intent(in)  , value :: spin_independent
  real    (c_double ) , intent(out)         :: factor_ee_gl(elec_num,4,walk_num)
  integer(qmckl_exit_code)                  :: info

  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

  if ((spin_independent < 0).or.(spin_independent > 1)) then
     info = QMCKL_INVALID_ARG_8
     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 (spin_independent == 1) then
              f = b_vector(1) * invdenom2
           else
              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
           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

4.1.6. Parameter Derivative of the Asymptotic component

Calculate the derivatives of the asymptotic component w.r.t. the jastrow parameters asymp_jasb_pderiv to be subtracted from the electron-electron jastrow factor \(J_{\text{ee}}\) derivatives. Two values are computed. The first one is for parallel spin pairs, and the second one for antiparallel spin pairs. If the spin_independent variable is set to 1, then \(\delta^{\uparrow \downarrow}\) is always equal to one.

\[ J_{\text{ee}}^{\infty} = \frac{\frac{1}{2}(1+\delta^{\uparrow \downarrow})\,B_0 \kappa_\text{ee}^{-1}}{1 + B_1\, \kappa_\text{ee}^{-1}} + \sum_{p=2}^{N_\text{ord}^B} B_{p+1}\, \kappa_\text{ee}^{-p} \]

\[ \partial_{B_0} J_{\text{ee}}^\infty = \frac{\frac{1}{2}(1 + \delta^{\uparrow\downarrow}_{ij}) \kappa_{\text{ee}}^{-1}}{1 + B_1 \kappa_{\text{ee}}^{-1}} \]

\[ \partial_{B_1}J_{\text{ee}}^\infty = -\frac{ \frac{1}{2}(1 + \delta_{ij}) \kappa_{\text{ee}}^{-2}} {(1 + B_1 \kappa_{\text{ee}}^{-1})^2 } \]

\[ \partial_{B_k} J_{\text{ee}}^\infty = \kappa_{\text{ee}}^{-(k-1)} \] for \(k > 1\)

4.1.6.1. Get

qmckl_exit_code
qmckl_get_jastrow_champ_asymp_jasb_pderiv(qmckl_context context,
                                          double* const asymp_jasb_pderiv,
                                          const int64_t size_max);

Fortran interface

interface
   integer(qmckl_exit_code) function qmckl_get_jastrow_champ_asymp_jasb_pderiv(context, &
        asymp_jasb_pderiv, 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
     real(c_double),   intent(out)               :: asymp_jasb_pderiv(size_max)
   end function qmckl_get_jastrow_champ_asymp_jasb_pderiv
end interface

4.1.6.2. Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`bord_num`	`int64_t`	in	Order of the polynomial
`b_vector`	`double[bord_num+1]`	in	Values of b
`rescale_factor_ee`	`double`	in	Electron coordinates
`spin_independent`	`int32_t`	in	If 1, same parameters for parallel and anti-parallel pairs
`asymp_jasb`	`double[2][bord_num+1]`	out	Asymptotic value

function qmckl_compute_jastrow_champ_asymp_jasb_pderiv_doc(context, &
     bord_num, b_vector, rescale_factor_ee, spin_independent, asymp_jasb_pderiv) &
     bind(C) result(info)
  use qmckl
  implicit none

  integer (qmckl_context) , intent(in)  , value :: context
  integer (c_int64_t) , intent(in)  , value :: bord_num
  real    (c_double ) , intent(in)          :: b_vector(bord_num+1)
  real    (c_double ) , intent(in)  , value :: rescale_factor_ee
  integer (c_int32_t) , intent(in)  , value :: spin_independent
  real    (c_double ) , intent(out)         :: asymp_jasb_pderiv(bord_num+1,2)
  integer(qmckl_exit_code)                  :: info

  integer*8 :: i, p
  double precision   :: kappa_inv, x, asym_one, asym_two
  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 = kappa_inv / (1.0d0 + b_vector(2) * kappa_inv)
  asym_two = -b_vector(1) * kappa_inv**2 / (1.0d0 + b_vector(2) * kappa_inv)**2
  if (spin_independent == 1) then
     asymp_jasb_pderiv(1,:) = (/asym_one, asym_one/)
     asymp_jasb_pderiv(2,:) = (/asym_two, asym_two/)
  else
     asymp_jasb_pderiv(1,:) = (/0.5d0*asym_one, asym_one/)
     asymp_jasb_pderiv(2,:) = (/0.5d0*asym_two, asym_two/)
  end if

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

end function qmckl_compute_jastrow_champ_asymp_jasb_pderiv_doc

4.1.7. Parameter Derivative

The derivative of factor_ee with respect to the jastrow parameters b_vector is computed. These derivatives are computed fron the ee_rescaled_distances and the parameter derivatives of the asymptotic part asymp_jasb_pderiv. If spin_independent is set to 1, then \(\delta^{\uparrow \downarrow}\) is always equal to 1.

\[ \partial_{B_0} f_\text{ee} = \sum_{j,i

\[ \partial_{B_1} f_\text{ee} = -\sum_{j,i

\[ \partial_{B_k} f_\text{ee} = \sum_{j, i 1 \]

4.1.7.1. Get

qmckl_exit_code
qmckl_get_jastrow_champ_factor_ee_pderiv(qmckl_context context,
                                         double* const factor_ee_pderiv,
                                         const int64_t size_max);

interface
   integer(qmckl_exit_code) function qmckl_get_jastrow_champ_factor_ee_pderiv (context, &
        factor_ee_pderiv, 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
     real(c_double),   intent(out)               :: factor_ee_pderiv(size_max)
   end function qmckl_get_jastrow_champ_factor_ee_pderiv
end interface

4.1.7.2. Compute

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

function qmckl_compute_jastrow_champ_factor_ee_pderiv_doc(context, &
     walk_num, elec_num, up_num, bord_num, b_vector, &
     ee_distance_rescaled, asymp_jasb_pderiv, spin_independent, factor_ee_pderiv) &
     bind(C) result(info)

  use, intrinsic :: iso_c_binding
  use qmckl
  implicit none

  integer (qmckl_context), intent(in), value :: context
  integer (c_int64_t)    , intent(in), value :: walk_num
  integer (c_int64_t)    , intent(in), value :: elec_num
  integer (c_int64_t)    , intent(in), value :: up_num
  integer (c_int64_t)    , intent(in), value :: bord_num
  real    (c_double )    , intent(in)        :: b_vector(bord_num+1)
  real    (c_double )    , intent(in)        :: ee_distance_rescaled(elec_num,elec_num,walk_num)
  real    (c_double )    , intent(in)        :: asymp_jasb_pderiv(bord_num+1, 2)
  integer (c_int32_t)    , intent(in), value :: spin_independent
  real    (c_double )    , intent(out)       :: factor_ee_pderiv(bord_num+1)
  integer(qmckl_exit_code)                   :: info

  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


  factor_ee_pderiv = 0.0d0

  do nw = 1, walk_num

     do j = 1, elec_num
        do i = 1, j - 1

           x = ee_distance_rescaled(i,j,nw)

           if (spin_independent == 1) then
              factor_ee_pderiv(1) = factor_ee_pderiv(1) + x / (1.d0 + b_vector(2) * x) - asymp_jasb_pderiv(1,2)
              factor_ee_pderiv(2) = factor_ee_pderiv(2) - b_vector(1) * x**2 / (1.d0 + b_vector(2) * x)**2 - asymp_jasb_pderiv(2,2)
           else
              if ( (j <= up_num).or.(i > up_num) ) then
                 factor_ee_pderiv(1) = factor_ee_pderiv(1) + 0.5d0 * x / (1.d0 + b_vector(2) * x) - asymp_jasb_pderiv(1,1)
                 factor_ee_pderiv(2) = factor_ee_pderiv(2) - 0.5d0 * b_vector(1) * x**2 / (1.d0 + b_vector(2) * x)**2 &
                      - asymp_jasb_pderiv(2,1)
              else
                 factor_ee_pderiv(1) = factor_ee_pderiv(1) + x / (1.d0 + b_vector(2) * x) - asymp_jasb_pderiv(2,1)
                 factor_ee_pderiv(2) = factor_ee_pderiv(2) - b_vector(1) * x**2 / (1.d0 + b_vector(2) * x)**2 &
                      - asymp_jasb_pderiv(2,2)
              endif
           endif
           xk = x
           do k = 2, bord_num
              xk = xk * x
              factor_ee_pderiv(k+1) = factor_ee_pderiv(k+1) + xk - asymp_jasb_pderiv(k+1, 1)
           end do
        end do
     end do
  end do

  factor_ee_pderiv = factor_ee_pderiv / dble(walk_num)

end function qmckl_compute_jastrow_champ_factor_ee_pderiv_doc

4.1.8. Parameter Derivative of the gradient and Laplacian

The derivative of factor_ee_gl with respect to all the jastrow parameters C is computed and stored into factor_ee_gl_pderiv

4.1.8.1. Get

qmckl_exit_code
qmckl_get_jastrow_champ_factor_ee_gl_pderiv(qmckl_context context,
                                            double* const factor_ee_gl_pderiv,
                                            const int64_t size_max);

Fortran interface

interface
   integer(qmckl_exit_code) function qmckl_get_jastrow_champ_factor_ee_gl_pderiv (context, &
        factor_ee_gl_pderiv, 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
     real(c_double),   intent(out)               :: factor_ee_gl_pderiv(size_max)
   end function qmckl_get_jastrow_champ_factor_ee_gl_pderiv
end interface

4.1.8.2. Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`walk_num`	`int64_t`	in	Number of walkers
`elec_num`	`int64_t`	in	Number of electrons
`up_num`	`int64_t`	in	Number of alpha electrons
`bord_num`	`int64_t`	in	Number of coefficients
`b_vector`	`double[bord_num+1]`	in	List of coefficients
`ee_distance_rescaled`	`double[walk_num][elec_num][elec_num]`	in	Electron-electron distances
`ee_distance_rescaled_gl`	`double[walk_num][elec_num][elec_num][4]`	in	Electron-electron distances
`spin_independent`	`int32_t`	in	If 1, same parameters for parallel and antiparallel spins
`factor_ee_gl_pderiv`	`double[4][elec_num][bord_num+1]`	out	Electron-electron distances

function qmckl_compute_jastrow_champ_factor_ee_gl_pderiv_doc( &
     context, walk_num, elec_num, up_num, bord_num, &
     b_vector, ee_distance_rescaled, ee_distance_rescaled_gl,  &
     spin_independent, factor_ee_gl_pderiv) &
     bind(C) result(info)
  use qmckl
  implicit none

  integer (qmckl_context), intent(in), value :: context
  integer (c_int64_t) , intent(in)  , value :: walk_num
  integer (c_int64_t) , intent(in)  , value :: elec_num
  integer (c_int64_t) , intent(in)  , value :: up_num
  integer (c_int64_t) , intent(in)  , value :: bord_num
  real    (c_double ) , intent(in)          :: b_vector(bord_num+1)
  real    (c_double ) , intent(in)          :: ee_distance_rescaled(elec_num,elec_num,walk_num)
  real    (c_double ) , intent(in)          :: ee_distance_rescaled_gl(4,elec_num,elec_num,walk_num)
  integer (c_int32_t) , intent(in)  , value :: spin_independent
  real    (c_double ) , intent(out)         :: factor_ee_gl_pderiv(4,elec_num,bord_num+1)
  integer(qmckl_exit_code)                  :: info

  integer*8 :: i, j, k, nw, ii
  double precision   :: x, x1, kf
  double precision   :: denom, invdenom, invdenom2, invdenom3, f, f1, f2
  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

  if ((spin_independent < 0).or.(spin_independent > 1)) then
     info = QMCKL_INVALID_ARG_8
     return
  endif

  factor_ee_gl_pderiv(:,:,:) = 0.0d0
  do nw = 1, walk_num
     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
           invdenom3     = invdenom2 * 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 (spin_independent == 1) then
              f1 = invdenom2
              f2 = -2.d0 * b_vector(1) * x * invdenom3
           else
              if((i <= up_num .and. j <= up_num ) .or. (i >  up_num .and. j >  up_num)) then
                 f1 = 0.5d0 * invdenom2
                 f2 = -1.d0 * b_vector(1) * x * invdenom3
              else
                 f1 = invdenom2
                 f2 = -2.d0 * b_vector(1) * x * invdenom3
              end if
           end if

           factor_ee_gl_pderiv(1,i,1) = factor_ee_gl_pderiv(1,i,1) + f1 * dx(1)
           factor_ee_gl_pderiv(2,i,1) = factor_ee_gl_pderiv(2,i,1) + f1 * dx(2)
           factor_ee_gl_pderiv(3,i,1) = factor_ee_gl_pderiv(3,i,1) + f1 * dx(3)
           factor_ee_gl_pderiv(4,i,1) = factor_ee_gl_pderiv(4,i,1) + f1 * dx(4) &
           + 2.d0 * b_vector(2) * grad_c2 * f1 * invdenom

           factor_ee_gl_pderiv(1,i,2) = factor_ee_gl_pderiv(1,i,2) + f2 * dx(1)
           factor_ee_gl_pderiv(2,i,2) = factor_ee_gl_pderiv(2,i,2) + f2 * dx(2)
           factor_ee_gl_pderiv(3,i,2) = factor_ee_gl_pderiv(3,i,2) + f2 * dx(3)
           factor_ee_gl_pderiv(4,i,2) = factor_ee_gl_pderiv(4,i,2) &
           + 2.d0 * f1 * b_vector(1) * (-invdenom * x * dx(4) - invdenom2 * grad_c2 * (1 - 2.d0 * x * b_vector(2)))


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

        end do
     end do
  end do
  factor_ee_gl_pderiv = factor_ee_gl_pderiv / dble(walk_num)

end function qmckl_compute_jastrow_champ_factor_ee_gl_pderiv_doc

4.2. Electron-nucleus component

4.2.1. Asymptotic component

Calculate the asymptotic component asymp_jasa to be subtracted 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}} \]

4.2.1.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
     real(c_double),   intent(out)               :: asymp_jasa(size_max)
   end function qmckl_get_jastrow_champ_asymp_jasa
end interface

4.2.1.2. Compute

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

function qmckl_compute_jastrow_champ_asymp_jasa(context, aord_num, type_nucl_num, a_vector, &
     rescale_factor_en, asymp_jasa) &
     result(info) bind(C)
  use, intrinsic :: iso_c_binding
  use qmckl
  implicit none
  integer(qmckl_context), intent(in), value :: context
  integer (c_int64_t) , intent(in)  , value :: aord_num
  integer (c_int64_t) , intent(in)  , value :: type_nucl_num
  real    (c_double ) , intent(in)          :: a_vector(aord_num+1,type_nucl_num)
  real    (c_double ) , intent(in)          :: rescale_factor_en(type_nucl_num)
  real    (c_double ) , intent(out)         :: asymp_jasa(type_nucl_num)
  integer(qmckl_exit_code)                  :: info

  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

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

4.2.2.1. Get

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

4.2.2.2. Compute

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

function qmckl_compute_en_distance_rescaled_doc(context, &
     elec_num, nucl_num, type_nucl_num, &
     type_nucl_vector, rescale_factor_en, walk_num, elec_coord, &
     nucl_coord, en_distance_rescaled) &
     bind(C) result(info)
  use qmckl
  implicit none
  integer (qmckl_context), intent(in), value :: context
  integer (c_int64_t) , intent(in)  , value :: elec_num
  integer (c_int64_t) , intent(in)  , value :: nucl_num
  integer (c_int64_t) , intent(in)  , value :: type_nucl_num
  integer (c_int64_t) , intent(in)          :: type_nucl_vector(nucl_num)
  real    (c_double ) , intent(in)          :: rescale_factor_en(type_nucl_num)
  integer (c_int64_t) , intent(in)  , value :: walk_num
  real    (c_double ) , intent(in)          :: elec_coord(elec_num,walk_num,3)
  real    (c_double ) , intent(in)          :: nucl_coord(nucl_num,3)
  real    (c_double ) , intent(out)         :: en_distance_rescaled(elec_num,nucl_num,walk_num)
  integer(qmckl_exit_code)                  :: info

  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_doc

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

4.2.3.1. Get

qmckl_exit_code
qmckl_get_jastrow_champ_en_distance_rescaled_gl(qmckl_context context,
                                                double* const distance_rescaled_gl,
                                                const int64_t size_max);

4.2.3.2. Compute

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

function qmckl_compute_en_distance_rescaled_gl_doc(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) bind(C)
  use, intrinsic :: iso_c_binding
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer (c_int64_t) , intent(in)  , value :: elec_num
  integer (c_int64_t) , intent(in)  , value :: nucl_num
  integer (c_int64_t) , intent(in)  , value :: type_nucl_num
  integer (c_int64_t) , intent(in)          :: type_nucl_vector(nucl_num)
  real    (c_double ) , intent(in)          :: rescale_factor_en(nucl_num)
  integer (c_int64_t) , intent(in)  , value :: walk_num
  real    (c_double ) , intent(in)          :: elec_coord(elec_num,walk_num,3)
  real    (c_double ) , intent(in)          :: nucl_coord(nucl_num,3)
  real    (c_double ) , intent(out)         :: en_distance_rescaled_gl(4,elec_num,nucl_num,walk_num)

  integer(qmckl_exit_code)                  :: info
  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_doc

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

4.2.4.1. 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
     real(c_double),   intent(out)               :: factor_en(size_max)
   end function qmckl_get_jastrow_champ_factor_en
end interface

4.2.4.2. Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`walk_num`	`int64_t`	in	Number of walkers
`elec_num`	`int64_t`	in	Number of electrons
`nucl_num`	`int64_t`	in	Number of 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[type_nucl_num][aord_num+1]`	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

function 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) &
     bind(C) result(info)
  use qmckl
  implicit none

  integer (qmckl_context), intent(in), value :: context
  integer (c_int64_t) , intent(in)  , value :: walk_num
  integer (c_int64_t) , intent(in)  , value :: elec_num
  integer (c_int64_t) , intent(in)  , value :: nucl_num
  integer (c_int64_t) , intent(in)  , value :: type_nucl_num
  integer (c_int64_t) , intent(in)          :: type_nucl_vector(nucl_num)
  integer (c_int64_t) , intent(in)  , value :: aord_num
  real    (c_double ) , intent(in)          :: a_vector(aord_num+1,type_nucl_num)
  real    (c_double ) , intent(in)          :: en_distance_rescaled(elec_num,nucl_num,walk_num)
  real    (c_double ) , intent(in)          :: asymp_jasa(type_nucl_num)
  real    (c_double ) , intent(out)         :: factor_en(walk_num)
  integer(qmckl_exit_code)                  :: info

  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

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

The derivative is calculated in the function qmckl_compute_jastrow_champ_factor_en_gl. The formula is given by: \[ \nabla_i f_{\alpha}(R_{i\alpha}) = \sum_{j=1}^{N_\text{elec}} \left[ \frac{ A_0\, \nabla_i C_{ij} }{(1 + A_1 C_{ij})^2} + \sum_{k=2}^{N^\alpha_{\text{ord}}} A_k\, k\, C_{ij}^{k-1}\,\nabla_i C_{ij} \right] \]

\[ \Delta_i f_{\alpha}(R_{i\alpha}) = \sum_{j=1}^{N_\text{elec}} \left[ \frac{ A_0\, \Delta_i C_{ij} }{(1 + A_1 C_{ij})^2} - \frac{ 2 A_0\, A_1 (\nabla_i C_{ij})^2}{(1 + A_1 C_{ij})^3} + \sum_{k=2}^{N^\alpha_{\text{ord}}} A_k\, k\, C_{ij}^{k-1} C_{ij}^{k-2} \left[ \Delta_i C_{ij} + (k-1)(\nabla_i C_{ij})^2 \right] \right] \]

4.2.5.1. 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
     real(c_double),   intent(out)               :: factor_en_gl(size_max)
   end function qmckl_get_jastrow_champ_factor_en_gl
end interface

4.2.5.2. Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`walk_num`	`int64_t`	in	Number of walkers
`elec_num`	`int64_t`	in	Number of electrons
`nucl_num`	`int64_t`	in	Number of 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[type_nucl_num][aord_num+1]`	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][nucl_num][elec_num][4]`	in	Electron-nucleus distance derivatives
`factor_en_gl`	`double[walk_num][4][elec_num]`	out	Electron-nucleus jastrow

function qmckl_compute_jastrow_champ_factor_en_gl_doc( &
     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) &
     bind(C) result(info)
  use qmckl
  implicit none

  integer (qmckl_context), intent(in), value :: context
  integer (c_int64_t) , intent(in)  , value :: walk_num
  integer (c_int64_t) , intent(in)  , value :: elec_num
  integer (c_int64_t) , intent(in)  , value :: nucl_num
  integer (c_int64_t) , intent(in)  , value :: type_nucl_num
  integer (c_int64_t) , intent(in)          :: type_nucl_vector(nucl_num)
  integer (c_int64_t) , intent(in)  , value :: aord_num
  real    (c_double ) , intent(in)          :: a_vector(aord_num+1,type_nucl_num)
  real    (c_double ) , intent(in)          :: en_distance_rescaled(elec_num,nucl_num,walk_num)
  real    (c_double ) , intent(in)          :: en_distance_rescaled_gl(4, elec_num,nucl_num,walk_num)
  real    (c_double ) , intent(out)         :: factor_en_gl(elec_num,4,walk_num)
  integer(qmckl_exit_code)                  :: info

  integer*8 :: i, a, 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 (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_gl(:,:,nw) = 0.0d0
     do a = 1, nucl_num
        do i = 1, elec_num

           x = en_distance_rescaled(i,a,nw)
           if(abs(x) < 1.d-12) continue

           denom = 1.0d0 + a_vector(2, type_nucl_vector(a)+1) * x
           invdenom = 1.0d0 / denom
           invdenom2 = invdenom*invdenom

           dx(1) = en_distance_rescaled_gl(1,i,a,nw)
           dx(2) = en_distance_rescaled_gl(2,i,a,nw)
           dx(3) = en_distance_rescaled_gl(3,i,a,nw)
           dx(4) = en_distance_rescaled_gl(4,i,a,nw)

           f = a_vector(1, type_nucl_vector(a)+1) * invdenom2
           grad_c2 = dx(1)*dx(1) + dx(2)*dx(2) + dx(3)*dx(3)

           factor_en_gl(i,1,nw) = factor_en_gl(i,1,nw) + f * dx(1)
           factor_en_gl(i,2,nw) = factor_en_gl(i,2,nw) + f * dx(2)
           factor_en_gl(i,3,nw) = factor_en_gl(i,3,nw) + f * dx(3)
           factor_en_gl(i,4,nw) = factor_en_gl(i,4,nw) &
                + f * (dx(4) - 2.d0 * a_vector(2, type_nucl_vector(a)+1) * grad_c2 * invdenom)


           kf = 2.d0
           x1 = x
           x = 1.d0
           do k=2, aord_num
              f = a_vector(k+1,type_nucl_vector(a)+1) * kf * x
              factor_en_gl(i,1,nw) = factor_en_gl(i,1,nw) + f * x1 * dx(1)
              factor_en_gl(i,2,nw) = factor_en_gl(i,2,nw) + f * x1 * dx(2)
              factor_en_gl(i,3,nw) = factor_en_gl(i,3,nw) + f * x1 * dx(3)
              factor_en_gl(i,4,nw) = factor_en_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_en_gl_doc

4.2.5.3. Test

4.2.6. Parameter Derivative of the Asymptotic component

Calculate the derivative of the asymptotic component asymp_jasa to be subtracted from the final electron-nucleus jastrow factor parameter derivatives \(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} = \frac{A_0 \kappa_\text{en}^{-1}}{1 + A_1\, \kappa_\text{en}^{-1}} + \sum_{p=2}^{N_\text{ord}^B} A_{p+1}\, \kappa_\text{en}^{-p} \]

\[ \partial_{B_0} J_{\text{ee}}^\infty = \frac{\frac{1}{2}(1 + \delta^{\uparrow\downarrow}_{ij}) \kappa_{\text{ee}}^{-1}}{1 + B_1 \kappa_{\text{ee}}^{-1}} \]

\[ \partial_{B_1}J_{\text{ee}}^\infty = -\frac{ \frac{1}{2}(1 + \delta_{ij}) \kappa_{\text{ee}}^{-2}} {(1 + B_1 \kappa_{\text{ee}}^{-1})^2 } \]

\[ \partial_{B_k} J_{\text{ee}}^\infty = \kappa_{\text{ee}}^{-(k-1)} \] for \(k > 1\)

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

4.2.6.1. Get

qmckl_exit_code
qmckl_get_jastrow_champ_asymp_jasa_pderiv(qmckl_context context,
                                          double* const asymp_jasa_pderiv,
                                          const int64_t size_max);

Fortran interface

interface
   integer(qmckl_exit_code) function qmckl_get_jastrow_champ_asymp_jasa_pderiv(context, &
        asymp_jasa_pderiv, 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
     real(c_double),   intent(out)               :: asymp_jasa_pderiv(size_max)
   end function qmckl_get_jastrow_champ_asymp_jasa_pderiv
end interface

4.2.6.2. Compute

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

function qmckl_compute_jastrow_champ_asymp_jasa_pderiv(context, aord_num, type_nucl_num, a_vector, &
     rescale_factor_en, asymp_jasa_pderiv) &
     result(info) bind(C)
  use, intrinsic :: iso_c_binding
  use qmckl
  implicit none
  integer(qmckl_context), intent(in), value :: context
  integer (c_int64_t) , intent(in)  , value :: aord_num
  integer (c_int64_t) , intent(in)  , value :: type_nucl_num
  real    (c_double ) , intent(in)          :: a_vector(aord_num+1,type_nucl_num)
  real    (c_double ) , intent(in)          :: rescale_factor_en(type_nucl_num)
  real    (c_double ) , intent(out)         :: asymp_jasa_pderiv(aord_num+1, type_nucl_num)
  integer(qmckl_exit_code)                  :: info

  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_pderiv(1,i) = kappa_inv / (1.0d0 + a_vector(2,i) * kappa_inv)
     asymp_jasa_pderiv(2,i) = - a_vector(1,i) * kappa_inv**2 / (1.0d0 + a_vector(2,i) * kappa_inv)**2

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

  end do

end function qmckl_compute_jastrow_champ_asymp_jasa_pderiv

4.2.7. Parameter Derivative

Calculate the parameter derivatives of the electron-nucleus jastrow component factor_en_pderiv using the a_vector coeffecients and the electron-nucleus rescaled distances en_distance_rescaled. The result is stored in factor_en_pderiv. If \(\alpha\) only runs over atoms of the same type then the expressions for the derivatives are: \[ \frac{\partial f}{\partial A_0} = \sum_{i, \alpha} \frac{ R_{i\alpha} }{1 + A_1 R_{i\alpha} } - \frac{\partial J_{en}^\infty}{\partial A_0} \] \[ \frac{\partial f}{\partial A_1} = \sum_{i, \alpha} \frac{ A_0 R_{i\alpha}^2 }{(1 + A_1 R_{i\alpha})^2 } - \frac{\partial J_{en}^\infty}{\partial A_1} \] \[ \frac{\partial f}{\partial A_k} = \sum_{i, \alpha} R_{i\alpha}^k - \frac{\partial J_{en}^\infty}{\partial A_k} \]

4.2.7.1. Get

qmckl_exit_code
qmckl_get_jastrow_champ_factor_en_pderiv(qmckl_context context,
                                         double* const factor_en_pderiv,
                                         const int64_t size_max);

interface
   integer(qmckl_exit_code) function qmckl_get_jastrow_champ_factor_en_pderiv (context, &
        factor_en_pderiv, 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
     real(c_double),   intent(out)               :: factor_en_pderiv(size_max)
   end function qmckl_get_jastrow_champ_factor_en_pderiv
end interface

4.2.7.2. Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`walk_num`	`int64_t`	in	Number of walkers
`elec_num`	`int64_t`	in	Number of electrons
`nucl_num`	`int64_t`	in	Number of 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[type_nucl_num][aord_num+1]`	in	List of coefficients
`en_distance_rescaled`	`double[walk_num][nucl_num][elec_num]`	in	Electron-nucleus distances
`asymp_jasa_pderiv`	`double[type_nucl_num][aord_num+1]`	in	Type of nuclei
`factor_en_pderiv`	`double[type_nucl_num][aord_num+1]`	out	Electron-nucleus jastrow

function qmckl_compute_jastrow_champ_factor_en_pderiv_doc( &
     context, walk_num, elec_num, nucl_num, type_nucl_num, &
     type_nucl_vector, aord_num, a_vector, &
     en_distance_rescaled, asymp_jasa_pderiv, factor_en_pderiv) &
     bind(C) result(info)
  use qmckl
  implicit none

  integer (qmckl_context), intent(in), value :: context
  integer (c_int64_t) , intent(in)  , value :: walk_num
  integer (c_int64_t) , intent(in)  , value :: elec_num
  integer (c_int64_t) , intent(in)  , value :: nucl_num
  integer (c_int64_t) , intent(in)  , value :: type_nucl_num
  integer (c_int64_t) , intent(in)          :: type_nucl_vector(nucl_num)
  integer (c_int64_t) , intent(in)  , value :: aord_num
  real    (c_double ) , intent(in)          :: a_vector(aord_num+1,type_nucl_num)
  real    (c_double ) , intent(in)          :: en_distance_rescaled(elec_num,nucl_num,walk_num)
  real    (c_double ) , intent(in)          :: asymp_jasa_pderiv(aord_num+1,type_nucl_num)
  real    (c_double ) , intent(out)         :: factor_en_pderiv(aord_num+1,type_nucl_num)
  integer(qmckl_exit_code)                  :: info

  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

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

           factor_en_pderiv(1,type_nucl_vector(a)+1) = factor_en_pderiv(1,type_nucl_vector(a)+1) + & 
           x / (1.0d0 + a_vector(2, type_nucl_vector(a)+1) * x) - asymp_jasa_pderiv(1,type_nucl_vector(a)+1)

           factor_en_pderiv(2,type_nucl_vector(a)+1) = factor_en_pderiv(2,type_nucl_vector(a)+1) - &
           a_vector(2, type_nucl_vector(a)+1) * x**2 / (1.0d0 + a_vector(2, type_nucl_vector(a)+1) * x)**2 &
           - asymp_jasa_pderiv(2,type_nucl_vector(a)+1)

           do p = 2, aord_num
              x = x * en_distance_rescaled(i, a, nw)
              factor_en_pderiv(p+1,type_nucl_vector(a)+1) = factor_en_pderiv(p+1,type_nucl_vector(a)+1) + x &
              - asymp_jasa_pderiv(p+1,type_nucl_vector(a)+1)
           end do
        end do
     end do
  end do
  factor_en_pderiv = factor_en_pderiv / dble(walk_num)
end function qmckl_compute_jastrow_champ_factor_en_pderiv_doc

4.2.8. Parameter Derivative of the gradient and Laplacian

The derivative of factor_en_gl with respect to all the jastrow parameters C is computed and stored into factor_en_gl_pderiv

4.2.8.1. Get

qmckl_exit_code
qmckl_get_jastrow_champ_factor_en_gl_pderiv(qmckl_context context,
                                    double* const factor_en_gl_pderiv,
                                    const int64_t size_max);

Fortran interface

interface
   integer(qmckl_exit_code) function qmckl_get_jastrow_champ_factor_en_gl_pderiv (context, &
        factor_en_gl_pderiv, 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
     real(c_double),   intent(out)               :: factor_en_gl_pderiv(size_max)
   end function qmckl_get_jastrow_champ_factor_en_gl_pderiv
end interface

4.2.8.2. Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`walk_num`	`int64_t`	in	Number of walkers
`elec_num`	`int64_t`	in	Number of electrons
`nucl_num`	`int64_t`	in	Number of 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[type_nucl_num][aord_num+1]`	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][nucl_num][elec_num][4]`	in	Electron-nucleus distance derivatives
`factor_en_gl_pderiv`	`double[type_nucl_num][aord_num+1][elec_num][4]`	out	Electron-nucleus jastrow

function qmckl_compute_jastrow_champ_factor_en_gl_pderiv_doc( &
     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_pderiv) &
     bind(C) result(info)
  use qmckl
  implicit none

  integer (qmckl_context), intent(in), value :: context
  integer (c_int64_t) , intent(in)  , value :: walk_num
  integer (c_int64_t) , intent(in)  , value :: elec_num
  integer (c_int64_t) , intent(in)  , value :: nucl_num
  integer (c_int64_t) , intent(in)  , value :: type_nucl_num
  integer (c_int64_t) , intent(in)          :: type_nucl_vector(nucl_num)
  integer (c_int64_t) , intent(in)  , value :: aord_num
  real    (c_double ) , intent(in)          :: a_vector(aord_num+1,type_nucl_num)
  real    (c_double ) , intent(in)          :: en_distance_rescaled(elec_num,nucl_num,walk_num)
  real    (c_double ) , intent(in)          :: en_distance_rescaled_gl(4, elec_num,nucl_num,walk_num)
  real    (c_double ) , intent(out)         :: factor_en_gl_pderiv(4,elec_num,aord_num+1,type_nucl_num)
  integer(qmckl_exit_code)                  :: info

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

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

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

  factor_en_gl_pderiv = 0.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.d-12) continue

           denom = 1.0d0 + a_vector(2, type_nucl_vector(a)+1) * x
           invdenom = 1.0d0 / denom
           invdenom2 = invdenom*invdenom
           invdenom3 = invdenom2*invdenom

           dx(1) = en_distance_rescaled_gl(1,i,a,nw)
           dx(2) = en_distance_rescaled_gl(2,i,a,nw)
           dx(3) = en_distance_rescaled_gl(3,i,a,nw)
           dx(4) = en_distance_rescaled_gl(4,i,a,nw)

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

           factor_en_gl_pderiv(1,i,1,type_nucl_vector(a)+1) = factor_en_gl_pderiv(1,i,1,type_nucl_vector(a)+1) + f * dx(1)
           factor_en_gl_pderiv(2,i,1,type_nucl_vector(a)+1) = factor_en_gl_pderiv(2,i,1,type_nucl_vector(a)+1) + f * dx(2)
           factor_en_gl_pderiv(3,i,1,type_nucl_vector(a)+1) = factor_en_gl_pderiv(3,i,1,type_nucl_vector(a)+1) + f * dx(3)
           factor_en_gl_pderiv(4,i,1,type_nucl_vector(a)+1) = factor_en_gl_pderiv(4,i,1,type_nucl_vector(a)+1) &
                + f * (dx(4) - 2.d0 * a_vector(2, type_nucl_vector(a)+1) * grad_c2 * invdenom)

           f = -2.d0 * a_vector(1, type_nucl_vector(a)+1) * invdenom3

           factor_en_gl_pderiv(1,i,2,type_nucl_vector(a)+1) = factor_en_gl_pderiv(1,i,2,type_nucl_vector(a)+1) + f * x * dx(1)
           factor_en_gl_pderiv(2,i,2,type_nucl_vector(a)+1) = factor_en_gl_pderiv(2,i,2,type_nucl_vector(a)+1) + f * x * dx(2)
           factor_en_gl_pderiv(3,i,2,type_nucl_vector(a)+1) = factor_en_gl_pderiv(3,i,2,type_nucl_vector(a)+1) + f * x * dx(3)
           factor_en_gl_pderiv(4,i,2,type_nucl_vector(a)+1) = factor_en_gl_pderiv(4,i,2,type_nucl_vector(a)+1) &
                + f * (x * dx(4) + invdenom * (grad_c2 - 2 * grad_c2 * x * a_vector(2, type_nucl_vector(a)+1)))

           kf = 2.d0
           x1 = x
           x = 1.d0
           do k=3, aord_num+1
              f = kf * x
              factor_en_gl_pderiv(1,i,k,type_nucl_vector(a)+1) = factor_en_gl_pderiv(1,i,k,type_nucl_vector(a)+1) + f * x1 * dx(1)
              factor_en_gl_pderiv(2,i,k,type_nucl_vector(a)+1) = factor_en_gl_pderiv(2,i,k,type_nucl_vector(a)+1) + f * x1 * dx(2)
              factor_en_gl_pderiv(3,i,k,type_nucl_vector(a)+1) = factor_en_gl_pderiv(3,i,k,type_nucl_vector(a)+1) + f * x1 * dx(3)
              factor_en_gl_pderiv(4,i,k,type_nucl_vector(a)+1) = factor_en_gl_pderiv(4,i,k,type_nucl_vector(a)+1) &
                   + f * (x1 * dx(4) + (kf-1.d0) * grad_c2)
              x = x*x1
              kf = kf + 1.d0
           end do

        end do
     end do
  end do
  factor_en_gl_pderiv = factor_en_gl_pderiv / dble(walk_num)
end function qmckl_compute_jastrow_champ_factor_en_gl_pderiv_doc

4.3. Electron-electron-nucleus component

4.3.1. 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:

\[ [g_e(r)_{ij}]^p = \begin{cases} \exp\left(-p\,\kappa_\text{e}\, r_{ij}\right) & \text{for } i \ne j \\ 0 & \text{for } i = j \]

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

4.3.1.1. Get

qmckl_exit_code
qmckl_get_jastrow_champ_een_rescaled_e(qmckl_context context,
                                                double* const een_rescaled_e,
                                                const int64_t size_max);

4.3.1.2. Compute

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

function qmckl_compute_een_rescaled_e_doc( &
     context, walk_num, elec_num, cord_num, rescale_factor_ee,  &
     ee_distance, een_rescaled_e) &
     result(info) bind(C)
  use, intrinsic :: iso_c_binding
  use qmckl
  implicit none
  integer(qmckl_context), intent(in), value :: context
  integer (c_int64_t) , intent(in)  , value :: walk_num
  integer (c_int64_t) , intent(in)  , value :: elec_num
  integer (c_int64_t) , intent(in)  , value :: cord_num
  real    (c_double ) , intent(in)  , value :: rescale_factor_ee
  real    (c_double ) , intent(in)          :: ee_distance(elec_num,elec_num,walk_num)
  real    (c_double ) , intent(out)         :: een_rescaled_e(elec_num,elec_num,0:cord_num,walk_num)
  integer(qmckl_exit_code)                  :: info

  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

  do nw = 1, walk_num
     do l = 0, cord_num
        do j = 1, elec_num
           do i = 1, elec_num
              een_rescaled_e(i, j, l, nw) = dexp(-rescale_factor_ee * ee_distance(i, j, nw))**l
           end do
           een_rescaled_e(j, j, l, nw) = 0.d0
        end do
     end do
  end do

end function qmckl_compute_een_rescaled_e_doc

4.3.1.3. Test

4.3.2. 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 = \kappa_\text{e}\, p\,\left[ - \frac{2}{r} + \kappa_\text{e}\, p \right] \left[ {g_\text{e}(r)} \right]^p \]

Derivatives are set to zero at \(r_{ii}\) to avoid NaNs.

4.3.2.1. Get

qmckl_exit_code
qmckl_get_jastrow_champ_een_rescaled_e_gl(qmckl_context context,
                                          double* const een_rescaled_e_gl,
                                          const int64_t size_max);

4.3.2.2. Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`walk_num`	`int64_t`	in	Number of walkers
`elec_num`	`int64_t`	in	Number of electrons
`cord_num`	`int64_t`	in	Order of polynomials
`rescale_factor_ee`	`double`	in	Factor to rescale ee distances
`coord_ee`	`double[3][walk_num][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

function qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl_doc( &
     context, walk_num, elec_num, cord_num, rescale_factor_ee,  &
     coord_ee, ee_distance, een_rescaled_e, een_rescaled_e_gl) &
     result(info) bind(C)
  use, intrinsic :: iso_c_binding
  use qmckl
  implicit none
  integer(qmckl_context), intent(in), value  :: context
  integer(c_int64_t)    , intent(in), value  :: walk_num
  integer(c_int64_t)    , intent(in), value  :: elec_num
  integer(c_int64_t)    , intent(in), value  :: cord_num
  real(c_double)        , intent(in), value  :: rescale_factor_ee
  real(c_double)        , intent(in)  :: coord_ee(elec_num,walk_num,3)
  real(c_double)        , intent(in)  :: ee_distance(elec_num,elec_num,walk_num)
  real(c_double)        , intent(in)  :: een_rescaled_e(elec_num,elec_num,0:cord_num,walk_num)
  real(c_double)        , intent(out) :: een_rescaled_e_gl(elec_num,4,elec_num,0:cord_num,walk_num)
  integer(qmckl_exit_code)            :: info

  double precision, allocatable       :: elec_dist_gl(:,:,:)
  double precision                    :: x, kappa_l
  integer*8                           :: i, j, k, l, nw, ii

  double precision  :: rij_inv(elec_num)


  allocate(elec_dist_gl(elec_num, 4, elec_num))
  elec_dist_gl = 0.d0
  een_rescaled_e_gl = 0.d0

  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

  do nw = 1, walk_num

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

     !   Not necessary: should be set to zero by qmckl_malloc
     een_rescaled_e_gl(:,:,:,0,nw) = 0.d0

     do l = 1, cord_num
        kappa_l = -dble(l) * rescale_factor_ee
        do j = 1, elec_num
           do i = 1, elec_num
              if (i /= j) then
                 een_rescaled_e_gl(i, 1, j, l, nw) = kappa_l *  elec_dist_gl(i, 1, j) * een_rescaled_e(i,j,l,nw)
                 een_rescaled_e_gl(i, 2, j, l, nw) = kappa_l *  elec_dist_gl(i, 2, j) * een_rescaled_e(i,j,l,nw)
                 een_rescaled_e_gl(i, 3, j, l, nw) = kappa_l *  elec_dist_gl(i, 3, j) * een_rescaled_e(i,j,l,nw)
                 een_rescaled_e_gl(i, 4, j, l, nw) = kappa_l * (elec_dist_gl(i, 4, j) + kappa_l) * een_rescaled_e(i,j,l,nw)
              else
                 een_rescaled_e_gl(i, 1, j, l, nw) = 0.d0
                 een_rescaled_e_gl(i, 2, j, l, nw) = 0.d0
                 een_rescaled_e_gl(i, 3, j, l, nw) = 0.d0
                 een_rescaled_e_gl(i, 4, j, l, nw) = 0.d0
              end if
           end do
        end do
     end do
  end do

end function qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl_doc

4.3.3. 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:

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

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

4.3.3.1. Get

qmckl_exit_code
qmckl_get_jastrow_champ_een_rescaled_n(qmckl_context context,
                                 double* const een_rescaled_n,
                                 const int64_t size_max);

4.3.3.2. Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`walk_num`	`int64_t`	in	Number of walkers
`elec_num`	`int64_t`	in	Number of electrons
`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

function qmckl_compute_een_rescaled_n( &
     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) bind(C)
  use, intrinsic :: iso_c_binding
  use qmckl
  implicit none
  integer(qmckl_context), intent(in), value :: context
  integer (c_int64_t) , intent(in)  , value :: walk_num
  integer (c_int64_t) , intent(in)  , value :: elec_num
  integer (c_int64_t) , intent(in)  , value :: nucl_num
  integer (c_int64_t) , intent(in)  , value :: type_nucl_num
  integer (c_int64_t) , intent(in)          :: type_nucl_vector(nucl_num)
  integer (c_int64_t) , intent(in)  , value :: cord_num
  real    (c_double ) , intent(in)          :: rescale_factor_en(nucl_num)
  real    (c_double ) , intent(in)          :: en_distance(nucl_num,elec_num,walk_num)
  real    (c_double ) , intent(out)         :: een_rescaled_n(elec_num,nucl_num,0:cord_num,walk_num)
  integer(qmckl_exit_code)                  :: info

  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

  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

4.3.4. 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 \]

4.3.4.1. Get

qmckl_exit_code
qmckl_get_jastrow_champ_een_rescaled_n_gl(qmckl_context context,
                                         double* const een_rescaled_n_gl,
                                         const int64_t size_max);

4.3.4.2. Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`walk_num`	`int64_t`	in	Number of walkers
`elec_num`	`int64_t`	in	Number of electrons
`nucl_num`	`int64_t`	in	Number of atoms
`type_nucl_num`	`int64_t`	in	Number of atom types
`type_nucl_vector`	`int64_t[nucl_num]`	in	Types of atoms
`cord_num`	`int64_t`	in	Order of polynomials
`rescale_factor_en`	`double[nucl_num]`	in	Factor to rescale ee distances
`coord_ee`	`double[3][walk_num][elec_num]`	in	Electron coordinates
`coord_n`	`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

function qmckl_compute_jastrow_champ_factor_een_rescaled_n_gl( &
     context, walk_num, elec_num, nucl_num, type_nucl_num, type_nucl_vector, &
     cord_num, rescale_factor_en, &
     coord_ee, coord_n, en_distance, een_rescaled_n, een_rescaled_n_gl) &
     result(info) bind(C)
  use, intrinsic :: iso_c_binding
  use qmckl
  implicit none
  integer(qmckl_context), intent(in), value :: context
  integer (c_int64_t) , intent(in)  , value :: walk_num
  integer (c_int64_t) , intent(in)  , value :: elec_num
  integer (c_int64_t) , intent(in)  , value :: nucl_num
  integer (c_int64_t) , intent(in)  , value :: type_nucl_num
  integer (c_int64_t) , intent(in)          :: type_nucl_vector(nucl_num)
  integer (c_int64_t) , intent(in)  , value :: cord_num
  real    (c_double ) , intent(in)          :: rescale_factor_en(nucl_num)
  real    (c_double ) , intent(in)          :: coord_ee(elec_num,walk_num,3)
  real    (c_double ) , intent(in)          :: coord_n(nucl_num,3)
  real    (c_double ) , intent(in)          :: en_distance(nucl_num,elec_num,walk_num)
  real    (c_double ) , intent(in)          :: een_rescaled_n(elec_num,nucl_num,0:cord_num,walk_num)
  real    (c_double ) , intent(out)         :: een_rescaled_n_gl(elec_num,4,nucl_num,0:cord_num,walk_num)
  integer(qmckl_exit_code)                  :: info

  double precision,allocatable        :: elnuc_dist_gl(:,:,:)
  double precision                    :: x, ria_inv, kappa_l
  integer*8                           :: i, a, k, l, nw, ii

  allocate(elnuc_dist_gl(elec_num, 4, 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(i, ii, a) = (coord_ee(i, nw, ii) - coord_n(a, ii)) * ria_inv
           end do
           elnuc_dist_gl(i, 4, 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(i, 1, a)
              een_rescaled_n_gl(i, 2, a, l, nw) = kappa_l * elnuc_dist_gl(i, 2, a)
              een_rescaled_n_gl(i, 3, a, l, nw) = kappa_l * elnuc_dist_gl(i, 3, a)
              een_rescaled_n_gl(i, 4, a, l, nw) = kappa_l * elnuc_dist_gl(i, 4, 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

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

The array tmp_c corresponds to the tensor \(P\) defined at the beginning of this section:

\[ P_{i\alpha}^{km} = \sum_{j=1}^{N_{\text{elec}}} \left[ g_\text{e}({r}_{ij}) \right]^k \left[ g_\alpha({R}_{j\alpha}) \right]^{m} \]

\[ \nabla_i P_{i\alpha}^{km} = \sum_{j=1}^{N_{\text{elec}}} \nabla_i \left[ g_\text{e}({r}_{ij}) \right]^k \left[ g_\alpha({R}_{j\alpha}) \right]^{m} \]

4.3.5.1. Compute dim_c_vector

Computes the dimension of the vector of parameters.

\(N_{ord}\)	Number of parameters

1	0
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

function qmckl_compute_dim_c_vector_doc( &
     context, cord_num, dim_c_vector) &
     result(info) bind(C)
  use, intrinsic :: iso_c_binding
  use qmckl
  implicit none
  integer(qmckl_context), intent(in), value :: context
  integer (c_int64_t) , intent(in)  , value :: cord_num
  integer (c_int64_t) , intent(out)         :: dim_c_vector
  integer(qmckl_exit_code)                  :: info

  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_doc

4.3.5.2. Get

4.3.5.3. Compute c_vector_full

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

function qmckl_compute_c_vector_full_doc( &
     context, nucl_num, dim_c_vector, type_nucl_num,  &
     type_nucl_vector, c_vector, c_vector_full) &
     result(info) bind(C)
  use, intrinsic :: iso_c_binding
  use qmckl
  implicit none
  integer(qmckl_context), intent(in), value :: context
  integer (c_int64_t) , intent(in)  , value :: nucl_num
  integer (c_int64_t) , intent(in)  , value :: dim_c_vector
  integer (c_int64_t) , intent(in)  , value :: type_nucl_num
  integer (c_int64_t) , intent(in)          :: type_nucl_vector(nucl_num)
  real    (c_double ) , intent(in)          :: c_vector(dim_c_vector, type_nucl_num)
  real    (c_double ) , intent(out)         :: c_vector_full(nucl_num,dim_c_vector)
  integer(qmckl_exit_code)                  :: info

  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(1:dim_c_vector, type_nucl_vector(a)+1)
  end do

end function qmckl_compute_c_vector_full_doc

4.3.5.4. Compute lkpm_combined_index

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

function qmckl_compute_lkpm_combined_index_doc( &
     context, cord_num, dim_c_vector,  lkpm_combined_index) &
     result(info) bind(C)
  use, intrinsic :: iso_c_binding
  use qmckl
  implicit none
  integer(qmckl_context), intent(in), value :: context
  integer (c_int64_t) , intent(in)  , value :: cord_num
  integer (c_int64_t) , intent(in)  , value :: dim_c_vector
  integer (c_int64_t) , intent(out)         :: lkpm_combined_index(dim_c_vector,4)
  integer(qmckl_exit_code)                  :: info

  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 /= 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_8) 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_doc

4.3.5.5. Compute tmp_c

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`cord_num`	`int64_t`	in	Order of polynomials
`elec_num`	`int64_t`	in	Number of electrons
`nucl_num`	`int64_t`	in	Number of 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

function qmckl_compute_tmp_c_doc( &
     context, cord_num, elec_num, nucl_num, &
     walk_num, een_rescaled_e, een_rescaled_n, tmp_c) &
     result(info) bind(C)
  use, intrinsic :: iso_c_binding
  use qmckl
  implicit none
  integer(qmckl_context), intent(in), value  :: context
  integer (c_int64_t) , intent(in)  , value :: cord_num
  integer (c_int64_t) , intent(in)  , value :: elec_num
  integer (c_int64_t) , intent(in)  , value :: nucl_num
  integer (c_int64_t) , intent(in)  , value :: walk_num
  real    (c_double ) , intent(in)          :: een_rescaled_e(elec_num,elec_num,0:cord_num,walk_num)
  real    (c_double ) , intent(in)          :: een_rescaled_n(elec_num,nucl_num,0:cord_num,walk_num)
  real    (c_double ) , intent(out)         :: tmp_c(elec_num,nucl_num,0:cord_num,0:cord_num-1,walk_num)
  integer(qmckl_exit_code)                  :: info

  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
  K = elec_num
  N = nucl_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
        do j=0, cord_num-i
           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,j,nw),LDB*1_8,        &
                beta,                                    &
                tmp_c(1,1,j,i,nw),LDC)
        end do
     end do
  end do

end function qmckl_compute_tmp_c_doc

4.3.5.6. Compute dtmp_c

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`cord_num`	`int64_t`	in	Order of polynomials
`elec_num`	`int64_t`	in	Number of electrons
`nucl_num`	`int64_t`	in	Number of 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][4][elec_num]`	out	vector of non-zero coefficients

function qmckl_compute_dtmp_c_doc( &
     context, cord_num, elec_num, nucl_num, &
     walk_num, een_rescaled_e_gl, een_rescaled_n, dtmp_c) &
     result(info) bind(C)
  use, intrinsic :: iso_c_binding
  use qmckl
  implicit none
  integer(qmckl_context), intent(in), value :: context
  integer (c_int64_t) , intent(in)  , value :: cord_num
  integer (c_int64_t) , intent(in)  , value :: elec_num
  integer (c_int64_t) , intent(in)  , value :: nucl_num
  integer (c_int64_t) , intent(in)  , value :: walk_num
  real    (c_double ) , intent(in)          :: een_rescaled_e_gl(elec_num,4,elec_num,0:cord_num,walk_num)
  real    (c_double ) , intent(in)          :: een_rescaled_n(elec_num,nucl_num,0:cord_num,walk_num)
  real    (c_double ) , intent(out)         :: dtmp_c(elec_num,4,nucl_num,0:cord_num,0:cord_num-1,walk_num)
  integer(qmckl_exit_code)                  :: info

  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
  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
        do j=0, cord_num-i
              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,j,nw),LDB*1_8,                      &
                   beta,                                                  &
                   dtmp_c(1,1,1,j,i,nw),LDC)
        end do
     end do
  end do

end function qmckl_compute_dtmp_c_doc

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

4.3.6.1. 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
     real(c_double),   intent(out)               :: factor_een(size_max)
   end function qmckl_get_jastrow_champ_factor_een
end interface

4.3.6.2. Compute naive

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`walk_num`	`int64_t`	in	Number of walkers
`elec_num`	`int64_t`	in	Number of electrons
`nucl_num`	`int64_t`	in	Number of 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][0:cord_num][elec_num][elec_num]`	in	Electron-nucleus rescaled
`een_rescaled_n`	`double[walk_num][0:cord_num][nucl_num][elec_num]`	in	Electron-nucleus rescaled factor
`factor_een`	`double[walk_num]`	out	Electron-nucleus jastrow

function qmckl_compute_jastrow_champ_factor_een_naive( &
     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) bind(C)
  use, intrinsic :: iso_c_binding
  use qmckl
  implicit none
  integer(qmckl_context), intent(in), value :: context
  integer (c_int64_t) , intent(in)  , value :: walk_num
  integer (c_int64_t) , intent(in)  , value :: elec_num
  integer (c_int64_t) , intent(in)  , value :: nucl_num
  integer (c_int64_t) , intent(in)  , value :: cord_num
  integer (c_int64_t) , intent(in)  , value :: dim_c_vector
  real    (c_double ) , intent(in)          :: c_vector_full(nucl_num,dim_c_vector)
  integer (c_int64_t) , intent(in)          :: lkpm_combined_index(dim_c_vector,4)
  real    (c_double ) , intent(in)          :: een_rescaled_e(elec_num,elec_num,0:cord_num,walk_num)
  real    (c_double ) , intent(in)          :: een_rescaled_n(elec_num,nucl_num,0:cord_num,walk_num)
  real    (c_double ) , intent(out)         :: factor_een(walk_num)
  integer(qmckl_exit_code)                  :: info

  integer*8 :: i, a, j, l, k, m, n, p, nw
  double precision :: 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

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

        do a = 1, nucl_num
           cn = c_vector_full(a, n)
           if (cn == 0.d0) cycle
           do j = 1, elec_num
              do i = 1, j-1
                 factor_een(nw) = factor_een(nw) + cn*( &
                      een_rescaled_e(i,j,k,nw) *       &
                      (een_rescaled_n(i,a,l,nw) + een_rescaled_n(j,a,l,nw)) * &
                      (een_rescaled_n(i,a,m,nw) * een_rescaled_n(j,a,m,nw)) )
              end do
           end do
        end do
     end do
  end do

end function qmckl_compute_jastrow_champ_factor_een_naive

4.3.6.3. Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`walk_num`	`int64_t`	in	Number of walkers
`elec_num`	`int64_t`	in	Number of electrons
`nucl_num`	`int64_t`	in	Number of 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

function qmckl_compute_jastrow_champ_factor_een_doc( &
     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) bind(C)
  use, intrinsic :: iso_c_binding
  use qmckl
  implicit none
  integer(qmckl_context), intent(in), value :: context
  integer (c_int64_t) , intent(in)  , value :: walk_num
  integer (c_int64_t) , intent(in)  , value :: elec_num
  integer (c_int64_t) , intent(in)  , value :: nucl_num
  integer (c_int64_t) , intent(in)  , value :: cord_num
  integer (c_int64_t) , intent(in)  , value :: dim_c_vector
  real    (c_double ) , intent(in)          :: c_vector_full(nucl_num,dim_c_vector)
  integer (c_int64_t) , intent(in)          :: lkpm_combined_index(dim_c_vector,4)
  real    (c_double ) , intent(in)          :: tmp_c(elec_num,nucl_num,0:cord_num,0:cord_num-1,walk_num)
  real    (c_double ) , intent(in)          :: een_rescaled_n(elec_num,nucl_num,0:cord_num,walk_num)
  real    (c_double ) , intent(out)         :: factor_een(walk_num)
  integer(qmckl_exit_code)                  :: info

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

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

4.3.7.1. Get

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

qmckl_exit_code
qmckl_get_jastrow_champ_factor_een_grad(qmckl_context context,
                                        double* const factor_een_grad,
                                        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
     real(c_double),   intent(out)               :: factor_een_gl(size_max)
   end function qmckl_get_jastrow_champ_factor_een_gl
end interface

interface
   integer(qmckl_exit_code) function qmckl_get_jastrow_champ_factor_een_grad (context, &
        factor_een_grad, 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
     real(c_double),   intent(out)               :: factor_een_grad(size_max)
   end function qmckl_get_jastrow_champ_factor_een_grad
end interface

4.3.7.2. Compute Naive

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`walk_num`	`int64_t`	in	Number of walkers
`elec_num`	`int64_t`	in	Number of electrons
`nucl_num`	`int64_t`	in	Number of 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][0:cord_num][elec_num][elec_num]`	in	Electron-nucleus rescaled
`een_rescaled_n`	`double[walk_num][0:cord_num][nucl_num][elec_num]`	in	Electron-nucleus rescaled factor
`een_rescaled_e_gl`	`double[walk_num][0:cord_num][elec_num][4][elec_num]`	in	Electron-nucleus rescaled
`een_rescaled_n_gl`	`double[walk_num][0:cord_num][nucl_num][4][elec_num]`	in	Electron-nucleus rescaled factor
`factor_een_gl`	`double[walk_num][4][elec_num]`	out	Electron-nucleus jastrow

function qmckl_compute_jastrow_champ_factor_een_gl_naive( &
     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) bind(C)
  use, intrinsic :: iso_c_binding
  use qmckl
  implicit none
  integer(qmckl_context), intent(in), value :: context
  integer (c_int64_t) , intent(in)  , value :: walk_num
  integer (c_int64_t) , intent(in)  , value :: elec_num
  integer (c_int64_t) , intent(in)  , value :: nucl_num
  integer (c_int64_t) , intent(in)  , value :: cord_num
  integer (c_int64_t) , intent(in)  , value :: dim_c_vector
  real    (c_double ) , intent(in)          :: c_vector_full(nucl_num,dim_c_vector)
  integer (c_int64_t) , intent(in)          :: lkpm_combined_index(dim_c_vector,4)
  real    (c_double ) , intent(in)          :: een_rescaled_e(elec_num,elec_num,0:cord_num,walk_num)
  real    (c_double ) , intent(in)          :: een_rescaled_n(elec_num,nucl_num,0:cord_num,walk_num)
  real    (c_double ) , intent(in)          :: een_rescaled_e_gl(elec_num,4,elec_num,0:cord_num,walk_num)
  real    (c_double ) , intent(in)          :: een_rescaled_n_gl(elec_num,4,nucl_num,0:cord_num,walk_num)
  real    (c_double ) , intent(out)         :: factor_een_gl(elec_num,4,walk_num)
  integer(qmckl_exit_code)                  :: info

  integer*8 :: i, a, j, l, k, 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)
        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(i, j, k, nw) * een_rescaled_n(i, a, m, nw)
                 accu2 = accu2 + een_rescaled_e(i, j, k, nw) * een_rescaled_n(i, a, m + l, nw)
                 daccu(1:4) = daccu(1:4) + een_rescaled_e_gl(j, 1:4, i, k, nw) *   &
                      een_rescaled_n(i, a, m, nw)
                 daccu2(1:4) = daccu2(1:4) + een_rescaled_e_gl(j, 1:4, i, k, nw) * &
                      een_rescaled_n(i, a, m + l, nw)
              end do
              factor_een_gl(j, 1:4, nw) = factor_een_gl(j, 1:4, nw) +   &
                   (accu * een_rescaled_n_gl(j, 1:4, a, m + l, nw)      &
                   + daccu(1:4) * een_rescaled_n(j, a, m + l, nw)       &
                   + daccu2(1:4) * een_rescaled_n(j, a, m, nw)          &
                   + accu2 * een_rescaled_n_gl(j, 1:4, a, m, nw)) * cn

              factor_een_gl(j, 4, nw) = factor_een_gl(j, 4, nw) + 2.0d0 * ( &
                   daccu (1) * een_rescaled_n_gl(j, 1, a, m + l, nw) +      &
                   daccu (2) * een_rescaled_n_gl(j, 2, a, m + l, nw) +      &
                   daccu (3) * een_rescaled_n_gl(j, 3, a, m + l, nw) +      &
                   daccu2(1) * een_rescaled_n_gl(j, 1, a, m, nw    ) +      &
                   daccu2(2) * een_rescaled_n_gl(j, 2, a, m, nw    ) +      &
                   daccu2(3) * 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_naive

4.3.7.3. Compute GL

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

function 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)&
     result(info) bind(C)
  use, intrinsic :: iso_c_binding
  use qmckl
  implicit none
  integer(qmckl_context), intent(in), value :: context
  integer (c_int64_t) , intent(in)  , value :: walk_num
  integer (c_int64_t) , intent(in)  , value :: elec_num
  integer (c_int64_t) , intent(in)  , value :: nucl_num
  integer (c_int64_t) , intent(in)  , value :: cord_num
  integer (c_int64_t) , intent(in)  , value :: dim_c_vector
  real    (c_double ) , intent(in)          :: c_vector_full(nucl_num,dim_c_vector)
  integer (c_int64_t) , intent(in)          :: lkpm_combined_index(dim_c_vector,4)
  real    (c_double ) , intent(in)          :: tmp_c(elec_num,nucl_num,0:cord_num,0:cord_num-1,walk_num)
  real    (c_double ) , intent(in)          :: dtmp_c(elec_num,4,nucl_num,0:cord_num,0:cord_num-1,walk_num)
  real    (c_double ) , intent(in)          :: een_rescaled_n(elec_num,nucl_num,0:cord_num,walk_num)
  real    (c_double ) , intent(in)          :: een_rescaled_n_gl(elec_num,4,nucl_num,0:cord_num,walk_num)
  real    (c_double ) , intent(out)         :: factor_een_gl(elec_num,4,walk_num)
  integer(qmckl_exit_code)                  :: info

  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

  if (cord_num == 0) then
     factor_een_gl = 0.0d0
     return
  end if

  do nw =1, walk_num
     factor_een_gl(:,:,nw) = 0.0d0
     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) + &
                      tmp_c (j,   a,m+l,k,nw) * een_rescaled_n_gl(j,ii,a,m  ,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)   &
                      ) * 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

4.3.7.4. Compute Gradient only

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_grad`	`double[walk_num][3][elec_num]`	out	Derivative of Electron-nucleus jastrow

function qmckl_compute_jastrow_champ_factor_een_grad_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_grad) &
     bind(C) result(info)

  use, intrinsic :: iso_c_binding
  use qmckl
  implicit none
  integer(qmckl_context), intent(in), value :: context
  integer (c_int64_t) , intent(in)  , value :: walk_num
  integer (c_int64_t) , intent(in)  , value :: elec_num
  integer (c_int64_t) , intent(in)  , value :: nucl_num
  integer (c_int64_t) , intent(in)  , value :: cord_num
  integer (c_int64_t) , intent(in)  , value :: dim_c_vector
  real    (c_double ) , intent(in)          :: c_vector_full(nucl_num,dim_c_vector)
  integer (c_int64_t) , intent(in)          :: lkpm_combined_index(dim_c_vector,4)
  real    (c_double ) , intent(in)          :: tmp_c(elec_num,nucl_num,0:cord_num,0:cord_num-1,walk_num)
  real    (c_double ) , intent(in)          :: dtmp_c(elec_num,4,nucl_num,0:cord_num,0:cord_num-1,walk_num)
  real    (c_double ) , intent(in)          :: een_rescaled_n(elec_num,nucl_num,0:cord_num,walk_num)
  real    (c_double ) , intent(in)          :: een_rescaled_n_gl(elec_num,4,nucl_num,0:cord_num,walk_num)
  real    (c_double ) , intent(out)         :: factor_een_grad(elec_num,3,walk_num)
  integer(qmckl_exit_code)                  :: info

  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


  if (cord_num == 0) then
     factor_een_grad = 0.0d0
     return
  end if

  do nw =1, walk_num
     factor_een_grad(:,:,nw) = 0.0d0
     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, 3
              do j = 1, elec_num
                 factor_een_grad(j,ii,nw) = factor_een_grad(j,ii,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  ,k,nw)     * een_rescaled_n_gl(j,ii,a,m+l,nw) + &
                      tmp_c(j,a,m+l,k,nw)     * een_rescaled_n_gl(j,ii,a,m  ,nw)   &
                      ) * cn
              end do
           end do

        end do
     end do
  end do

end function qmckl_compute_jastrow_champ_factor_een_grad_doc

Test

4.3.8. Electron-electron-nucleus Jastrow Parameter derivatives

Calculate the derivatives of electron-electron-nuclear three-body jastrow component factor_een wit respect to the parameters.

TODO: write equations.

4.3.8.1. Get

qmckl_exit_code
qmckl_get_jastrow_champ_factor_een_pderiv(qmckl_context context,
                                          double* const factor_een_pderiv,
                                          const int64_t size_max);

Fortran interface

interface
   integer(qmckl_exit_code) function qmckl_get_jastrow_champ_factor_een_pderiv (context, &
        factor_een_pderiv, 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
     real(c_double),   intent(out)               :: factor_een_pderiv(size_max)
   end function qmckl_get_jastrow_champ_factor_een_pderiv
end interface

4.3.8.2. Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`walk_num`	`int64_t`	in	Number of walkers
`elec_num`	`int64_t`	in	Number of electrons
`nucl_num`	`int64_t`	in	Number of nuclei
`type_nucl_num`	`int64_t`	in	Number of unique nuclei
`type_nucl_vector`	`int64_t[nucl_num]`	in	IDs of unique 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_pderiv`	`double[type_nucl_num][dim_c_vector]`	out	Electron-nucleus jastrow

function qmckl_compute_jastrow_champ_factor_een_pderiv_doc( &
     context, walk_num, elec_num, nucl_num, type_nucl_num, &
     type_nucl_vector, cord_num, dim_c_vector, c_vector_full, &
     lkpm_combined_index, tmp_c, een_rescaled_n, factor_een_pderiv) &
     result(info) bind(C)
  use, intrinsic :: iso_c_binding
  use qmckl
  implicit none
  integer(qmckl_context), intent(in), value :: context
  integer (c_int64_t) , intent(in)  , value :: walk_num
  integer (c_int64_t) , intent(in)  , value :: elec_num
  integer (c_int64_t) , intent(in)  , value :: nucl_num
  integer (c_int64_t) , intent(in)  , value :: type_nucl_num
  integer (c_int64_t) , intent(in)          :: type_nucl_vector(nucl_num)
  integer (c_int64_t) , intent(in)  , value :: cord_num
  integer (c_int64_t) , intent(in)  , value :: dim_c_vector
  real    (c_double ) , intent(in)          :: c_vector_full(nucl_num,dim_c_vector)
  integer (c_int64_t) , intent(in)          :: lkpm_combined_index(dim_c_vector,4)
  real    (c_double ) , intent(in)          :: tmp_c(elec_num,nucl_num,0:cord_num,0:cord_num-1,walk_num)
  real    (c_double ) , intent(in)          :: een_rescaled_n(elec_num,nucl_num,0:cord_num,walk_num)
  real    (c_double ) , intent(out)         :: factor_een_pderiv(dim_c_vector, type_nucl_num)
  integer(qmckl_exit_code)                  :: info

  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_pderiv = 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
           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_pderiv(n, type_nucl_vector(a)+1) = factor_een_pderiv(n, type_nucl_vector(a)+1) + accu
        end do
     end do
  end do

end function qmckl_compute_jastrow_champ_factor_een_pderiv_doc

4.3.9. Electron-electron-nucleus Parameter Derivative of the gradient and Laplacian

The derivative of factor_een_gl with respect to all the jastrow parameters C is computed and stored into factor_een_gl_pderiv

TODO: write equations

4.3.9.1. Get

qmckl_exit_code
qmckl_get_jastrow_champ_factor_een_gl_pderiv(qmckl_context context,
                                     double* const factor_een_gl_pderiv,
                                     const int64_t size_max);

Fortran interface

interface
   integer(qmckl_exit_code) function qmckl_get_jastrow_champ_factor_een_gl_pderiv (context, &
        factor_een_gl_pderiv, 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
     real(c_double),   intent(out)               :: factor_een_gl_pderiv(size_max)
   end function qmckl_get_jastrow_champ_factor_een_gl_pderiv
end interface

4.3.9.2. Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`walk_num`	`int64_t`	in	Number of walkers
`elec_num`	`int64_t`	in	Number of electrons
`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	The nucl type index of each nucleus
`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_pderiv`	`double[type_nucl_num][dim_c_vector][elec_num][4]`	out	Derivative of Electron-nucleus jastrow

function qmckl_compute_jastrow_champ_factor_een_gl_pderiv_doc( &
     context, walk_num, elec_num, nucl_num, type_nucl_num, type_nucl_vector, &
     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_pderiv)&
     result(info) bind(C)
  use, intrinsic :: iso_c_binding
  use qmckl
  implicit none
  integer(qmckl_context), intent(in), value :: context
  integer (c_int64_t) , intent(in)  , value :: walk_num
  integer (c_int64_t) , intent(in)  , value :: elec_num
  integer (c_int64_t) , intent(in)  , value :: nucl_num
  integer (c_int64_t) , intent(in)  , value :: type_nucl_num
  integer (c_int64_t) , intent(in)          :: type_nucl_vector(nucl_num)
  integer (c_int64_t) , intent(in)  , value :: cord_num
  integer (c_int64_t) , intent(in)  , value :: dim_c_vector
  real    (c_double ) , intent(in)          :: c_vector_full(nucl_num,dim_c_vector)
  integer (c_int64_t) , intent(in)          :: lkpm_combined_index(dim_c_vector,4)
  real    (c_double ) , intent(in)          :: tmp_c(elec_num,nucl_num,0:cord_num,0:cord_num-1,walk_num)
  real    (c_double ) , intent(in)          :: dtmp_c(elec_num,4,nucl_num,0:cord_num,0:cord_num-1,walk_num)
  real    (c_double ) , intent(in)          :: een_rescaled_n(elec_num,nucl_num,0:cord_num,walk_num)
  real    (c_double ) , intent(in)          :: een_rescaled_n_gl(elec_num,4,nucl_num,0:cord_num,walk_num)
  real    (c_double ) , intent(out)         :: factor_een_gl_pderiv(4,elec_num,dim_c_vector,type_nucl_num)
  integer(qmckl_exit_code)                  :: info

  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_pderiv = 0.0d0
  if (cord_num == 0) then
     return
  end if

  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_pderiv(ii,j,n,type_nucl_vector(a)+1) = factor_een_gl_pderiv(ii,j,n,type_nucl_vector(a)+1) + ( &
                      tmp_c (j,   a,m  ,k,nw) * een_rescaled_n_gl(j,ii,a,m+l,nw) + &
                      tmp_c (j,   a,m+l,k,nw) * een_rescaled_n_gl(j,ii,a,m  ,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)   &
                      )
              end do
           end do

           do j = 1, elec_num
              factor_een_gl_pderiv(4,j,n,type_nucl_vector(a)+1) = factor_een_gl_pderiv(4,j,n,type_nucl_vector(a)+1) +  ( &
                   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)    &
                   ) * 2
           end do
        end do
     end do
  end do

end function qmckl_compute_jastrow_champ_factor_een_gl_pderiv_doc

4.4. Total Jastrow

4.4.1. Value

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

4.4.1.1. 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
     real(c_double),   intent(out)               :: value(size_max)
   end function qmckl_get_jastrow_champ_value
end interface

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

function qmckl_compute_jastrow_champ_value_doc(context, &
     walk_num, f_ee, f_en, f_een, value) &
     result(info) bind(C)
  use, intrinsic :: iso_c_binding
  use qmckl
  implicit none
  integer(qmckl_context), intent(in), value :: context
  integer (c_int64_t) , intent(in)  , value :: walk_num
  real    (c_double ) , intent(in)          :: f_ee(walk_num)
  real    (c_double ) , intent(in)          :: f_en(walk_num)
  real    (c_double ) , intent(in)          :: f_een(walk_num)
  real    (c_double ) , intent(out)         :: value(walk_num)
  integer(qmckl_exit_code)                  :: info

  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

4.4.2. 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] \]

4.4.2.1. Get

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

qmckl_exit_code
qmckl_get_jastrow_champ_grad(qmckl_context context,
                            double* const grad,
                            const int64_t size_max);

qmckl_exit_code
qmckl_get_jastrow_champ_grad(qmckl_context context,
                            double* const grad,
                            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
     real(c_double),   intent(out)               :: gl(size_max)
   end function qmckl_get_jastrow_champ_gl
end interface

4.4.2.2. Compute GL

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

function qmckl_compute_jastrow_champ_gl_doc(context, &
     walk_num, elec_num, value, gl_ee, gl_en, gl_een, gl) &
     bind(C) result(info)
  use qmckl
  use, intrinsic :: iso_c_binding
  implicit none
  integer(qmckl_context), intent(in), value :: context
  integer (c_int64_t) , intent(in)  , value :: walk_num
  integer (c_int64_t) , intent(in)  , value :: elec_num
  real    (c_double ) , intent(in)          :: value(walk_num)
  real    (c_double ) , intent(in)          :: gl_ee(elec_num,4,walk_num)
  real    (c_double ) , intent(in)          :: gl_en(elec_num,4,walk_num)
  real    (c_double ) , intent(in)          :: gl_een(elec_num,4,walk_num)
  real    (c_double ) , intent(out)         :: gl(elec_num,4,walk_num)

  integer(qmckl_exit_code)                 :: info
  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

Test

4.4.2.3. Compute Gradient only

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
`grad_een`	`double[walk_num][3][elec_num]`	in	eeN component
`grad`	`double[walk_num][3][elec_num]`	out	Total Jastrow factor

function qmckl_compute_jastrow_champ_grad_doc(context, &
     walk_num, elec_num, value, gl_ee, gl_en, grad_een, grad) &
     result(info) bind(C)
  use, intrinsic :: iso_c_binding
  use qmckl
  implicit none
  integer(qmckl_context), intent(in), value :: context
  integer (c_int64_t) , intent(in)  , value :: walk_num
  integer (c_int64_t) , intent(in)  , value :: elec_num
  real    (c_double ) , intent(in)          :: value(walk_num)
  real    (c_double ) , intent(in)          :: gl_ee(elec_num,4,walk_num)
  real    (c_double ) , intent(in)          :: gl_en(elec_num,4,walk_num)
  real    (c_double ) , intent(in)          :: grad_een(elec_num,3,walk_num)
  real    (c_double ) , intent(out)         :: grad(elec_num,3,walk_num)

  integer(qmckl_exit_code) :: info
  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,3
        do i = 1, elec_num
           grad(i,j,k) = gl_ee(i,j,k) + gl_en(i,j,k) + grad_een(i,j,k)
        end do
     end do
     grad(:,:,k) = grad(:,:,k) * value(k)
  end do


end function qmckl_compute_jastrow_champ_grad_doc

Test

CHAMP Jastrow Factor

Table of Contents

1. Introduction

1.1. Reformulation of the three-body part

2.

3. Context

3.1. Data structure

3.2. Initialization functions

3.2.0.1. Fortran interface

3.3. Access functions

3.3.0.1. Fortran interface

4. Computation

4.1. Electron-electron component

4.1.1. Asymptotic component

4.1.1.1. Get

4.1.1.2. Compute

4.1.2. Electron-electron rescaled distances

4.1.2.1. Get

4.1.2.2. Compute

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

4.1.3.1. Get

4.1.3.2. Compute

4.1.4. Electron-electron component

4.1.4.1. Get

4.1.4.2. Compute

4.1.5. Derivative

4.1.5.1. Get

4.1.5.2. Compute

4.1.6. Parameter Derivative of the Asymptotic component

4.1.6.1. Get

4.1.6.2. Compute

4.1.7. Parameter Derivative

4.1.7.1. Get

4.1.7.2. Compute

4.1.8. Parameter Derivative of the gradient and Laplacian

4.1.8.1. Get

4.1.8.2. Compute

4.2. Electron-nucleus component

4.2.1. Asymptotic component

4.2.1.1. Get

4.2.1.2. Compute

4.2.2. Electron-nucleus rescaled distances

4.2.2.1. Get

4.2.2.2. Compute

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

4.2.3.1. Get

4.2.3.2. Compute

4.2.4. Electron-nucleus component

4.2.4.1. Get

4.2.4.2. Compute

4.2.5. Derivative

4.2.5.1. Get

4.2.5.2. Compute

4.2.5.3. Test

4.2.6. Parameter Derivative of the Asymptotic component

4.2.6.1. Get

4.2.6.2. Compute

4.2.7. Parameter Derivative

4.2.7.1. Get

4.2.7.2. Compute

4.2.8. Parameter Derivative of the gradient and Laplacian

4.2.8.1. Get

4.2.8.2. Compute

4.3. Electron-electron-nucleus component

4.3.1. Electron-electron rescaled distances in \(J_\text{eeN}\)

4.3.1.1. Get

4.3.1.2. Compute

4.3.1.3. Test

4.3.2. Electron-electron rescaled distances derivatives in \(J_\text{eeN}\)

4.3.2.1. Get

4.3.2.2. Compute

4.3.3. Electron-nucleus rescaled distances in \(J_\text{eeN}\)

4.3.3.1. Get

4.3.3.2. Compute

4.3.4. Electron-nucleus rescaled distances derivatives in \(J_\text{eeN}\)

4.3.4.1. Get

4.3.4.2. Compute

4.3.5. Temporary arrays for electron-electron-nucleus Jastrow \(f_{een}\)

4.3.5.1. Compute dimcvector

4.3.5.2. Get

4.3.5.1. Compute dim_c_vector

4.3.5.3. Compute c_vector_full

4.3.5.4. Compute lkpm_combined_index

4.3.5.5. Compute tmp_c

4.3.5.6. Compute dtmp_c