CHAMP Jastrow Factor
Table of Contents
- 1. Introduction
- 2. Context
- 3. Computation
- 3.1. Electron-electron component
- 3.2. Electron-nucleus component
- 3.3. Electron-electron-nucleus component
- 3.3.1. Electron-electron rescaled distances in \(J_\text{eeN}\)
- 3.3.2. Electron-electron rescaled distances derivatives in \(J_\text{eeN}\)
- 3.3.3. Electron-nucleus rescaled distances in \(J_\text{eeN}\)
- 3.3.4. Electron-nucleus rescaled distances derivatives in \(J_\text{eeN}\)
- 3.3.5. Temporary arrays for electron-electron-nucleus Jastrow \(f_{een}\)
- 3.3.6. Electron-electron-nucleus Jastrow \(f_{een}\)
- 3.3.7. Electron-electron-nucleus Jastrow \(f_{een}\) derivative
- 3.4. Total Jastrow
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.
2 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 |
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 |
Number of unique C coefficients |
asymp_jasa |
double[type_nucl_num] |
Asymptotic component |
asymp_jasa_date |
uint64_t |
Ladt modification of the asymptotic component |
asymp_jasb |
double[2] |
Asymptotic component (up- or down-spin) |
asymp_jasb_date |
uint64_t |
Ladt modification of the asymptotic component |
c_vector_full |
double[dim_c_vector][nucl_num] |
vector of non-zero coefficients |
c_vector_full_date |
uint64_t |
Keep track of changes here |
lkpm_combined_index |
int64_t[4][dim_c_vector] |
Transform l,k,p, and m into consecutive indices |
lkpm_combined_index_date |
uint64_t |
Transform l,k,p, and m into consecutive indices |
tmp_c |
double[walk_num][cord_num][cord_num+1][nucl_num][elec_num] |
vector of non-zero coefficients |
dtmp_c |
double[walk_num][elec_num][4][nucl_num][cord_num+1][cord_num] |
vector of non-zero coefficients |
ee_distance_rescaled |
double[walk_num][num][num] |
Electron-electron rescaled distances |
ee_distance_rescaled_date |
uint64_t |
Last modification date of the electron-electron distances |
ee_distance_rescaled_gl |
double[walk_num][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_en |
double[walk_num] |
Jastrow factor: electron-nucleus part |
factor_en_date |
uint64_t |
Jastrow factor: electron-nucleus part |
factor_een |
double[walk_num] |
Jastrow factor: electron-electron-nucleus part |
factor_een_date |
uint64_t |
Jastrow factor: electron-electron-nucleus part |
factor_ee_gl |
double[walk_num][4][elec_num] |
Derivative of the Jastrow factor: electron-electron-nucleus part |
factor_ee_gl_date |
uint64_t |
Keep track of the date for the derivative |
factor_en_gl |
double[walk_num][4][elec_num] |
Derivative of the Jastrow factor: electron-electron-nucleus part |
factor_en_gl_date |
uint64_t |
Keep track of the date for the en derivative |
factor_een_gl |
double[walk_num][4][elec_num] |
Derivative of the Jastrow factor: electron-electron-nucleus part |
factor_een_gl_date |
uint64_t |
Keep track of the date for the een derivative |
value |
double[walk_num] |
Value of the Jastrow factor |
value_date |
uint64_t |
Keep track of the date |
gl |
double[walk_num][4][elec_num] |
Gradient and Laplacian of the Jastrow factor |
value_date |
uint64_t |
Keep track of the date |
2.1 Data structure
typedef struct qmckl_jastrow_champ_struct{ int64_t * restrict lkpm_combined_index; int64_t * restrict type_nucl_vector; double * restrict asymp_jasa; double asymp_jasb[2]; double * restrict a_vector; double * restrict b_vector; double * restrict c_vector; double * restrict c_vector_full; double * restrict dtmp_c; double * restrict ee_distance_rescaled; double * restrict ee_distance_rescaled_gl; double * restrict een_rescaled_e; double * restrict een_rescaled_e_gl; double * restrict een_rescaled_n; double * restrict een_rescaled_n_gl; double * restrict en_distance_rescaled; double * restrict en_distance_rescaled_gl; double * restrict factor_ee; double * restrict factor_ee_gl; double * restrict factor_een; double * restrict factor_een_gl; double * restrict factor_en; double * restrict factor_en_gl; double * restrict rescale_factor_en; double * restrict tmp_c; double * restrict value; double * restrict gl; int64_t aord_num; int64_t bord_num; int64_t cord_num; int64_t dim_c_vector; int64_t type_nucl_num; uint64_t asymp_jasa_date; uint64_t asymp_jasb_date; uint64_t c_vector_full_date; uint64_t dim_c_vector_date; uint64_t dtmp_c_date; uint64_t ee_distance_rescaled_date; uint64_t ee_distance_rescaled_gl_date; uint64_t een_rescaled_e_date; uint64_t een_rescaled_e_gl_date; uint64_t een_rescaled_n_date; uint64_t een_rescaled_n_gl_date; uint64_t en_distance_rescaled_date; uint64_t en_distance_rescaled_gl_date; uint64_t factor_ee_date; uint64_t factor_ee_gl_date; uint64_t factor_een_date; uint64_t factor_een_gl_date; uint64_t factor_en_date; uint64_t factor_en_gl_date; uint64_t lkpm_combined_index_date; uint64_t tmp_c_date; uint64_t value_date; uint64_t gl_date; double rescale_factor_ee; int32_t uninitialized; int32_t spin_independent; bool provided; } qmckl_jastrow_champ_struct;
The uninitialized integer contains one bit set to one for each
initialization function which has not been called. It becomes equal
to zero after all initialization functions have been called. The
struct is then initialized and provided == true.
Some values are initialized by default, and are not concerned by
this mechanism.
qmckl_exit_code qmckl_init_jastrow_champ(qmckl_context context);
qmckl_exit_code qmckl_init_jastrow_champ(qmckl_context context) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return false; } qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); ctx->jastrow_champ.uninitialized = (1 << 11) - 1; /* Default values */ ctx->jastrow_champ.aord_num = -1; ctx->jastrow_champ.bord_num = -1; ctx->jastrow_champ.cord_num = -1; ctx->jastrow_champ.dim_c_vector = -1; ctx->jastrow_champ.type_nucl_num = -1; ctx->jastrow_champ.spin_independent = -1; return QMCKL_SUCCESS; }
2.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.
2.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
2.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);
2.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
2.4 Test
/* Reference input data */ int64_t walk_num = n2_walk_num; int64_t elec_num = n2_elec_num; int64_t elec_up_num = n2_elec_up_num; int64_t elec_dn_num = n2_elec_dn_num; int64_t nucl_num = n2_nucl_num; double rescale_factor_ee = 0.6; double rescale_factor_en[2] = { 0.6, 0.6 }; double* elec_coord = &(n2_elec_coord[0][0][0]); const double* nucl_charge = n2_charge; double* nucl_coord = &(n2_nucl_coord[0][0]); int64_t size_max; /* Provide Electron data */ qmckl_exit_code rc; assert(!qmckl_electron_provided(context)); rc = qmckl_check(context, qmckl_set_electron_num (context, elec_up_num, elec_dn_num) ); assert(rc == QMCKL_SUCCESS); assert(qmckl_electron_provided(context)); rc = qmckl_check(context, qmckl_set_electron_coord (context, 'N', walk_num, elec_coord, walk_num*3*elec_num) ); assert(rc == QMCKL_SUCCESS); double elec_coord2[walk_num*3*elec_num]; rc = qmckl_check(context, qmckl_get_electron_coord (context, 'N', elec_coord2, walk_num*3*elec_num) ); assert(rc == QMCKL_SUCCESS); for (int64_t i=0 ; i<3*elec_num ; ++i) { assert( elec_coord[i] == elec_coord2[i] ); } /* Provide Nucleus data */ assert(!qmckl_nucleus_provided(context)); rc = qmckl_check(context, qmckl_set_nucleus_num (context, nucl_num) ); assert(rc == QMCKL_SUCCESS); assert(!qmckl_nucleus_provided(context)); double nucl_coord2[3*nucl_num]; rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2, 3*nucl_num); assert(rc == QMCKL_NOT_PROVIDED); rc = qmckl_check(context, qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), 3*nucl_num) ); assert(rc == QMCKL_SUCCESS); assert(!qmckl_nucleus_provided(context)); rc = qmckl_check(context, qmckl_get_nucleus_coord (context, 'N', nucl_coord2, nucl_num*3) ); assert(rc == QMCKL_SUCCESS); for (int64_t k=0 ; k<3 ; ++k) { for (int64_t i=0 ; i<nucl_num ; ++i) { assert( nucl_coord[nucl_num*k+i] == nucl_coord2[3*i+k] ); } } rc = qmckl_check(context, qmckl_get_nucleus_coord (context, 'T', nucl_coord2, nucl_num*3) ); assert(rc == QMCKL_SUCCESS); for (int64_t i=0 ; i<3*nucl_num ; ++i) { assert( nucl_coord[i] == nucl_coord2[i] ); } double nucl_charge2[nucl_num]; rc = qmckl_get_nucleus_charge(context, nucl_charge2, nucl_num); assert(rc == QMCKL_NOT_PROVIDED); rc = qmckl_check(context, qmckl_set_nucleus_charge(context, nucl_charge, nucl_num) ); assert(rc == QMCKL_SUCCESS); rc = qmckl_check(context, qmckl_get_nucleus_charge(context, nucl_charge2, nucl_num) ); assert(rc == QMCKL_SUCCESS); for (int64_t i=0 ; i<nucl_num ; ++i) { assert( nucl_charge[i] == nucl_charge2[i] ); } assert(qmckl_nucleus_provided(context));
3 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.
3.1 Electron-electron component
3.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} \]
3.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
3.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, intrinsic :: iso_c_binding 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
qmckl_exit_code qmckl_compute_jastrow_champ_asymp_jasb_hpc (const qmckl_context context, const int64_t bord_num, const double* b_vector, const double rescale_factor_ee, const int32_t spin_independent, double* const asymp_jasb ) { if (context == QMCKL_NULL_CONTEXT) { return QMCKL_INVALID_CONTEXT; } if (bord_num < 0) { return QMCKL_INVALID_ARG_2; } const double kappa_inv = 1.0 / rescale_factor_ee; const double asym_one = b_vector[0] * kappa_inv / (1.0 + b_vector[1] * kappa_inv); double f = 0.; double x = kappa_inv; for (int k = 2; k <= bord_num; ++k) { x *= kappa_inv; f = f + b_vector[k]*x; } asymp_jasb[0] = spin_independent == 1 ? asym_one + f : 0.5 * asym_one + f; asymp_jasb[1] = asym_one + f; return QMCKL_SUCCESS; }
qmckl_exit_code qmckl_compute_jastrow_champ_asymp_jasb (const qmckl_context context, const int64_t bord_num, const double* b_vector, const double rescale_factor_ee, const int32_t spin_independent, double* const asymp_jasb ) { #ifdef HAVE_HPC return qmckl_compute_jastrow_champ_asymp_jasb_hpc #else return qmckl_compute_jastrow_champ_asymp_jasb_doc #endif (context, bord_num, b_vector, rescale_factor_ee, spin_independent, asymp_jasb); }
3.1.1.3 Test
assert(qmckl_electron_provided(context)); int64_t type_nucl_num = n2_type_nucl_num; int64_t* type_nucl_vector = &(n2_type_nucl_vector[0]); int64_t aord_num = n2_aord_num; int64_t bord_num = n2_bord_num; int64_t cord_num = n2_cord_num; double* a_vector = &(n2_a_vector[0][0]); double* b_vector = &(n2_b_vector[0]); double* c_vector = &(n2_c_vector[0][0]); int64_t dim_c_vector=0; assert(!qmckl_jastrow_champ_provided(context)); /* Set the data */ rc = qmckl_check(context, qmckl_set_jastrow_champ_spin_independent(context, 0) ); rc = qmckl_check(context, qmckl_set_jastrow_champ_aord_num(context, aord_num) ); rc = qmckl_check(context, qmckl_set_jastrow_champ_bord_num(context, bord_num) ); rc = qmckl_check(context, qmckl_set_jastrow_champ_cord_num(context, cord_num) ); assert(rc == QMCKL_SUCCESS); rc = qmckl_check(context, qmckl_set_jastrow_champ_type_nucl_num(context, type_nucl_num) ); assert(rc == QMCKL_SUCCESS); rc = qmckl_check(context, qmckl_set_jastrow_champ_type_nucl_vector(context, type_nucl_vector, nucl_num) ); assert(rc == QMCKL_SUCCESS); rc = qmckl_check(context, qmckl_set_jastrow_champ_a_vector(context, a_vector,(aord_num+1)*type_nucl_num) ); assert(rc == QMCKL_SUCCESS); rc = qmckl_check(context, qmckl_set_jastrow_champ_b_vector(context, b_vector,(bord_num+1)) ); assert(rc == QMCKL_SUCCESS); rc = qmckl_check(context, qmckl_get_jastrow_champ_dim_c_vector(context, &dim_c_vector) ); assert(rc == QMCKL_SUCCESS); rc = qmckl_check(context, qmckl_set_jastrow_champ_c_vector(context, c_vector, dim_c_vector*type_nucl_num) ); assert(rc == QMCKL_SUCCESS); double k_ee = 0.; double k_en[2] = { 0., 0. }; rc = qmckl_check(context, qmckl_set_jastrow_champ_rescale_factor_en(context, rescale_factor_en, type_nucl_num) ); assert(rc == QMCKL_SUCCESS); rc = qmckl_check(context, qmckl_set_jastrow_champ_rescale_factor_ee(context, rescale_factor_ee) ); assert(rc == QMCKL_SUCCESS); rc = qmckl_check(context, qmckl_get_jastrow_champ_rescale_factor_ee (context, &k_ee) ); assert(rc == QMCKL_SUCCESS); assert(k_ee == rescale_factor_ee); rc = qmckl_check(context, qmckl_get_jastrow_champ_rescale_factor_en (context, &(k_en[0]), type_nucl_num) ); assert(rc == QMCKL_SUCCESS); for (int i=0 ; i<type_nucl_num ; ++i) { assert(k_en[i] == rescale_factor_en[i]); } /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); double asymp_jasb[2]; rc = qmckl_check(context, qmckl_get_jastrow_champ_asymp_jasb(context, asymp_jasb,2) ); // calculate asymp_jasb assert(fabs(asymp_jasb[0]-0.7115733522582638) < 1.e-12); assert(fabs(asymp_jasb[1]-1.043287918508297 ) < 1.e-12);
3.1.2 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]
\]
3.1.2.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
3.1.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 |
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, 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(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
qmckl_exit_code qmckl_compute_jastrow_champ_factor_ee_hpc (const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t up_num, const int64_t bord_num, const double* b_vector, const double* ee_distance_rescaled, const double* asymp_jasb, const int32_t spin_independent, double* const factor_ee ) { if (context == QMCKL_NULL_CONTEXT) { return QMCKL_INVALID_CONTEXT; } if (walk_num <= 0) { return QMCKL_INVALID_ARG_2; } if (elec_num <= 0) { return QMCKL_INVALID_ARG_3; } if (bord_num < 0) { return QMCKL_INVALID_ARG_4; } const int64_t dn_num = elec_num - up_num; const double fshift = 0.5 * (double) ((dn_num-1)*dn_num + (up_num-1)*up_num) * asymp_jasb[0] + (float) (up_num*dn_num) * asymp_jasb[1]; #ifdef HAVE_OPENMP #pragma omp parallel #endif for (int nw = 0; nw < walk_num; ++nw) { double result = 0.; size_t ishift = nw * elec_num * elec_num; if (spin_independent == 1) { for (int j = 0; j < elec_num; ++j ) { const double* xj = &(ee_distance_rescaled[j * elec_num + ishift]); for (int i = 0; i < j ; ++i) { result = result + b_vector[0]*xj[i] / (1. + b_vector[1]*xj[i]); } } } else { for (int j = 0; j < up_num; ++j ) { const double* xj = &(ee_distance_rescaled[j * elec_num + ishift]); for (int i = 0; i < j ; ++i) { result = result + 0.5 * b_vector[0]*xj[i] / (1. + b_vector[1]*xj[i]); } } for (int j = up_num ; j < elec_num; ++j ) { const double* xj = &(ee_distance_rescaled[j * elec_num + ishift]); for (int i = 0; i < up_num; ++i) { result = result + b_vector[0]*xj[i] / (1. + b_vector[1]*xj[i]); } for (int i = up_num ; i < j ; ++i) { result = result + 0.5 * b_vector[0]*xj[i] / (1. + b_vector[1]*xj[i]); } } } result = result - fshift; for (int j=0; j < elec_num; ++j ) { const double* xj = &(ee_distance_rescaled[j * elec_num + ishift]); for (int i=0; i < j ; ++i) { const double x = xj[i]; double xk = x; for (int k = 2; k <= bord_num; ++k) { xk *= x; result = result + b_vector[k] * xk; } } } factor_ee[nw] = result; } return QMCKL_SUCCESS; }
qmckl_exit_code qmckl_compute_jastrow_champ_factor_ee (const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t up_num, const int64_t bord_num, const double* b_vector, const double* ee_distance_rescaled, const double* asymp_jasb, const int32_t spin_independent, double* const factor_ee ) { #ifdef HAVE_HPC return qmckl_compute_jastrow_champ_factor_ee_hpc #else return qmckl_compute_jastrow_champ_factor_ee_doc #endif (context, walk_num, elec_num, up_num, bord_num, b_vector, ee_distance_rescaled, asymp_jasb, spin_independent, factor_ee); }
3.1.2.3 Test
/* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); double factor_ee[walk_num]; rc = qmckl_check(context, qmckl_get_jastrow_champ_factor_ee(context, factor_ee, walk_num) ); // calculate factor_ee printf("%e\n%e\n\n",factor_ee[0],-16.83886184243964); assert(fabs(factor_ee[0]+16.83886184243964) < 1.e-12);
3.1.3 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] \]
3.1.3.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
3.1.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 |
up_num |
int64_t |
in | Number of alpha electrons |
bord_num |
int64_t |
in | Number of coefficients |
b_vector |
double[bord_num+1] |
in | List of coefficients |
ee_distance_rescaled |
double[walk_num][elec_num][elec_num] |
in | Electron-electron distances |
ee_distance_rescaled_gl |
double[walk_num][4][elec_num][elec_num] |
in | Electron-electron distances |
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, 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) :: 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
qmckl_exit_code qmckl_compute_jastrow_champ_factor_ee_gl_hpc(const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t up_num, const int64_t bord_num, const double* b_vector, const double* ee_distance_rescaled, const double* ee_distance_rescaled_gl, const int32_t spin_independent, double* const factor_ee_gl ) { if (context == QMCKL_NULL_CONTEXT) return QMCKL_INVALID_CONTEXT; if (walk_num <= 0) return QMCKL_INVALID_ARG_2; if (elec_num <= 0) return QMCKL_INVALID_ARG_3; if (up_num <= 0) return QMCKL_INVALID_ARG_4; if (bord_num < 0) return QMCKL_INVALID_ARG_5; if (b_vector == NULL) return QMCKL_INVALID_ARG_6; if (ee_distance_rescaled == NULL) return QMCKL_INVALID_ARG_7; if (ee_distance_rescaled_gl == NULL) return QMCKL_INVALID_ARG_8; if (spin_independent & (int32_t) (-2)) return QMCKL_INVALID_ARG_8; if (factor_ee_gl == NULL) return QMCKL_INVALID_ARG_9; double kf[bord_num+1]; for (int k=0 ; k<=bord_num ; ++k) { kf[k] = (double) k; } #ifdef HAVE_OPENMP #pragma omp parallel for #endif for (int nw = 0; nw < walk_num; ++nw) { double xk[bord_num+1]; bool touched = false; for (int j = 0; j < elec_num; ++j) { const double* dxj = &ee_distance_rescaled_gl[4*elec_num*(j+nw*elec_num)]; const double* xj = &ee_distance_rescaled [ elec_num*(j+nw*elec_num)]; double * restrict factor_ee_gl_0 = &(factor_ee_gl[nw*elec_num*4]); double * restrict factor_ee_gl_1 = factor_ee_gl_0 + elec_num; double * restrict factor_ee_gl_2 = factor_ee_gl_1 + elec_num; double * restrict factor_ee_gl_3 = factor_ee_gl_2 + elec_num; for (int i = 0; i < elec_num; ++i) { if (j == i) continue; double x = xj[i]; const double denom = 1.0 + b_vector[1]*x; const double invdenom = 1.0 / denom; const double invdenom2 = invdenom * invdenom; const double* restrict dx = dxj + 4*i; const double grad_c2 = dx[0]*dx[0] + dx[1]*dx[1] + dx[2]*dx[2]; double f = b_vector[0] * invdenom2; if ((spin_independent == 0) && ( ((i < up_num) && (j < up_num)) || ((i >= up_num) && (j >= up_num))) ) { f *= 0.5; } if (touched) { factor_ee_gl_0[i] = factor_ee_gl_0[i] + f*dx[0]; factor_ee_gl_1[i] = factor_ee_gl_1[i] + f*dx[1]; factor_ee_gl_2[i] = factor_ee_gl_2[i] + f*dx[2]; factor_ee_gl_3[i] = factor_ee_gl_3[i] + f*dx[3]; } else { touched = true; factor_ee_gl_0[i] = f*dx[0]; factor_ee_gl_1[i] = f*dx[1]; factor_ee_gl_2[i] = f*dx[2]; factor_ee_gl_3[i] = f*dx[3]; } factor_ee_gl_3[i] = factor_ee_gl_3[i] - f*grad_c2*invdenom*2.0 * b_vector[1]; xk[0] = 1.0; for (int k=1 ; k<= bord_num ; ++k) { xk[k] = xk[k-1]*x; } for (int k=2 ; k<= bord_num ; ++k) { const double f1 = b_vector[k] * kf[k] * xk[k-2]; const double f2 = f1*xk[1]; factor_ee_gl_0[i] = factor_ee_gl_0[i] + f2*dx[0]; factor_ee_gl_1[i] = factor_ee_gl_1[i] + f2*dx[1]; factor_ee_gl_2[i] = factor_ee_gl_2[i] + f2*dx[2]; factor_ee_gl_3[i] = factor_ee_gl_3[i] + f2*dx[3]; factor_ee_gl_3[i] = factor_ee_gl_3[i] + f1*kf[k-1]*grad_c2; } } } if (!touched) { memset(&(factor_ee_gl[nw*4*elec_num]), 0, elec_num*4*sizeof(double)); } } return QMCKL_SUCCESS; }
qmckl_exit_code qmckl_compute_jastrow_champ_factor_ee_gl_hpc (const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t up_num, const int64_t bord_num, const double* b_vector, const double* ee_distance_rescaled, const double* ee_distance_rescaled_gl, const int32_t spin_independent, double* const factor_ee_gl );
qmckl_exit_code qmckl_compute_jastrow_champ_factor_ee_gl_doc (const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t up_num, const int64_t bord_num, const double* b_vector, const double* ee_distance_rescaled, const double* ee_distance_rescaled_gl, const int32_t spin_independent, double* const factor_ee_gl );
qmckl_exit_code qmckl_compute_jastrow_champ_factor_ee_gl (const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t up_num, const int64_t bord_num, const double* b_vector, const double* ee_distance_rescaled, const double* ee_distance_rescaled_gl, const int32_t spin_independent, double* const factor_ee_gl ) { #ifdef HAVE_HPC return qmckl_compute_jastrow_champ_factor_ee_gl_hpc #else return qmckl_compute_jastrow_champ_factor_ee_gl_doc #endif (context, walk_num, elec_num, up_num, bord_num, b_vector, ee_distance_rescaled, ee_distance_rescaled_gl, spin_independent, factor_ee_gl ); }
3.1.3.3 Test
/* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); // calculate factor_ee_gl double factor_ee_gl[walk_num][4][elec_num]; rc = qmckl_get_jastrow_champ_factor_ee_gl(context, &(factor_ee_gl[0][0][0]),walk_num*4*elec_num); // check factor_ee_gl printf("%f %f\n", factor_ee_gl[0][0][0], -0.39319353942687446); assert(fabs(factor_ee_gl[0][0][0]+0.39319353942687446) < 1.e-12); printf("%f %f\n", factor_ee_gl[0][1][0], 1.0535615450668214); assert(fabs(factor_ee_gl[0][1][0]-1.0535615450668214) < 1.e-12); printf("%f %f\n", factor_ee_gl[0][2][0],-0.39098406960784515); assert(fabs(factor_ee_gl[0][2][0]+0.39098406960784515) < 1.e-12); printf("%f %f\n", factor_ee_gl[0][3][0],2.8650469630854483); assert(fabs(factor_ee_gl[0][3][0]-2.8650469630854483) < 1.e-12);
3.1.4 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.
3.1.4.1 Get
qmckl_exit_code qmckl_get_jastrow_champ_ee_distance_rescaled(qmckl_context context, double* const distance_rescaled);
3.1.4.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, intrinsic :: iso_c_binding 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
3.1.4.3 Test
assert(qmckl_electron_provided(context)); double ee_distance_rescaled[walk_num * elec_num * elec_num]; rc = qmckl_get_jastrow_champ_ee_distance_rescaled(context, ee_distance_rescaled); // (e1,e2,w) // (0,0,0) == 0. assert(ee_distance_rescaled[0] == 0.); // (1,0,0) == (0,1,0) assert(ee_distance_rescaled[1] == ee_distance_rescaled[elec_num]); // value of (1,0,0) assert(fabs(ee_distance_rescaled[1]-0.6347507420688708) < 1.e-12); // (0,0,1) == 0. assert(ee_distance_rescaled[5*elec_num + 5] == 0.); // (1,0,1) == (0,1,1) assert(ee_distance_rescaled[5*elec_num+6] == ee_distance_rescaled[6*elec_num+5]); // value of (1,0,1) assert(fabs(ee_distance_rescaled[5*elec_num+6]-0.3941735387855409) < 1.e-12);
3.1.5 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.
3.1.5.1 Get
qmckl_exit_code qmckl_get_jastrow_champ_ee_distance_rescaled_gl(qmckl_context context, double* const distance_rescaled_gl);
3.1.5.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, intrinsic :: iso_c_binding 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
3.1.5.3 Test
assert(qmckl_electron_provided(context)); double ee_distance_rescaled_gl[4 * walk_num * elec_num * elec_num]; rc = qmckl_get_jastrow_champ_ee_distance_rescaled_gl(context, ee_distance_rescaled_gl); // TODO: Get exact values //// (e1,e2,w) //// (0,0,0) == 0. //assert(ee_distance[0] == 0.); // //// (1,0,0) == (0,1,0) //assert(ee_distance[1] == ee_distance[elec_num]); // //// value of (1,0,0) //assert(fabs(ee_distance[1]-7.152322512964209) < 1.e-12); // //// (0,0,1) == 0. //assert(ee_distance[elec_num*elec_num] == 0.); // //// (1,0,1) == (0,1,1) //assert(ee_distance[elec_num*elec_num+1] == ee_distance[elec_num*elec_num+elec_num]); // //// value of (1,0,1) //assert(fabs(ee_distance[elec_num*elec_num+1]-6.5517646321055665) < 1.e-12);
3.2 Electron-nucleus component
3.2.1 Asymptotic component for
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}} \]
3.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
3.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 |
integer function qmckl_compute_jastrow_champ_asymp_jasa_f(context, aord_num, type_nucl_num, a_vector, & rescale_factor_en, asymp_jasa) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: aord_num integer*8 , intent(in) :: type_nucl_num double precision , intent(in) :: a_vector(aord_num + 1, type_nucl_num) double precision , intent(in) :: rescale_factor_en(type_nucl_num) double precision , intent(out) :: asymp_jasa(type_nucl_num) integer*8 :: i, j, p double precision :: kappa_inv, x, asym_one info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif if (aord_num < 0) then info = QMCKL_INVALID_ARG_2 return endif do i=1,type_nucl_num kappa_inv = 1.0d0 / rescale_factor_en(i) asymp_jasa(i) = a_vector(1,i) * kappa_inv / (1.0d0 + a_vector(2,i) * kappa_inv) x = kappa_inv do p = 2, aord_num x = x * kappa_inv asymp_jasa(i) = asymp_jasa(i) + a_vector(p+1, i) * x end do end do end function qmckl_compute_jastrow_champ_asymp_jasa_f
qmckl_exit_code qmckl_compute_jastrow_champ_asymp_jasa ( const qmckl_context context, const int64_t aord_num, const int64_t type_nucl_num, const double* a_vector, const double* rescale_factor_en, double* const asymp_jasa );
3.2.1.3 Test
double asympjasa[2]; rc = qmcklgetjastrowchampasympjasa(context, asympjasa, typenuclnum);
// calculate asympjasb printf("%e %e\n", asympjasa[0], -1.75529774); assert(fabs(-1.75529774 - asympjasa[0]) < 1.e-8);
#+endsrc
3.2.2 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
3.2.2.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
3.2.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 |
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, 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 :: 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
qmckl_exit_code qmckl_compute_jastrow_champ_factor_en_doc ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t type_nucl_num, const int64_t* type_nucl_vector, const int64_t aord_num, const double* a_vector, const double* en_distance_rescaled, const double* asymp_jasa, double* const factor_en );
qmckl_exit_code qmckl_compute_jastrow_champ_factor_en ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t type_nucl_num, const int64_t* type_nucl_vector, const int64_t aord_num, const double* a_vector, const double* en_distance_rescaled, const double* asymp_jasa, double* const factor_en ) { #ifdef HAVE_HPC return qmckl_compute_jastrow_champ_factor_en_hpc #else return qmckl_compute_jastrow_champ_factor_en_doc #endif (context, walk_num, elec_num, nucl_num, type_nucl_num, type_nucl_vector, aord_num, a_vector, en_distance_rescaled, asymp_jasa, factor_en ); }
3.2.2.3 Test
/* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); double factor_en[walk_num]; rc = qmckl_get_jastrow_champ_factor_en(context, factor_en,walk_num); // calculate factor_en printf("%f %f\n", factor_en[0], 22.781375792083587); assert(fabs(22.781375792083587 - factor_en[0]) < 1.e-12);
3.2.3 Derivative
Calculate the electron-electron jastrow component factor_en_gl derivative
with respect to the electron coordinates using the en_distance_rescaled and
en_distance_rescaled_gl which are already calculated previously.
TODO: write equations.
3.2.3.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
3.2.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 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, 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 :: 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
3.2.3.3 Test
/* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); // calculate factor_en_gl double factor_en_gl[walk_num][4][elec_num]; rc = qmckl_get_jastrow_champ_factor_en_gl(context, &(factor_en_gl[0][0][0]),walk_num*4*elec_num); // check factor_en_gl assert(fabs( 0.19656663796630847 - factor_en_gl[0][0][0]) < 1.e-12); assert(fabs( -0.3945140890522283 - factor_en_gl[0][1][0]) < 1.e-12); assert(fabs( 0.5082964671286118 - factor_en_gl[0][2][0]) < 1.e-12); assert(fabs( -1.8409460670666289 - factor_en_gl[0][3][0]) < 1.e-12);
3.2.4 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.
3.2.4.1 Get
qmckl_exit_code qmckl_get_electron_en_distance_rescaled(qmckl_context context, double* distance_rescaled);
3.2.4.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][elec_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, intrinsic :: iso_c_binding 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
3.2.4.3 Test
assert(qmckl_electron_provided(context)); assert(qmckl_nucleus_provided(context)); double en_distance_rescaled[walk_num][nucl_num][elec_num]; rc = qmckl_check(context, qmckl_get_electron_en_distance_rescaled(context, &(en_distance_rescaled[0][0][0])) ); assert (rc == QMCKL_SUCCESS); // (e,n,w) in Fortran notation // (1,1,1) assert(fabs(en_distance_rescaled[0][0][0] - 0.4942158656729477) < 1.e-12); // (1,2,1) assert(fabs(en_distance_rescaled[0][1][0] - 1.2464137498005765) < 1.e-12); // (2,1,1) assert(fabs(en_distance_rescaled[0][0][1] - 0.5248654474756858) < 1.e-12); // (1,1,2) assert(fabs(en_distance_rescaled[0][0][5] - 0.19529459944794733) < 1.e-12); // (1,2,2) assert(fabs(en_distance_rescaled[0][1][5] - 1.2091967687767369) < 1.e-12); // (2,1,2) assert(fabs(en_distance_rescaled[0][0][6] - 0.4726452953409436) < 1.e-12);
3.2.5 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.
3.2.5.1 Get
qmckl_exit_code qmckl_get_electron_en_distance_rescaled_gl(qmckl_context context, double* distance_rescaled_gl);
3.2.5.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][elec_num] |
in | Nuclear coordinates |
en_distance_rescaled_gl |
double[walk_num][nucl_num][elec_num][4] |
out | Electron-nucleus distance derivatives |
integer function qmckl_compute_en_distance_rescaled_gl_doc_f(context, elec_num, nucl_num, & type_nucl_num, type_nucl_vector, rescale_factor_en, walk_num, elec_coord, & nucl_coord, en_distance_rescaled_gl) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: elec_num integer*8 , intent(in) :: nucl_num integer*8 , intent(in) :: type_nucl_num integer*8 , intent(in) :: type_nucl_vector(nucl_num) double precision , intent(in) :: rescale_factor_en(nucl_num) integer*8 , intent(in) :: walk_num double precision , intent(in) :: elec_coord(elec_num,walk_num,3) double precision , intent(in) :: nucl_coord(nucl_num,3) double precision , intent(out) :: en_distance_rescaled_gl(4,elec_num,nucl_num,walk_num) integer*8 :: i, k double precision :: coord(3) info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif if (elec_num <= 0) then info = QMCKL_INVALID_ARG_2 return endif if (nucl_num <= 0) then info = QMCKL_INVALID_ARG_3 return endif if (walk_num <= 0) then info = QMCKL_INVALID_ARG_5 return endif do i=1, nucl_num coord(1:3) = nucl_coord(i,1:3) do k=1,walk_num info = qmckl_distance_rescaled_gl(context, 'T', 'T', elec_num, 1_8, & elec_coord(1,k,1), elec_num*walk_num, coord, 1_8, & en_distance_rescaled_gl(1,1,i,k), elec_num, rescale_factor_en(type_nucl_vector(i)+1)) if (info /= QMCKL_SUCCESS) then return endif end do end do end function qmckl_compute_en_distance_rescaled_gl_doc_f
3.2.5.3 Test
assert(qmckl_electron_provided(context)); assert(qmckl_nucleus_provided(context)); double en_distance_rescaled_gl[walk_num][4][nucl_num][elec_num]; rc = qmckl_check(context, qmckl_get_electron_en_distance_rescaled_gl(context, &(en_distance_rescaled_gl[0][0][0][0])) ); assert (rc == QMCKL_SUCCESS); // TODO: check exact values //// (e,n,w) in Fortran notation //// (1,1,1) //assert(fabs(en_distance_rescaled[0][0][0] - 7.546738741619978) < 1.e-12); // //// (1,2,1) //assert(fabs(en_distance_rescaled[0][1][0] - 8.77102435246984) < 1.e-12); // //// (2,1,1) //assert(fabs(en_distance_rescaled[0][0][1] - 3.698922010513608) < 1.e-12); // //// (1,1,2) //assert(fabs(en_distance_rescaled[1][0][0] - 5.824059436060509) < 1.e-12); // //// (1,2,2) //assert(fabs(en_distance_rescaled[1][1][0] - 7.080482110317645) < 1.e-12); // //// (2,1,2) //assert(fabs(en_distance_rescaled[1][0][1] - 3.1804527583077356) < 1.e-12);
3.3 Electron-electron-nucleus component
3.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:
\[ C_{ij,p} = \left[ \exp\left(-\kappa_\text{e}\, r_{ij}\right) \right]^p \]
where \(r_{ij}\) is the matrix of electron-electron distances.
3.3.1.1 Get
qmckl_exit_code qmckl_get_jastrow_champ_een_rescaled_e(qmckl_context context, double* const distance_rescaled, const int64_t size_max);
3.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 |
integer function qmckl_compute_een_rescaled_e_doc_f( & context, walk_num, elec_num, cord_num, rescale_factor_ee, & ee_distance, een_rescaled_e) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: walk_num integer*8 , intent(in) :: elec_num integer*8 , intent(in) :: cord_num double precision , intent(in) :: rescale_factor_ee double precision , intent(in) :: ee_distance(elec_num,elec_num,walk_num) double precision , intent(out) :: een_rescaled_e(elec_num,elec_num,0:cord_num,walk_num) double precision,dimension(:,:),allocatable :: een_rescaled_e_ij double precision :: x integer*8 :: i, j, k, l, nw info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif if (walk_num <= 0) then info = QMCKL_INVALID_ARG_2 return endif if (elec_num <= 0) then info = QMCKL_INVALID_ARG_3 return endif if (cord_num < 0) then info = QMCKL_INVALID_ARG_4 return endif allocate(een_rescaled_e_ij(elec_num * (elec_num - 1) / 2, cord_num + 1)) ! Prepare table of exponentiated distances raised to appropriate power do nw = 1, walk_num een_rescaled_e_ij(:, 1) = 1.0d0 k = 0 do j = 1, elec_num do i = 1, j - 1 k = k + 1 een_rescaled_e_ij(k, 2) = dexp(-rescale_factor_ee * ee_distance(i, j, nw)) end do end do do l = 2, cord_num do k = 1, elec_num * (elec_num - 1)/2 een_rescaled_e_ij(k, l + 1) = een_rescaled_e_ij(k, l) * een_rescaled_e_ij(k, 2) end do end do ! prepare the actual een table een_rescaled_e(:, :, 0, nw) = 1.0d0 do l = 1, cord_num k = 0 do j = 1, elec_num do i = 1, j - 1 k = k + 1 x = een_rescaled_e_ij(k, l + 1) een_rescaled_e(i, j, l, nw) = x een_rescaled_e(j, i, l, nw) = x end do end do end do do l = 0, cord_num do j = 1, elec_num een_rescaled_e(j, j, l, nw) = 0.0d0 end do end do end do end function qmckl_compute_een_rescaled_e_doc_f
qmckl_exit_code qmckl_compute_een_rescaled_e_hpc ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t cord_num, const double rescale_factor_ee, const double* ee_distance, double* const een_rescaled_e ) { if (context == QMCKL_NULL_CONTEXT) { return QMCKL_INVALID_CONTEXT; } if (walk_num <= 0) { return QMCKL_INVALID_ARG_2; } if (elec_num <= 0) { return QMCKL_INVALID_ARG_3; } if (cord_num < 0) { return QMCKL_INVALID_ARG_4; } // Prepare table of exponentiated distances raised to appropriate power // init const size_t elec_pairs = (size_t) (elec_num * (elec_num - 1)) / 2; const size_t len_een_ij = (size_t) elec_pairs * (cord_num + 1); // number of elements for the een_rescaled_e_ij[N_e*(N_e-1)/2][cord+1] // probably in C is better [cord+1, Ne*(Ne-1)/2] // elec_pairs = (elec_num * (elec_num - 1)) / 2; // len_een_ij = elec_pairs * (cord_num + 1); const size_t e2 = elec_num*elec_num; #ifdef HAVE_OPENMP #pragma omp parallel #endif { double* restrict een_rescaled_e_ij = calloc(len_een_ij,sizeof(double)); for (size_t kk = 0; kk < elec_pairs ; ++kk) { een_rescaled_e_ij[kk]= 1.0; } #ifdef HAVE_OPENMP #pragma omp for #endif for (size_t nw = 0; nw < (size_t) walk_num; ++nw) { size_t kk = 0; for (size_t i = 0; i < (size_t) elec_num; ++i) { double* restrict ee1 = &een_rescaled_e_ij[kk + elec_pairs]; const double* restrict ee2 = &ee_distance[i*elec_num + nw*e2]; #ifdef HAVE_OPENMP #pragma omp simd #endif for (size_t j = 0; j < i; ++j) { // een_rescaled_e_ij[j + kk + elec_pairs] = -rescale_factor_ee * ee_distance[j + i*elec_num + nw*e2]; ee1[j] = -rescale_factor_ee * ee2[j]; } kk += i; } #ifdef HAVE_OPENMP #pragma omp simd #endif for (size_t k = elec_pairs; k < 2*elec_pairs; ++k) { een_rescaled_e_ij[k] = exp(een_rescaled_e_ij[k]); } const double* const ee3 = &een_rescaled_e_ij[elec_pairs]; for (size_t l = 2; l < (size_t) (cord_num+1); ++l) { double* restrict ee1 = &een_rescaled_e_ij[l*elec_pairs]; const double* restrict ee2 = &een_rescaled_e_ij[(l-1)*elec_pairs]; #ifdef HAVE_OPENMP #pragma omp simd #endif for (size_t k = 0; k < elec_pairs; ++k) { // een_rescaled_e_ij(k, l + 1) = een_rescaled_e_ij(k, l + 1 - 1) * een_rescaled_e_ij(k, 2) // een_rescaled_e_ij[k+l*elec_pairs] = een_rescaled_e_ij[k + (l - 1)*elec_pairs] * // een_rescaled_e_ij[k + elec_pairs]; ee1[k] = ee2[k] * ee3[k]; } } double* restrict const een_rescaled_e_ = &(een_rescaled_e[nw*(cord_num+1)*e2]); // prepare the actual een table #ifdef HAVE_OPENMP #pragma omp simd #endif for (size_t i = 0; i < e2; ++i){ een_rescaled_e_[i] = 1.0; } for ( size_t l = 1; l < (size_t) (cord_num+1); ++l) { double* x = een_rescaled_e_ij + l*elec_pairs; double* const een_rescaled_e__ = &(een_rescaled_e_[l*e2]); double* een_rescaled_e_i = een_rescaled_e__; for (size_t i = 0; i < (size_t) elec_num; ++i) { for (size_t j = 0; j < i; ++j) { een_rescaled_e_i[j] = *x; een_rescaled_e__[i + j*elec_num] = *x; x += 1; } een_rescaled_e_i += elec_num; } } double* const x0 = &(een_rescaled_e[nw*e2*(cord_num+1)]); for (size_t l = 0; l < (size_t) (cord_num + 1); ++l) { double* x1 = &(x0[l*e2]); for (size_t j = 0; j < (size_t) elec_num; ++j) { *x1 = 0.0; x1 += 1+elec_num; } } } free(een_rescaled_e_ij); } // OpenMP return QMCKL_SUCCESS; }
qmckl_exit_code qmckl_compute_een_rescaled_e_doc ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t cord_num, const double rescale_factor_ee, const double* ee_distance, double* const een_rescaled_e );
qmckl_exit_code qmckl_compute_een_rescaled_e_hpc ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t cord_num, const double rescale_factor_ee, const double* ee_distance, double* const een_rescaled_e );
qmckl_exit_code qmckl_compute_een_rescaled_e ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t cord_num, const double rescale_factor_ee, const double* ee_distance, double* const een_rescaled_e ) { #ifdef HAVE_HPC return qmckl_compute_een_rescaled_e_hpc #else return qmckl_compute_een_rescaled_e_doc #endif (context, walk_num, elec_num, cord_num, rescale_factor_ee, ee_distance, een_rescaled_e); }
3.3.1.3 Test
assert(qmckl_electron_provided(context)); double een_rescaled_e[walk_num][(cord_num + 1)][elec_num][elec_num]; rc = qmckl_get_jastrow_champ_een_rescaled_e(context, &(een_rescaled_e[0][0][0][0]),elec_num*elec_num*(cord_num+1)*walk_num); // value of (0,2,1) assert(fabs(een_rescaled_e[0][1][0][2]- 0.2211015082992776 ) < 1.e-12); assert(fabs(een_rescaled_e[0][1][0][3]- 0.2611178387068169 ) < 1.e-12); assert(fabs(een_rescaled_e[0][1][0][4]- 0.0884012350763747 ) < 1.e-12); assert(fabs(een_rescaled_e[0][2][1][3]- 0.1016685507354656 ) < 1.e-12); assert(fabs(een_rescaled_e[0][2][1][4]- 0.0113118073246869 ) < 1.e-12); assert(fabs(een_rescaled_e[0][2][1][5]- 0.5257156022077619 ) < 1.e-12);
3.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 jastrowchamp.
\[ \frac{\partial}{\partial x} \left[ {g_\text{e}(r)}\right]^p = -\frac{x}{r} \kappa_\text{e}\, p\,\left[ {g_\text{e}(r)}\right]^p \] \[ \Delta \left[ {g_\text{e}(r)}\right]^p = \frac{2}{r} \kappa_\text{e}\, p\,\left[ {g_\text{e}(r)}\right]^p \right] + \left(\frac{\partial}{\partial x}\left[ {g_\text{e}(r)}\right]^p \right)^2 + \left(\frac{\partial}{\partial y}\left[ {g_\text{e}(r)}\right]^p \right)^2 + \left(\frac{\partial}{\partial z}\left[ {g_\text{e}(r)}\right]^p \right)^2 \]
3.3.2.1 Get
qmckl_exit_code qmckl_get_jastrow_champ_een_rescaled_e_gl(qmckl_context context, double* const distance_rescaled, const int64_t size_max);
3.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[walk_num][3][elec_num] |
in | Electron coordinates |
ee_distance |
double[walk_num][elec_num][elec_num] |
in | Electron-electron distances |
een_rescaled_e |
double[walk_num][0:cord_num][elec_num][elec_num] |
in | Electron-electron distances |
een_rescaled_e_gl |
double[walk_num][0:cord_num][elec_num][4][elec_num] |
out | Electron-electron rescaled distances |
integer function qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl_f( & context, walk_num, elec_num, cord_num, rescale_factor_ee, & coord_ee, ee_distance, een_rescaled_e, een_rescaled_e_gl) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: walk_num integer*8 , intent(in) :: elec_num integer*8 , intent(in) :: cord_num double precision , intent(in) :: rescale_factor_ee double precision , intent(in) :: coord_ee(elec_num,3,walk_num) double precision , intent(in) :: ee_distance(elec_num,elec_num,walk_num) double precision , intent(in) :: een_rescaled_e(elec_num,elec_num,0:cord_num,walk_num) double precision , intent(out) :: een_rescaled_e_gl(elec_num,4,elec_num,0:cord_num,walk_num) double precision,dimension(:,:,:),allocatable :: elec_dist_gl double precision :: x, rij_inv, kappa_l integer*8 :: i, j, k, l, nw, ii allocate(elec_dist_gl(elec_num, 4, elec_num)) info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif if (walk_num <= 0) then info = QMCKL_INVALID_ARG_2 return endif if (elec_num <= 0) then info = QMCKL_INVALID_ARG_3 return endif if (cord_num < 0) then info = QMCKL_INVALID_ARG_4 return endif ! Prepare table of exponentiated distances raised to appropriate power do nw = 1, walk_num do j = 1, elec_num do i = 1, j-1 rij_inv = 1.0d0 / ee_distance(i, j, nw) do ii = 1, 3 elec_dist_gl(i, ii, j) = (coord_ee(i, ii, nw) - coord_ee(j, ii, nw)) * rij_inv end do elec_dist_gl(i, 4, j) = 2.0d0 * rij_inv end do elec_dist_gl(j, :, j) = 0.0d0 do i = j+1, elec_num rij_inv = 1.0d0 / ee_distance(i, j, nw) do ii = 1, 3 elec_dist_gl(i, ii, j) = (coord_ee(i, ii, nw) - coord_ee(j, ii, nw)) * rij_inv end do elec_dist_gl(i, 4, j) = 2.0d0 * rij_inv 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 een_rescaled_e_gl(i, 1, j, l, nw) = kappa_l * elec_dist_gl(i, 1, j) een_rescaled_e_gl(i, 2, j, l, nw) = kappa_l * elec_dist_gl(i, 2, j) een_rescaled_e_gl(i, 3, j, l, nw) = kappa_l * elec_dist_gl(i, 3, j) een_rescaled_e_gl(i, 4, j, l, nw) = kappa_l * elec_dist_gl(i, 4, j) een_rescaled_e_gl(i, 4, j, l, nw) = een_rescaled_e_gl(i, 4, j, l, nw) & + een_rescaled_e_gl(i, 1, j, l, nw) * een_rescaled_e_gl(i, 1, j, l, nw) & + een_rescaled_e_gl(i, 2, j, l, nw) * een_rescaled_e_gl(i, 2, j, l, nw) & + een_rescaled_e_gl(i, 3, j, l, nw) * een_rescaled_e_gl(i, 3, j, l, nw) een_rescaled_e_gl(i,1,j,l,nw) = een_rescaled_e_gl(i,1,j,l,nw) * een_rescaled_e(i,j,l,nw) een_rescaled_e_gl(i,2,j,l,nw) = een_rescaled_e_gl(i,2,j,l,nw) * een_rescaled_e(i,j,l,nw) een_rescaled_e_gl(i,3,j,l,nw) = een_rescaled_e_gl(i,3,j,l,nw) * een_rescaled_e(i,j,l,nw) een_rescaled_e_gl(i,4,j,l,nw) = een_rescaled_e_gl(i,4,j,l,nw) * een_rescaled_e(i,j,l,nw) end do end do end do end do end function qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl_f
qmckl_exit_code qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl_hpc ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t cord_num, const double rescale_factor_ee, const double* coord_ee, const double* ee_distance, const double* een_rescaled_e, double* const een_rescaled_e_gl ) { if (context == QMCKL_NULL_CONTEXT) return QMCKL_INVALID_CONTEXT; if (walk_num <= 0) return QMCKL_INVALID_ARG_2; if (elec_num <= 0) return QMCKL_INVALID_ARG_3; if (cord_num < 0) return QMCKL_INVALID_ARG_4; double* restrict elec_dist_gl0 = (double*) calloc(elec_num * elec_num, sizeof(double)); double* restrict elec_dist_gl1 = (double*) calloc(elec_num * elec_num, sizeof(double)); double* restrict elec_dist_gl2 = (double*) calloc(elec_num * elec_num, sizeof(double)); double* restrict elec_dist_gl3 = (double*) calloc(elec_num * elec_num, sizeof(double)); assert (elec_dist_gl0 != NULL); assert (elec_dist_gl1 != NULL); assert (elec_dist_gl2 != NULL); assert (elec_dist_gl3 != NULL); #pragma omp parallel for for (int64_t nw = 0; nw < walk_num; ++nw) { double rij_inv[elec_num]; for (int64_t j = 0; j < elec_num; ++j) { #ifdef HAVE_OPENMP #pragma omp simd #endif for (int64_t i = 0; i < elec_num ; ++i) { rij_inv[i] = ee_distance[i + j * elec_num + nw * elec_num * elec_num] + 1.e-30; } #ifdef HAVE_OPENMP #pragma omp simd #endif for (int64_t i = 0; i < elec_num ; ++i) { rij_inv[i] = 1.0/rij_inv[i]; } rij_inv[j] = 0.; const double xj = coord_ee[j + nw * elec_num * 3]; const double yj = coord_ee[j + elec_num + nw * elec_num * 3]; const double zj = coord_ee[j + 2 * elec_num + nw * elec_num * 3]; #ifdef HAVE_OPENMP #pragma omp simd #endif for (int64_t i = 0; i < elec_num ; ++i) { const double xi = coord_ee[i + nw * elec_num * 3]; const double yi = coord_ee[i + elec_num + nw * elec_num * 3]; const double zi = coord_ee[i + 2 * elec_num + nw * elec_num * 3]; elec_dist_gl0[i + j * elec_num] = rij_inv[i] * (xi-xj); elec_dist_gl1[i + j * elec_num] = rij_inv[i] * (yi-yj); elec_dist_gl2[i + j * elec_num] = rij_inv[i] * (zi-zj); elec_dist_gl3[i + j * elec_num] = rij_inv[i] + rij_inv[i]; } } for (int64_t j = 0; j < elec_num; ++j) { double* restrict eegl = &een_rescaled_e_gl[ elec_num * 4 * (j + elec_num * (cord_num + 1) * nw)]; #ifdef HAVE_OPENMP #pragma omp simd #endif for (int64_t i = 0; i < 4*elec_num; ++i) { eegl[i] = 0.0; } } for (int64_t l = 1; l <= cord_num; ++l) { double kappa_l = - (double)l * rescale_factor_ee; for (int64_t j = 0; j < elec_num; ++j) { double* restrict eegl = &een_rescaled_e_gl[ elec_num * 4 * (j + elec_num * (l + (cord_num + 1) * nw))]; const double* restrict ee = &een_rescaled_e [ elec_num * (j + elec_num * (l + (cord_num + 1) * nw))]; #ifdef HAVE_OPENMP #pragma omp simd #endif for (int64_t i = 0; i < elec_num; ++i) { eegl[i ] = kappa_l * elec_dist_gl0[i + j * elec_num]; eegl[i + elec_num ] = kappa_l * elec_dist_gl1[i + j * elec_num]; eegl[i + elec_num * 2] = kappa_l * elec_dist_gl2[i + j * elec_num]; eegl[i + elec_num * 3] = kappa_l * elec_dist_gl3[i + j * elec_num]; } #ifdef HAVE_OPENMP #pragma omp simd #endif for (int64_t i = 0; i < elec_num; ++i) { eegl[i + elec_num*3] = eegl[i + elec_num*3] + eegl[i] * eegl[i] + eegl[i + elec_num*1] * eegl[i + elec_num*1] + eegl[i + elec_num*2] * eegl[i + elec_num*2]; } #ifdef HAVE_OPENMP #pragma omp simd #endif for (int64_t i = 0; i < elec_num; ++i) { eegl[i ] *= ee[i]; eegl[i + elec_num * 1] *= ee[i]; eegl[i + elec_num * 2] *= ee[i]; eegl[i + elec_num * 3] *= ee[i]; } } } } free(elec_dist_gl0); free(elec_dist_gl1); free(elec_dist_gl2); free(elec_dist_gl3); return QMCKL_SUCCESS; }
qmckl_exit_code qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t cord_num, const double rescale_factor_ee, const double* coord_ee, const double* ee_distance, const double* een_rescaled_e, double* const een_rescaled_e_gl ) { #ifdef HAVE_HPC return qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl_hpc #else return qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl_doc #endif (context, walk_num, elec_num, cord_num, rescale_factor_ee, coord_ee, ee_distance, een_rescaled_e, een_rescaled_e_gl ); }
3.3.2.3 Test
double een_rescaled_e_gl[walk_num][(cord_num + 1)][elec_num][4][elec_num]; size_max=walk_num*(cord_num + 1)*elec_num*4*elec_num; rc = qmckl_get_jastrow_champ_een_rescaled_e_gl(context, &(een_rescaled_e_gl[0][0][0][0][0]),size_max); // value of (0,0,0,2,1) assert(fabs(een_rescaled_e_gl[0][1][0][0][2] + 0.09831391870751387 ) < 1.e-12); assert(fabs(een_rescaled_e_gl[0][1][0][0][3] + 0.017204157459682526 ) < 1.e-12); assert(fabs(een_rescaled_e_gl[0][1][0][0][4] + 0.013345768421098641 ) < 1.e-12); assert(fabs(een_rescaled_e_gl[0][2][1][0][3] + 0.03733086358273962 ) < 1.e-12); assert(fabs(een_rescaled_e_gl[0][2][1][0][4] + 0.004922634822943517 ) < 1.e-12); assert(fabs(een_rescaled_e_gl[0][2][1][0][5] + 0.5416751547830984 ) < 1.e-12);
3.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:
\[ C_{i\alpha,p} = \left[ \exp\left(-\kappa_\alpha\, R_{i\alpha}\right) \right]^p \]
where \(R_{i\alpha}\) is the matrix of electron-nucleus distances.
3.3.3.1 Get
qmckl_exit_code qmckl_get_jastrow_champ_een_rescaled_n(qmckl_context context, double* const distance_rescaled, const int64_t size_max);
3.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 |
integer function qmckl_compute_een_rescaled_n_f( & context, walk_num, elec_num, nucl_num, & type_nucl_num, type_nucl_vector, cord_num, rescale_factor_en, & en_distance, een_rescaled_n) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: walk_num integer*8 , intent(in) :: elec_num integer*8 , intent(in) :: nucl_num integer*8 , intent(in) :: type_nucl_num integer*8 , intent(in) :: type_nucl_vector(nucl_num) integer*8 , intent(in) :: cord_num double precision , intent(in) :: rescale_factor_en(type_nucl_num) double precision , intent(in) :: en_distance(nucl_num,elec_num,walk_num) double precision , intent(out) :: een_rescaled_n(elec_num,nucl_num,0:cord_num,walk_num) double precision :: x integer*8 :: i, a, k, l, nw info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif if (walk_num <= 0) then info = QMCKL_INVALID_ARG_2 return endif if (elec_num <= 0) then info = QMCKL_INVALID_ARG_3 return endif if (nucl_num <= 0) then info = QMCKL_INVALID_ARG_4 return endif if (cord_num < 0) then info = QMCKL_INVALID_ARG_5 return endif do nw = 1, walk_num ! prepare the actual een table een_rescaled_n(:, :, 0, nw) = 1.0d0 do a = 1, nucl_num do i = 1, elec_num een_rescaled_n(i, a, 1, nw) = dexp(-rescale_factor_en(type_nucl_vector(a)+1) * en_distance(a, i, nw)) end do end do do l = 2, cord_num do a = 1, nucl_num do i = 1, elec_num een_rescaled_n(i, a, l, nw) = een_rescaled_n(i, a, l - 1, nw) * een_rescaled_n(i, a, 1, nw) end do end do end do end do end function qmckl_compute_een_rescaled_n_f
/* qmckl_exit_code qmckl_compute_een_rescaled_n ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t type_nucl_num, int64_t* const type_nucl_vector, const int64_t cord_num, const double* rescale_factor_en, const double* en_distance, double* const een_rescaled_n ) { if (context == QMCKL_NULL_CONTEXT) { return QMCKL_INVALID_CONTEXT; } if (walk_num <= 0) { return QMCKL_INVALID_ARG_2; } if (elec_num <= 0) { return QMCKL_INVALID_ARG_3; } if (nucl_num <= 0) { return QMCKL_INVALID_ARG_4; } if (cord_num < 0) { return QMCKL_INVALID_ARG_5; } // Prepare table of exponentiated distances raised to appropriate power for (int i = 0; i < (walk_num*(cord_num+1)*nucl_num*elec_num); ++i) { een_rescaled_n[i] = 1.0; } for (int nw = 0; nw < walk_num; ++nw) { for (int a = 0; a < nucl_num; ++a) { for (int i = 0; i < elec_num; ++i) { een_rescaled_n[i + a*elec_num + nw * elec_num*nucl_num*(cord_num+1)] = 1.0; een_rescaled_n[i + a*elec_num + elec_num*nucl_num + nw*elec_num*nucl_num*(cord_num+1)] = exp(-rescale_factor_en[type_nucl_vector[a]] * en_distance[a + i*nucl_num + nw*elec_num*nucl_num]); } } for (int l = 2; l < (cord_num+1); ++l){ for (int a = 0; a < nucl_num; ++a) { for (int i = 0; i < elec_num; ++i) { een_rescaled_n[i + a*elec_num + l*elec_num*nucl_num + nw*elec_num*nucl_num*(cord_num+1)] = een_rescaled_n[i + a*elec_num + (l-1)*elec_num*nucl_num + nw*elec_num*nucl_num*(cord_num+1)] * een_rescaled_n[i + a*elec_num + elec_num*nucl_num + nw*elec_num*nucl_num*(cord_num+1)]; } } } } return QMCKL_SUCCESS; } */
3.3.3.3 Test
assert(qmckl_electron_provided(context)); double een_rescaled_n[walk_num][(cord_num + 1)][nucl_num][elec_num]; size_max=walk_num*(cord_num + 1)*nucl_num*elec_num; rc = qmckl_get_jastrow_champ_een_rescaled_n(context, &(een_rescaled_n[0][0][0][0]),size_max); // value of (0,2,1) assert(fabs(een_rescaled_n[0][1][0][2]-0.2603169838750542 )< 1.e-12); assert(fabs(een_rescaled_n[0][1][0][3]-0.3016180139679065 )< 1.e-12); assert(fabs(een_rescaled_n[0][1][0][4]-0.10506023826192266)< 1.e-12); assert(fabs(een_rescaled_n[0][2][1][3]-0.9267719759374164 )< 1.e-12); assert(fabs(een_rescaled_n[0][2][1][4]-0.11497585238132658)< 1.e-12); assert(fabs(een_rescaled_n[0][2][1][5]-0.07534033469115217)< 1.e-12);
3.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 jastrowchamp.
\[ \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 \]
3.3.4.1 Get
qmckl_exit_code qmckl_get_jastrow_champ_een_rescaled_n_gl(qmckl_context context, double* const distance_rescaled, const int64_t size_max);
3.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[walk_num][3][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 |
integer function qmckl_compute_jastrow_champ_factor_een_rescaled_n_gl_f( & context, walk_num, elec_num, nucl_num, type_nucl_num, type_nucl_vector, & cord_num, rescale_factor_en, & coord_ee, coord_n, en_distance, een_rescaled_n, een_rescaled_n_gl) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: walk_num integer*8 , intent(in) :: elec_num integer*8 , intent(in) :: nucl_num integer*8 , intent(in) :: type_nucl_num integer*8 , intent(in) :: type_nucl_vector(nucl_num) integer*8 , intent(in) :: cord_num double precision , intent(in) :: rescale_factor_en(type_nucl_num) double precision , intent(in) :: coord_ee(elec_num,3,walk_num) double precision , intent(in) :: coord_n(nucl_num,3) double precision , intent(in) :: en_distance(nucl_num,elec_num,walk_num) double precision , intent(in) :: een_rescaled_n(elec_num,nucl_num,0:cord_num,walk_num) double precision , intent(out) :: een_rescaled_n_gl(elec_num,4,nucl_num,0:cord_num,walk_num) double precision,dimension(:,:,:),allocatable :: elnuc_dist_gl double precision :: x, ria_inv, kappa_l integer*8 :: i, a, k, l, nw, ii allocate(elnuc_dist_gl(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, ii, nw) - 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_f
3.3.4.3 Test
assert(qmckl_electron_provided(context)); double een_rescaled_n_gl[walk_num][(cord_num + 1)][nucl_num][4][elec_num]; size_max=walk_num*(cord_num + 1)*nucl_num*4*elec_num; rc = qmckl_get_jastrow_champ_een_rescaled_n_gl(context, &(een_rescaled_n_gl[0][0][0][0][0]),size_max); // value of (0,2,1) assert(fabs( -0.11234061209936878 - een_rescaled_n_gl[0][1][0][0][2]) < 1.e-12); assert(fabs( 0.0004440109367151707 - een_rescaled_n_gl[0][1][0][0][3]) < 1.e-12); assert(fabs( -0.012868642597346566 - een_rescaled_n_gl[0][1][0][0][4]) < 1.e-12); assert(fabs( 0.08601122289922644 - een_rescaled_n_gl[0][2][1][0][3]) < 1.e-12); assert(fabs( -0.058681563677207206 - een_rescaled_n_gl[0][2][1][0][4]) < 1.e-12); assert(fabs( 0.005359281880312882 - een_rescaled_n_gl[0][2][1][0][5]) < 1.e-12);
3.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.
3.3.5.1 Compute dimcvector
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 |
integer function qmckl_compute_dim_c_vector_f( & context, cord_num, dim_c_vector) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: cord_num integer*8 , intent(out) :: dim_c_vector double precision :: x integer*8 :: i, a, k, l, p, lmax info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif if (cord_num < 0) then info = QMCKL_INVALID_ARG_2 return endif dim_c_vector = 0 do p = 2, cord_num do k = p - 1, 0, -1 if (k .ne. 0) then lmax = p - k else lmax = p - k - 2 endif do l = lmax, 0, -1 if (iand(p - k - l, 1_8) == 1) cycle dim_c_vector = dim_c_vector + 1 end do end do end do end function qmckl_compute_dim_c_vector_f
3.3.5.2 Get
qmckl_exit_code qmckl_get_jastrow_champ_tmp_c(qmckl_context context, double* const tmp_c); qmckl_exit_code qmckl_get_jastrow_champ_dtmp_c(qmckl_context context, double* const dtmp_c);
3.3.5.3 Compute cvectorfull
| 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[nucl_num][dim_c_vector] |
out | Full list of coefficients |
integer function qmckl_compute_c_vector_full_doc_f( & context, nucl_num, dim_c_vector, type_nucl_num, & type_nucl_vector, c_vector, c_vector_full) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: nucl_num integer*8 , intent(in) :: dim_c_vector integer*8 , intent(in) :: type_nucl_num integer*8 , intent(in) :: type_nucl_vector(nucl_num) double precision , intent(in) :: c_vector(dim_c_vector, type_nucl_num) double precision , intent(out) :: c_vector_full(nucl_num, dim_c_vector) double precision :: x integer*8 :: i, a, k, l, nw info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT if (nucl_num <= 0) info = QMCKL_INVALID_ARG_2 if (dim_c_vector < 0) info = QMCKL_INVALID_ARG_3 if (type_nucl_num <= 0) info = QMCKL_INVALID_ARG_4 if (info /= QMCKL_SUCCESS) return do a = 1, nucl_num c_vector_full(a,1:dim_c_vector) = c_vector(1:dim_c_vector, type_nucl_vector(a)+1) end do end function qmckl_compute_c_vector_full_doc_f
qmckl_exit_code qmckl_compute_c_vector_full_hpc ( const qmckl_context context, const int64_t nucl_num, const int64_t dim_c_vector, const int64_t type_nucl_num, const int64_t* type_nucl_vector, const double* c_vector, double* const c_vector_full ) { if (context == QMCKL_NULL_CONTEXT) return QMCKL_INVALID_CONTEXT; if (nucl_num <= 0) return QMCKL_INVALID_ARG_2; if (dim_c_vector < 0) return QMCKL_INVALID_ARG_3; if (type_nucl_num <= 0) return QMCKL_INVALID_ARG_4; if (type_nucl_vector == NULL) return QMCKL_INVALID_ARG_5; if (c_vector == NULL) return QMCKL_INVALID_ARG_6; if (c_vector_full == NULL) return QMCKL_INVALID_ARG_7; for (int i=0; i < dim_c_vector; ++i) { for (int a=0; a < nucl_num; ++a){ c_vector_full[a + i*nucl_num] = c_vector[i + type_nucl_vector[a]*dim_c_vector]; } } return QMCKL_SUCCESS; }
qmckl_exit_code qmckl_compute_c_vector_full_doc ( const qmckl_context context, const int64_t nucl_num, const int64_t dim_c_vector, const int64_t type_nucl_num, const int64_t* type_nucl_vector, const double* c_vector, double* const c_vector_full );
qmckl_exit_code qmckl_compute_c_vector_full_hpc ( const qmckl_context context, const int64_t nucl_num, const int64_t dim_c_vector, const int64_t type_nucl_num, const int64_t* type_nucl_vector, const double* c_vector, double* const c_vector_full );
qmckl_exit_code qmckl_compute_c_vector_full ( const qmckl_context context, const int64_t nucl_num, const int64_t dim_c_vector, const int64_t type_nucl_num, const int64_t* type_nucl_vector, const double* c_vector, double* const c_vector_full ) { #ifdef HAVE_HPC return qmckl_compute_c_vector_full_hpc #else return qmckl_compute_c_vector_full_doc #endif (context, nucl_num, dim_c_vector, type_nucl_num, type_nucl_vector, c_vector, c_vector_full); }
3.3.5.4 Compute lkpmcombinedindex
| Variable | Type | In/Out | Description |
|---|---|---|---|
context |
qmckl_context |
in | Global state |
cord_num |
int64_t |
in | Order of polynomials |
dim_c_vector |
int64_t |
in | dimension of cord full table |
lkpm_combined_index |
int64_t[4][dim_c_vector] |
out | Full list of combined indices |
integer function qmckl_compute_lkpm_combined_index_doc_f( & context, cord_num, dim_c_vector, lkpm_combined_index) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: cord_num integer*8 , intent(in) :: dim_c_vector integer*8 , intent(out) :: lkpm_combined_index(dim_c_vector, 4) double precision :: x integer*8 :: i, a, k, l, kk, p, lmax, m info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT if (cord_num < 0) info = QMCKL_INVALID_ARG_2 if (dim_c_vector < 0) info = QMCKL_INVALID_ARG_3 if (info /= QMCKL_SUCCESS) return kk = 0 do p = 2, cord_num do k = p - 1, 0, -1 if (k /= 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_f
qmckl_exit_code qmckl_compute_lkpm_combined_index_hpc ( const qmckl_context context, const int64_t cord_num, const int64_t dim_c_vector, int64_t* const lkpm_combined_index ) { int kk, lmax, m; if (context == QMCKL_NULL_CONTEXT) return QMCKL_INVALID_CONTEXT; if (cord_num < 0) return QMCKL_INVALID_ARG_2; if (dim_c_vector < 0) return QMCKL_INVALID_ARG_3; kk = 0; for (int p = 2; p <= cord_num; ++p) { for (int k=(p-1); k >= 0; --k) { if (k != 0) { lmax = p - k; } else { lmax = p - k - 2; } for (int l=lmax; l >= 0; --l) { if (((p - k - l) & 1) == 1) continue; m = (p - k - l)/2; lkpm_combined_index[kk ] = l; lkpm_combined_index[kk + dim_c_vector] = k; lkpm_combined_index[kk + 2*dim_c_vector] = p; lkpm_combined_index[kk + 3*dim_c_vector] = m; kk = kk + 1; } } } return QMCKL_SUCCESS; }
qmckl_exit_code qmckl_compute_lkpm_combined_index ( const qmckl_context context, const int64_t cord_num, const int64_t dim_c_vector, int64_t* const lkpm_combined_index ) { #ifdef HAVE_HPC return qmckl_compute_lkpm_combined_index_hpc #else return qmckl_compute_lkpm_combined_index_doc #endif (context, cord_num, dim_c_vector, lkpm_combined_index); }
qmckl_exit_code qmckl_compute_lkpm_combined_index ( const qmckl_context context, const int64_t cord_num, const int64_t dim_c_vector, int64_t* const lkpm_combined_index );
qmckl_exit_code qmckl_compute_lkpm_combined_index_doc ( const qmckl_context context, const int64_t cord_num, const int64_t dim_c_vector, int64_t* const lkpm_combined_index );
qmckl_exit_code qmckl_compute_lkpm_combined_index_hpc ( const qmckl_context context, const int64_t cord_num, const int64_t dim_c_vector, int64_t* const lkpm_combined_index );
qmckl_exit_code qmckl_compute_lkpm_combined_index ( const qmckl_context context, const int64_t cord_num, const int64_t dim_c_vector, int64_t* const lkpm_combined_index );
3.3.5.5 Compute tmpc
| Variable | Type | In/Out | Description |
|---|---|---|---|
context |
qmckl_context |
in | Global state |
cord_num |
int64_t |
in | Order of polynomials |
elec_num |
int64_t |
in | Number of electrons |
nucl_num |
int64_t |
in | Number of nuclei |
walk_num |
int64_t |
in | Number of walkers |
een_rescaled_e |
double[walk_num][0:cord_num][elec_num][elec_num] |
in | Electron-electron rescaled factor |
een_rescaled_n |
double[walk_num][0:cord_num][nucl_num][elec_num] |
in | Electron-nucleus rescaled factor |
tmp_c |
double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num] |
out | vector of non-zero coefficients |
qmckl_exit_code qmckl_compute_tmp_c (const qmckl_context context, const int64_t cord_num, const int64_t elec_num, const int64_t nucl_num, const int64_t walk_num, const double* een_rescaled_e, const double* een_rescaled_n, double* const tmp_c ) { #ifdef HAVE_HPC return qmckl_compute_tmp_c_hpc #else return qmckl_compute_tmp_c_doc #endif (context, cord_num, elec_num, nucl_num, walk_num, een_rescaled_e, een_rescaled_n, tmp_c); }
integer function qmckl_compute_tmp_c_doc_f( & context, cord_num, elec_num, nucl_num, & walk_num, een_rescaled_e, een_rescaled_n, tmp_c) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: cord_num integer*8 , intent(in) :: elec_num integer*8 , intent(in) :: nucl_num integer*8 , intent(in) :: walk_num double precision , intent(in) :: een_rescaled_e(elec_num, elec_num, 0:cord_num, walk_num) double precision , intent(in) :: een_rescaled_n(elec_num, nucl_num, 0:cord_num, walk_num) double precision , intent(out) :: tmp_c(elec_num, nucl_num,0:cord_num, 0:cord_num-1, walk_num) double precision :: x integer*8 :: i, j, a, l, kk, p, lmax, nw character :: TransA, TransB double precision :: alpha, beta integer*8 :: M, N, K, LDA, LDB, LDC TransA = 'N' TransB = 'N' alpha = 1.0d0 beta = 0.0d0 info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT if (cord_num < 0) info = QMCKL_INVALID_ARG_2 if (elec_num <= 0) info = QMCKL_INVALID_ARG_3 if (nucl_num <= 0) info = QMCKL_INVALID_ARG_4 if (walk_num <= 0) info = QMCKL_INVALID_ARG_5 if (info /= QMCKL_SUCCESS) return M = elec_num N = nucl_num*(cord_num + 1) K = elec_num LDA = size(een_rescaled_e,1) LDB = size(een_rescaled_n,1) LDC = size(tmp_c,1) do nw=1, walk_num do i=0, cord_num-1 info = qmckl_dgemm(context, TransA, TransB, M, N, K, alpha, & een_rescaled_e(1,1,i,nw),LDA*1_8, & een_rescaled_n(1,1,0,nw),LDB*1_8, & beta, & tmp_c(1,1,0,i,nw),LDC) end do end do end function qmckl_compute_tmp_c_doc_f
qmckl_exit_code qmckl_compute_tmp_c_doc ( const qmckl_context context, const int64_t cord_num, const int64_t elec_num, const int64_t nucl_num, const int64_t walk_num, const double* een_rescaled_e, const double* een_rescaled_n, double* const tmp_c );
3.3.5.6 Compute dtmpc
| Variable | Type | In/Out | Description |
|---|---|---|---|
context |
qmckl_context |
in | Global state |
cord_num |
int64_t |
in | Order of polynomials |
elec_num |
int64_t |
in | Number of electrons |
nucl_num |
int64_t |
in | Number of nuclei |
walk_num |
int64_t |
in | Number of walkers |
een_rescaled_e_gl |
double[walk_num][0:cord_num][elec_num][4][elec_num] |
in | Electron-electron rescaled factor derivatives |
een_rescaled_n |
double[walk_num][0:cord_num][nucl_num][elec_num] |
in | Electron-nucleus rescaled factor |
dtmp_c |
double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num] |
out | vector of non-zero coefficients |
qmckl_exit_code qmckl_compute_dtmp_c (const qmckl_context context, const int64_t cord_num, const int64_t elec_num, const int64_t nucl_num, const int64_t walk_num, const double* een_rescaled_e_gl, const double* een_rescaled_n, double* const dtmp_c ) { #ifdef HAVE_HPC return qmckl_compute_dtmp_c_hpc #else return qmckl_compute_dtmp_c_doc #endif (context, cord_num, elec_num, nucl_num, walk_num, een_rescaled_e_gl, een_rescaled_n, dtmp_c ); }
integer function qmckl_compute_dtmp_c_doc_f( & context, cord_num, elec_num, nucl_num, & walk_num, een_rescaled_e_gl, een_rescaled_n, dtmp_c) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: cord_num integer*8 , intent(in) :: elec_num integer*8 , intent(in) :: nucl_num integer*8 , intent(in) :: walk_num double precision , intent(in) :: een_rescaled_e_gl(elec_num, 4, elec_num, 0:cord_num, walk_num) double precision , intent(in) :: een_rescaled_n(elec_num, nucl_num, 0:cord_num, walk_num) double precision , intent(out) :: dtmp_c(elec_num, 4, nucl_num,0:cord_num, 0:cord_num-1, walk_num) double precision :: x integer*8 :: i, j, a, l, kk, p, lmax, nw, ii character :: TransA, TransB double precision :: alpha, beta integer*8 :: M, N, K, LDA, LDB, LDC info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT if (cord_num < 0) info = QMCKL_INVALID_ARG_2 if (elec_num <= 0) info = QMCKL_INVALID_ARG_3 if (nucl_num <= 0) info = QMCKL_INVALID_ARG_4 if (walk_num <= 0) info = QMCKL_INVALID_ARG_5 if (info /= QMCKL_SUCCESS) return TransA = 'N' TransB = 'N' alpha = 1.0d0 beta = 0.0d0 M = 4*elec_num N = nucl_num*(cord_num + 1) K = elec_num LDA = 4*size(een_rescaled_e_gl,1) LDB = size(een_rescaled_n,1) LDC = 4*size(dtmp_c,1) do nw=1, walk_num do i=0, cord_num-1 info = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & een_rescaled_e_gl(1,1,1,i,nw),LDA*1_8, & een_rescaled_n(1,1,0,nw),LDB*1_8, & beta, & dtmp_c(1,1,1,0,i,nw),LDC) end do end do end function qmckl_compute_dtmp_c_doc_f
3.3.5.7 Test
assert(qmckl_electron_provided(context)); double tmp_c[walk_num][cord_num][cord_num+1][nucl_num][elec_num]; rc = qmckl_get_jastrow_champ_tmp_c(context, &(tmp_c[0][0][0][0][0])); double dtmp_c[walk_num][cord_num][cord_num+1][nucl_num][4][elec_num]; rc = qmckl_get_jastrow_champ_dtmp_c(context, &(dtmp_c[0][0][0][0][0][0])); printf("%e\n%e\n", tmp_c[0][0][1][0][0], 3.954384); assert(fabs(tmp_c[0][0][1][0][0] - 3.954384) < 1e-6); printf("%e\n%e\n", dtmp_c[0][1][0][0][0][0],3.278657e-01); assert(fabs(dtmp_c[0][1][0][0][0][0] - 3.278657e-01 ) < 1e-6);
3.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.
3.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
3.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 |
integer function qmckl_compute_jastrow_champ_factor_een_naive_f( & context, walk_num, elec_num, nucl_num, cord_num,& dim_c_vector, c_vector_full, lkpm_combined_index, & een_rescaled_e, een_rescaled_n, factor_een) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: walk_num, elec_num, cord_num, nucl_num, dim_c_vector integer*8 , intent(in) :: lkpm_combined_index(dim_c_vector,4) double precision , intent(in) :: c_vector_full(nucl_num, dim_c_vector) double precision , intent(in) :: een_rescaled_e(elec_num, elec_num, 0:cord_num, walk_num) double precision , intent(in) :: een_rescaled_n(elec_num, nucl_num, 0:cord_num, walk_num) double precision , intent(out) :: factor_een(walk_num) integer*8 :: i, a, j, l, k, p, m, n, nw double precision :: accu, accu2, cn info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT if (walk_num <= 0) info = QMCKL_INVALID_ARG_2 if (elec_num <= 0) info = QMCKL_INVALID_ARG_3 if (nucl_num <= 0) info = QMCKL_INVALID_ARG_4 if (cord_num < 0) info = QMCKL_INVALID_ARG_5 if (info /= QMCKL_SUCCESS) return ! do nw =1, walk_num ! factor_een(nw) = 0.0d0 ! do n = 1, dim_c_vector ! l = lkpm_combined_index(n, 1) ! k = lkpm_combined_index(n, 2) ! p = lkpm_combined_index(n, 3) ! m = lkpm_combined_index(n, 4) ! do a = 1, nucl_num ! accu2 = 0.0d0 ! cn = c_vector_full(a, n) ! do j = 1, elec_num ! accu = 0.0d0 ! do i = 1, elec_num ! accu = accu + een_rescaled_e(i,j,k,nw) * & ! een_rescaled_n(i,a,m,nw) ! end do ! accu2 = accu2 + accu * een_rescaled_n(j,a,m + l,nw) ! end do ! factor_een(nw) = factor_een(nw) + accu2 * cn ! end do ! end do ! end do 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 accu2 = 0.0d0 cn = c_vector_full(a, n) print *, a, l, k, p, cn do j = 1, elec_num accu = 0.0d0 do i = 1, j-1 accu = accu + 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 accu2 = accu2 + accu end do factor_een(nw) = factor_een(nw) + accu2 * cn end do end do end do end function qmckl_compute_jastrow_champ_factor_een_naive_f
3.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 |
integer function qmckl_compute_jastrow_champ_factor_een_doc_f( & context, walk_num, elec_num, nucl_num, cord_num, & dim_c_vector, c_vector_full, lkpm_combined_index, & tmp_c, een_rescaled_n, factor_een) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: walk_num, elec_num, cord_num, nucl_num, dim_c_vector integer*8 , intent(in) :: lkpm_combined_index(dim_c_vector,4) double precision , intent(in) :: c_vector_full(nucl_num, dim_c_vector) double precision , intent(in) :: tmp_c(elec_num, nucl_num,0:cord_num, 0:cord_num-1, walk_num) double precision , intent(in) :: een_rescaled_n(elec_num, nucl_num, 0:cord_num, walk_num) double precision , intent(out) :: factor_een(walk_num) integer*8 :: i, a, j, l, k, p, m, n, nw double precision :: accu, accu2, cn info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT if (walk_num <= 0) info = QMCKL_INVALID_ARG_2 if (elec_num <= 0) info = QMCKL_INVALID_ARG_3 if (nucl_num <= 0) info = QMCKL_INVALID_ARG_4 if (cord_num < 0) info = QMCKL_INVALID_ARG_5 if (info /= QMCKL_SUCCESS) return factor_een = 0.0d0 if (cord_num == 0) return do nw =1, walk_num do n = 1, dim_c_vector l = lkpm_combined_index(n, 1) k = lkpm_combined_index(n, 2) p = lkpm_combined_index(n, 3) m = lkpm_combined_index(n, 4) do a = 1, nucl_num cn = c_vector_full(a, n) if(cn == 0.d0) cycle accu = 0.0d0 do j = 1, elec_num accu = accu + een_rescaled_n(j,a,m,nw) * tmp_c(j,a,m+l,k,nw) end do factor_een(nw) = factor_een(nw) + accu * cn end do end do end do end function qmckl_compute_jastrow_champ_factor_een_doc_f
3.3.6.4 Test
/* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); double factor_een[walk_num]; rc = qmckl_get_jastrow_champ_factor_een(context, &(factor_een[0]),walk_num); assert(fabs(factor_een[0] + 0.382580260174321) < 1e-12);
3.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.
3.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);
- 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
3.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 |
integer function qmckl_compute_jastrow_champ_factor_een_gl_naive_f( & context, walk_num, elec_num, nucl_num, cord_num, dim_c_vector, & c_vector_full, lkpm_combined_index, een_rescaled_e, een_rescaled_n, & een_rescaled_e_gl, een_rescaled_n_gl, factor_een_gl)& result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: walk_num, elec_num, cord_num, nucl_num, dim_c_vector integer*8 , intent(in) :: lkpm_combined_index(dim_c_vector, 4) double precision , intent(in) :: c_vector_full(nucl_num, dim_c_vector) double precision , intent(in) :: een_rescaled_e(elec_num, elec_num, 0:cord_num, walk_num) double precision , intent(in) :: een_rescaled_n(elec_num, nucl_num, 0:cord_num, walk_num) double precision , intent(in) :: een_rescaled_e_gl(elec_num, 4, elec_num, 0:cord_num, walk_num) double precision , intent(in) :: een_rescaled_n_gl(elec_num, 4, nucl_num, 0:cord_num, walk_num) double precision , intent(out) :: factor_een_gl(elec_num, 4, walk_num) integer*8 :: i, a, j, l, k, p, m, n, nw double precision :: accu, accu2, cn double precision :: daccu(1:4), daccu2(1:4) info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT if (walk_num <= 0) info = QMCKL_INVALID_ARG_2 if (elec_num <= 0) info = QMCKL_INVALID_ARG_3 if (nucl_num <= 0) info = QMCKL_INVALID_ARG_4 if (cord_num < 0) info = QMCKL_INVALID_ARG_5 if (info /= QMCKL_SUCCESS) return factor_een_gl = 0.0d0 do nw =1, walk_num do n = 1, dim_c_vector l = lkpm_combined_index(n, 1) k = lkpm_combined_index(n, 2) p = lkpm_combined_index(n, 3) m = lkpm_combined_index(n, 4) do a = 1, nucl_num cn = c_vector_full(a, n) do j = 1, elec_num accu = 0.0d0 accu2 = 0.0d0 daccu = 0.0d0 daccu2 = 0.0d0 do i = 1, elec_num accu = accu + een_rescaled_e(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_f
3.3.7.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 | Temporary intermediate tensor |
dtmp_c |
double[walk_num][0:cord_num-1][0:cord_num][nucl_num][4][elec_num] |
in | vector of non-zero coefficients |
een_rescaled_n |
double[walk_num][0:cord_num][nucl_num][elec_num] |
in | Electron-nucleus rescaled factor |
een_rescaled_n_gl |
double[walk_num][0:cord_num][nucl_num][4][elec_num] |
in | Derivative of Electron-nucleus rescaled factor |
factor_een_gl |
double[walk_num][4][elec_num] |
out | Derivative of Electron-nucleus jastrow |
integer function qmckl_compute_jastrow_champ_factor_een_gl_doc_f( & context, walk_num, elec_num, nucl_num, & cord_num, dim_c_vector, c_vector_full, lkpm_combined_index, & tmp_c, dtmp_c, een_rescaled_n, een_rescaled_n_gl, factor_een_gl)& result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: walk_num, elec_num, cord_num, nucl_num, dim_c_vector integer*8 , intent(in) :: lkpm_combined_index(dim_c_vector,4) double precision , intent(in) :: c_vector_full(nucl_num, dim_c_vector) double precision , intent(in) :: tmp_c(elec_num, nucl_num,0:cord_num, 0:cord_num-1, walk_num) double precision , intent(in) :: dtmp_c(elec_num, 4, nucl_num,0:cord_num, 0:cord_num-1, walk_num) double precision , intent(in) :: een_rescaled_n(elec_num, nucl_num, 0:cord_num, walk_num) double precision , intent(in) :: een_rescaled_n_gl(elec_num, 4, nucl_num, 0:cord_num, walk_num) double precision , intent(out) :: factor_een_gl(elec_num,4,walk_num) integer*8 :: i, a, j, l, k, m, n, nw, ii double precision :: accu, accu2, cn info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT if (walk_num <= 0) info = QMCKL_INVALID_ARG_2 if (elec_num <= 0) info = QMCKL_INVALID_ARG_3 if (nucl_num <= 0) info = QMCKL_INVALID_ARG_4 if (cord_num < 0) info = QMCKL_INVALID_ARG_5 if (info /= QMCKL_SUCCESS) return factor_een_gl = 0.0d0 if (cord_num == 0) return do nw =1, walk_num do n = 1, dim_c_vector l = lkpm_combined_index(n, 1) k = lkpm_combined_index(n, 2) m = lkpm_combined_index(n, 4) do a = 1, nucl_num cn = c_vector_full(a, n) if(cn == 0.d0) cycle do ii = 1, 4 do j = 1, elec_num factor_een_gl(j,ii,nw) = factor_een_gl(j,ii,nw) + ( & tmp_c(j,a,m,k,nw) * een_rescaled_n_gl(j,ii,a,m+l,nw) + & (dtmp_c(j,ii,a,m,k,nw)) * een_rescaled_n(j,a,m+l,nw) + & (dtmp_c(j,ii,a,m+l,k,nw)) * een_rescaled_n(j,a,m ,nw) + & tmp_c(j,a,m+l,k,nw) * een_rescaled_n_gl(j,ii,a,m,nw) & ) * cn end do end do cn = cn + cn do j = 1, elec_num factor_een_gl(j,4,nw) = factor_een_gl(j,4,nw) + ( & (dtmp_c(j,1,a,m ,k,nw)) * een_rescaled_n_gl(j,1,a,m+l,nw) + & (dtmp_c(j,2,a,m ,k,nw)) * een_rescaled_n_gl(j,2,a,m+l,nw) + & (dtmp_c(j,3,a,m ,k,nw)) * een_rescaled_n_gl(j,3,a,m+l,nw) + & (dtmp_c(j,1,a,m+l,k,nw)) * een_rescaled_n_gl(j,1,a,m ,nw) + & (dtmp_c(j,2,a,m+l,k,nw)) * een_rescaled_n_gl(j,2,a,m ,nw) + & (dtmp_c(j,3,a,m+l,k,nw)) * een_rescaled_n_gl(j,3,a,m ,nw) & ) * cn end do end do end do end do end function qmckl_compute_jastrow_champ_factor_een_gl_doc_f
qmckl_exit_code qmckl_compute_jastrow_champ_factor_een_gl_doc ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, const int64_t dim_c_vector, const double* c_vector_full, const int64_t* lkpm_combined_index, const double* tmp_c, const double* dtmp_c, const double* een_rescaled_n, const double* een_rescaled_n_gl, double* const factor_een_gl );
qmckl_exit_code qmckl_compute_jastrow_champ_factor_een_gl ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, const int64_t dim_c_vector, const double* c_vector_full, const int64_t* lkpm_combined_index, const double* tmp_c, const double* dtmp_c, const double* een_rescaled_n, const double* een_rescaled_n_gl, double* const factor_een_gl );
qmckl_exit_code qmckl_compute_jastrow_champ_factor_een_gl(const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, const int64_t dim_c_vector, const double *c_vector_full, const int64_t *lkpm_combined_index, const double *tmp_c, const double *dtmp_c, const double *een_rescaled_n, const double *een_rescaled_n_gl, double* const factor_een_gl) { #ifdef HAVE_HPC return qmckl_compute_jastrow_champ_factor_een_gl_hpc #else return qmckl_compute_jastrow_champ_factor_een_gl_doc #endif (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); }
3.3.7.4 Test
/* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); double factor_een_gl[4][walk_num][elec_num]; rc = qmckl_get_jastrow_champ_factor_een_gl(context, &(factor_een_gl[0][0][0]),4*walk_num*elec_num); printf("%20.15e\n", factor_een_gl[0][0][0]); assert(fabs(8.967809309100624e-02 - factor_een_gl[0][0][0]) < 1e-12); printf("%20.15e\n", factor_een_gl[1][0][1]); assert(fabs(3.543090132452453e-02 - factor_een_gl[1][0][1]) < 1e-12); printf("%20.15e\n", factor_een_gl[2][0][2]); assert(fabs(8.996044894431991e-04 - factor_een_gl[2][0][2]) < 1e-12); printf("%20.15e\n", factor_een_gl[3][0][3]); assert(fabs(-1.175028308456619e+00 - factor_een_gl[3][0][3]) < 1e-12);
3.4 Total Jastrow
3.4.1 Value
Value of the total Jastrow factor: \(\exp(J)\)
3.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
3.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 |
integer function qmckl_compute_jastrow_champ_value_doc_f(context, & walk_num, f_ee, f_en, f_een, value) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: walk_num double precision , intent(in) :: f_ee(walk_num), f_en(walk_num), f_een(walk_num) double precision , intent(out) :: value(walk_num) integer*8 :: i info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif if (walk_num <= 0) then info = QMCKL_INVALID_ARG_2 return endif do i = 1, walk_num value(i) = f_ee(i) + f_en(i) + f_een(i) end do do i = 1, walk_num ! Flush to zero to avoid floating-point exception if (value(i) < -100.d0) then value(i) = 0.d0 else value(i) = dexp(value(i)) endif end do end function qmckl_compute_jastrow_champ_value_doc_f
qmckl_exit_code qmckl_compute_jastrow_champ_value ( const qmckl_context context, const int64_t walk_num, const double* f_ee, const double* f_en, const double* f_een, double* const value );
qmckl_exit_code qmckl_compute_jastrow_champ_value_doc ( const qmckl_context context, const int64_t walk_num, const double* f_ee, const double* f_en, const double* f_een, double* const value );
qmckl_exit_code qmckl_compute_jastrow_champ_value_hpc ( const qmckl_context context, const int64_t walk_num, const double* f_ee, const double* f_en, const double* f_een, double* const value );
qmckl_exit_code qmckl_compute_jastrow_champ_value ( const qmckl_context context, const int64_t walk_num, const double* factor_ee, const double* factor_en, const double* factor_een, double* const value) { #ifdef HAVE_HPC return qmckl_compute_jastrow_champ_value_hpc #else return qmckl_compute_jastrow_champ_value_doc #endif (context, walk_num, factor_ee, factor_en, factor_een, value); }
3.4.1.3 Test
printf("Total Jastrow value\n"); /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); rc = qmckl_check(context, qmckl_get_jastrow_champ_factor_ee(context, &(factor_ee[0]), walk_num) ); assert(rc == QMCKL_SUCCESS); rc = qmckl_check(context, qmckl_get_jastrow_champ_factor_en(context, &(factor_en[0]), walk_num) ); assert(rc == QMCKL_SUCCESS); rc = qmckl_check(context, qmckl_get_jastrow_champ_factor_een(context, &(factor_een[0]), walk_num) ); assert(rc == QMCKL_SUCCESS); double total_j[walk_num]; rc = qmckl_check(context, qmckl_get_jastrow_champ_value(context, &(total_j[0]), walk_num) ); assert(rc == QMCKL_SUCCESS); for (int64_t i=0 ; i< walk_num ; ++i) { assert (total_j[i] - exp(factor_ee[i] + factor_en[i] + factor_een[i]) < 1.e-12); }
3.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] \]
3.4.2.1 Get
qmckl_exit_code qmckl_get_jastrow_champ_gl(qmckl_context context, double* const gl, const int64_t size_max);
- Fortran interface
interface integer(qmckl_exit_code) function qmckl_get_jastrow_champ_gl (context, & gl, size_max) bind(C) use, intrinsic :: iso_c_binding import implicit none integer (qmckl_context) , intent(in), value :: context integer(c_int64_t), intent(in), value :: size_max real(c_double), intent(out) :: gl(size_max) end function qmckl_get_jastrow_champ_gl end interface
3.4.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 |
value |
double[walk_num] |
in | Total Jastrow |
gl_ee |
double[walk_num][4][elec_num] |
in | ee component |
gl_en |
double[walk_num][4][elec_num] |
in | eN component |
gl_een |
double[walk_num][4][elec_num] |
in | eeN component |
gl |
double[walk_num][4][elec_num] |
out | Total Jastrow factor |
integer function qmckl_compute_jastrow_champ_gl_doc_f(context, & walk_num, elec_num, value, gl_ee, gl_en, gl_een, gl) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: walk_num, elec_num double precision , intent(in) :: value (walk_num) double precision , intent(in) :: gl_ee (elec_num,4,walk_num) double precision , intent(in) :: gl_en (elec_num,4,walk_num) double precision , intent(in) :: gl_een(elec_num,4,walk_num) double precision , intent(out) :: gl (elec_num,4,walk_num) integer*8 :: i, j, k info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif if (walk_num <= 0) then info = QMCKL_INVALID_ARG_2 return endif do k = 1, walk_num do j=1,4 do i = 1, elec_num gl(i,j,k) = gl_ee(i,j,k) + gl_en(i,j,k) + gl_een(i,j,k) end do end do do i = 1, elec_num gl(i,4,k) = gl(i,4,k) + & gl(i,1,k) * gl(i,1,k) + & gl(i,2,k) * gl(i,2,k) + & gl(i,3,k) * gl(i,3,k) end do gl(:,:,k) = gl(:,:,k) * value(k) end do end function qmckl_compute_jastrow_champ_gl_doc_f
qmckl_exit_code qmckl_compute_jastrow_champ_gl ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const double* value, const double* gl_ee, const double* gl_en, const double* gl_een, double* const gl );
qmckl_exit_code qmckl_compute_jastrow_champ_gl_doc ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const double* value, const double* gl_ee, const double* gl_en, const double* gl_een, double* const gl );
qmckl_exit_code qmckl_compute_jastrow_champ_gl_hpc ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const double* value, const double* gl_ee, const double* gl_en, const double* gl_een, double* const gl );
qmckl_exit_code qmckl_compute_jastrow_champ_gl ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const double* value, const double* gl_ee, const double* gl_en, const double* gl_een, double* const gl) { #ifdef HAVE_HPC return qmckl_compute_jastrow_champ_gl_hpc #else return qmckl_compute_jastrow_champ_gl_doc #endif (context, walk_num, elec_num, value, gl_ee, gl_en, gl_een, gl); }
3.4.2.3 Test
printf("Total Jastrow derivatives\n"); /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); rc = qmckl_check(context, qmckl_get_jastrow_champ_factor_ee_gl(context, &(factor_ee_gl[0][0][0]), walk_num*elec_num*4) ); assert(rc == QMCKL_SUCCESS); rc = qmckl_check(context, qmckl_get_jastrow_champ_factor_en_gl(context, &(factor_en_gl[0][0][0]), walk_num*elec_num*4) ); assert(rc == QMCKL_SUCCESS); rc = qmckl_check(context, qmckl_get_jastrow_champ_factor_een_gl(context, &(factor_een_gl[0][0][0]), walk_num*elec_num*4) ); assert(rc == QMCKL_SUCCESS); double total_j_deriv[walk_num][4][elec_num]; rc = qmckl_check(context, qmckl_get_jastrow_champ_gl(context, &(total_j_deriv[0][0][0]), walk_num*elec_num*4) ); assert(rc == QMCKL_SUCCESS); rc = qmckl_check(context, qmckl_get_jastrow_champ_value(context, &(total_j[0]), walk_num) ); assert(rc == QMCKL_SUCCESS); for (int64_t k=0 ; k< walk_num ; ++k) { for (int64_t m=0 ; m<4; ++m) { for (int64_t e=0 ; e<elec_num; ++e) { if (m < 3) { /* test only gradients */ assert (total_j_deriv[k][m][e]/total_j[k] - (factor_ee_gl[k][m][e] + factor_en_gl[k][m][e] + factor_een_gl[k][m][e]) < 1.e-12); } } } }