Jastrow Factor
Table of Contents
- 1. Introduction
- 2. Context
- 3. Computation
- 3.1. Asymptotic component for \(J_{ee}\)
- 3.2. Electron-electron component \(f_{ee}\)
- 3.3. Electron-electron component derivative \(f'_{ee}\)
- 3.4. Electron-nucleus component \(f_{en}\)
- 3.5. Electron-nucleus component derivative \(f'_{en}\)
- 3.6. Electron-electron rescaled distances for each order
- 3.7. Electron-electron rescaled distances for each order and derivatives
- 3.8. Electron-nucleus rescaled distances for each order
- 3.9. Electron-nucleus rescaled distances for each order and derivatives
- 3.10. Prepare for electron-electron-nucleus Jastrow \(f_{een}\)
- 3.11. Electron-electron-nucleus Jastrow \(f_{een}\)
- 3.12. Electron-electron-nucleus Jastrow \(f_{een}\) derivative
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 us the notations \(r_{ij} = |\mathbf{r}_i - \mathbf{r}_j|\) and \(R_{i\alpha} = |\mathbf{r}_i - \mathbf{R}_\alpha|\).
\(J_{\text{eN}}\) contains electron-nucleus terms:
\[ J_{\text{eN}}(\mathbf{r},\mathbf{R}) = \sum_{i=1}^{N_\text{elec}} \sum_{\alpha=1}^{N_\text{nucl}} \frac{a_1\, f(R_{i\alpha})}{1+a_2\, f(R_{i\alpha})} + \sum_{p=2}^{N_\text{ord}^a} a_{p+1}\, [f(R_{i\alpha})]^p - J_{eN}^\infty \]
\(J_{\text{ee}}\) contains electron-electron terms: \[ J_{\text{ee}}(\mathbf{r}) = \sum_{i=1}^{N_\text{elec}} \sum_{j=1}^{i-1} \frac{b_1\, f(r_{ij})}{1+b_2\, f(r_{ij})} + \sum_{p=2}^{N_\text{ord}^b} a_{p+1}\, [f(r_{ij})]^p - J_{ee}^\infty \]
and \(J_{\text{eeN}}\) contains electron-electron-Nucleus terms:
\[ J_{\text{eeN}}(\mathbf{r},\mathbf{R}) = \sum_{\alpha=1}^{N_{\text{nucl}}} \sum_{i=1}^{N_{\text{elec}}} \sum_{j=1}^{i-1} \sum_{p=2}^{N_{\text{ord}}} \sum_{k=0}^{p-1} \sum_{l=0}^{p-k-2\delta_{k,0}} c_{lkp\alpha} \left[ g({r}_{ij}) \right]^k \left[ \left[ g({R}_{i\alpha}) \right]^l + \left[ g({R}_{j\alpha}) \right]^l \right] \left[ g({R}_{i\,\alpha}) \, g({R}_{j\alpha}) \right]^{(p-k-l)/2} \]
\(c_{lkp\alpha}\) are non-zero only when \(p-k-l\) is even.
\(f\) and \(g\) are scaling function defined as
\[ f(r) = \frac{1-e^{-\kappa\, r}}{\kappa} \text{ and } g(r) = e^{-\kappa\, r}. \]
The terms \(J_{\text{ee}}^\infty\) and \(J_{\text{eN}}^\infty\) are shifts to ensure that \(J_{\text{ee}}\) and \(J_{\text{eN}}\) have an asymptotic value of zero.
2 Context
The following data stored in the context:
| Variable | Type | In/Out | Description |
|---|---|---|---|
uninitialized |
int32_t |
in | Keeps bits set for uninitialized data |
aord_num |
int64_t |
in | The number of a coeffecients |
bord_num |
int64_t |
in | The number of b coeffecients |
cord_num |
int64_t |
in | The number of c coeffecients |
type_nucl_num |
int64_t |
in | Number of Nucleii types |
type_nucl_vector |
int64_t[nucl_num] |
in | IDs of types of Nucleii |
aord_vector |
double[aord_num + 1][type_nucl_num] |
in | Order of a polynomial coefficients |
bord_vector |
double[bord_num + 1] |
in | Order of b polynomial coefficients |
cord_vector |
double[cord_num][type_nucl_num] |
in | Order of c polynomial coefficients |
factor_ee |
double[walk_num] |
out | Jastrow factor: electron-electron part |
factor_ee_date |
uint64_t |
out | Jastrow factor: electron-electron part |
factor_en |
double[walk_num] |
out | Jastrow factor: electron-nucleus part |
factor_en_date |
uint64_t |
out | Jastrow factor: electron-nucleus part |
factor_een |
double[walk_num] |
out | Jastrow factor: electron-electron-nucleus part |
factor_een_date |
uint64_t |
out | Jastrow factor: electron-electron-nucleus part |
factor_ee_deriv_e |
double[4][nelec][walk_num] |
out | Derivative of the Jastrow factor: electron-electron-nucleus part |
factor_ee_deriv_e_date |
uint64_t |
out | Keep track of the date for the derivative |
factor_en_deriv_e |
double[4][nelec][walk_num] |
out | Derivative of the Jastrow factor: electron-electron-nucleus part |
factor_en_deriv_e_date |
uint64_t |
out | Keep track of the date for the en derivative |
factor_een_deriv_e |
double[4][nelec][walk_num] |
out | Derivative of the Jastrow factor: electron-electron-nucleus part |
factor_een_deriv_e_date |
uint64_t |
out | Keep track of the date for the een derivative |
computed data:
| Variable | Type | In/Out | Description |
|---|---|---|---|
dim_cord_vect |
int64_t |
Number of unique C coefficients | |
dim_cord_vect_date |
uint64_t |
Number of unique C coefficients | |
asymp_jasb |
double[2] |
Asymptotic component | |
asymp_jasb_date |
uint64_t |
Asymptotic component | |
cord_vect_full |
double[dim_cord_vect][nucl_num] |
vector of non-zero coefficients | |
cord_vect_full_date |
uint64_t |
Keep track of changes here | |
lkpm_combined_index |
int64_t[4][dim_cord_vect] |
Transform l,k,p, and m into consecutive indices | |
lkpm_combined_index_date |
uint64_t |
Transform l,k,p, and m into consecutive indices | |
tmp_c |
double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num] |
out | vector of non-zero coefficients |
dtmp_c |
double[walk_num][elec_num][4][nucl_num][0:cord_num][0:cord_num-1] |
vector of non-zero coefficients |
een_rescaled_n |
double[walk_num][elec_num][nucl_num][0:cord_num] |
The electron-electron rescaled distances raised to the powers defined by cord | |
een_rescaled_n_date |
uint64_t |
Keep track of the date of creation | |
een_rescaled_e_deriv_e |
double[walk_num][elec_num][4][elec_num][0:cord_num] |
The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons | |
een_rescaled_e_deriv_e_date |
uint64_t |
Keep track of the date of creation | |
een_rescaled_n_deriv_e |
double[walk_num][elec_num][4][nucl_num][0:cord_num] |
The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons | |
een_rescaled_n_deriv_e_date |
uint64_t |
Keep track of the date of creation |
2.1 Data structure
typedef struct qmckl_jastrow_struct{ int32_t uninitialized; int64_t aord_num; int64_t bord_num; int64_t cord_num; int64_t type_nucl_num; uint64_t asymp_jasb_date; uint64_t tmp_c_date; uint64_t dtmp_c_date; uint64_t factor_ee_date; uint64_t factor_en_date; uint64_t factor_een_date; uint64_t factor_ee_deriv_e_date; uint64_t factor_en_deriv_e_date; uint64_t factor_een_deriv_e_date; int64_t* type_nucl_vector; double * aord_vector; double * bord_vector; double * cord_vector; double * asymp_jasb; double * factor_ee; double * factor_en; double * factor_een; double * factor_ee_deriv_e; double * factor_en_deriv_e; double * factor_een_deriv_e; int64_t dim_cord_vect; uint64_t dim_cord_vect_date; double * cord_vect_full; uint64_t cord_vect_full_date; int64_t* lkpm_combined_index; uint64_t lkpm_combined_index_date; double * tmp_c; double * dtmp_c; double * een_rescaled_e; double * een_rescaled_n; uint64_t een_rescaled_e_date; uint64_t een_rescaled_n_date; double * een_rescaled_e_deriv_e; double * een_rescaled_n_deriv_e; uint64_t een_rescaled_e_deriv_e_date; uint64_t een_rescaled_n_deriv_e_date; bool provided; char * type; } qmckl_jastrow_struct;
The uninitialized integer contains one bit set to one for each
initialization function which has not been called. It becomes equal
to zero after all initialization functions have been called. The
struct is then initialized and provided == true.
Some values are initialized by default, and are not concerned by
this mechanism.
qmckl_exit_code qmckl_init_jastrow(qmckl_context context);
qmckl_exit_code qmckl_init_jastrow(qmckl_context context) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return false; } qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); ctx->jastrow.uninitialized = (1 << 5) - 1; /* Default values */ return QMCKL_SUCCESS; }
2.2 Access functions
Along with these core functions, calculation of the jastrow factor requires the following additional information to be set:
When all the data for the AOs have been provided, the following
function returns true.
bool qmckl_jastrow_provided (const qmckl_context context);
2.3 Initialization functions
To prepare for the Jastrow and its derivative, all the following functions need to be called.
qmckl_exit_code qmckl_set_jastrow_ord_num (qmckl_context context, const int64_t aord_num, const int64_t bord_num, const int64_t cord_num); qmckl_exit_code qmckl_set_jastrow_type_nucl_num (qmckl_context context, const int64_t type_nucl_num); qmckl_exit_code qmckl_set_jastrow_type_nucl_vector (qmckl_context context, const int64_t* type_nucl_vector, const int64_t nucl_num); qmckl_exit_code qmckl_set_jastrow_aord_vector (qmckl_context context, const double * aord_vector, const int64_t size_max); qmckl_exit_code qmckl_set_jastrow_bord_vector (qmckl_context context, const double * bord_vector, const int64_t size_max); qmckl_exit_code qmckl_set_jastrow_cord_vector (qmckl_context context, const double * cord_vector, const int64_t size_max);
When the required information is completely entered, other data structures are computed to accelerate the calculations. The intermediates factors are precontracted using BLAS LEVEL 3 operations for an optimal flop count.
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; double rescale_factor_kappa_ee = 1.0; double rescale_factor_kappa_en = 1.0; double nucl_rescale_factor_kappa = 1.0; double* elec_coord = &(n2_elec_coord[0][0][0]); const double* nucl_charge = n2_charge; int64_t nucl_num = n2_nucl_num; double* nucl_coord = &(n2_nucl_coord[0][0]); int64_t size_max; /* Provide Electron data */ qmckl_exit_code rc; assert(!qmckl_electron_provided(context)); int64_t n; rc = qmckl_get_electron_num (context, &n); assert(rc == QMCKL_NOT_PROVIDED); rc = qmckl_get_electron_up_num (context, &n); assert(rc == QMCKL_NOT_PROVIDED); rc = qmckl_get_electron_down_num (context, &n); assert(rc == QMCKL_NOT_PROVIDED); rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num); assert(rc == QMCKL_SUCCESS); assert(!qmckl_electron_provided(context)); rc = qmckl_get_electron_up_num (context, &n); assert(rc == QMCKL_SUCCESS); assert(n == elec_up_num); rc = qmckl_get_electron_down_num (context, &n); assert(rc == QMCKL_SUCCESS); assert(n == elec_dn_num); rc = qmckl_get_electron_num (context, &n); assert(rc == QMCKL_SUCCESS); assert(n == elec_num); double k_ee = 0.; double k_en = 0.; rc = qmckl_get_electron_rescale_factor_ee (context, &k_ee); assert(rc == QMCKL_SUCCESS); assert(k_ee == 1.0); rc = qmckl_get_electron_rescale_factor_en (context, &k_en); assert(rc == QMCKL_SUCCESS); assert(k_en == 1.0); rc = qmckl_set_electron_rescale_factor_en(context, rescale_factor_kappa_en); assert(rc == QMCKL_SUCCESS); rc = qmckl_set_electron_rescale_factor_ee(context, rescale_factor_kappa_ee); assert(rc == QMCKL_SUCCESS); rc = qmckl_get_electron_rescale_factor_ee (context, &k_ee); assert(rc == QMCKL_SUCCESS); assert(k_ee == rescale_factor_kappa_ee); rc = qmckl_get_electron_rescale_factor_en (context, &k_en); assert(rc == QMCKL_SUCCESS); assert(k_en == rescale_factor_kappa_en); int64_t w; rc = qmckl_get_electron_walk_num (context, &w); assert(rc == QMCKL_NOT_PROVIDED); rc = qmckl_set_electron_walk_num (context, walk_num); assert(rc == QMCKL_SUCCESS); rc = qmckl_get_electron_walk_num (context, &w); assert(rc == QMCKL_SUCCESS); assert(w == walk_num); assert(qmckl_electron_provided(context)); rc = qmckl_set_electron_coord (context, 'N', elec_coord, walk_num*3*elec_num); assert(rc == QMCKL_SUCCESS); double elec_coord2[walk_num*3*elec_num]; rc = qmckl_get_electron_coord (context, 'N', elec_coord2, walk_num*3*elec_num); assert(rc == QMCKL_SUCCESS); for (int64_t i=0 ; i<3*elec_num ; ++i) { assert( elec_coord[i] == elec_coord2[i] ); } /* Provide Nucleus data */ assert(!qmckl_nucleus_provided(context)); rc = qmckl_get_nucleus_num (context, &n); assert(rc == QMCKL_NOT_PROVIDED); rc = qmckl_set_nucleus_num (context, nucl_num); assert(rc == QMCKL_SUCCESS); assert(!qmckl_nucleus_provided(context)); rc = qmckl_get_nucleus_num (context, &n); assert(rc == QMCKL_SUCCESS); assert(n == nucl_num); double k; rc = qmckl_get_nucleus_rescale_factor (context, &k); assert(rc == QMCKL_SUCCESS); assert(k == 1.0); rc = qmckl_set_nucleus_rescale_factor (context, nucl_rescale_factor_kappa); assert(rc == QMCKL_SUCCESS); rc = qmckl_get_nucleus_rescale_factor (context, &k); assert(rc == QMCKL_SUCCESS); assert(k == nucl_rescale_factor_kappa); double nucl_coord2[3*nucl_num]; rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2, 3*nucl_num); assert(rc == QMCKL_NOT_PROVIDED); rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), 3*nucl_num); assert(rc == QMCKL_SUCCESS); assert(!qmckl_nucleus_provided(context)); rc = qmckl_get_nucleus_coord (context, 'N', nucl_coord2, nucl_num*3); assert(rc == QMCKL_SUCCESS); for (int64_t k=0 ; k<3 ; ++k) { for (int64_t i=0 ; i<nucl_num ; ++i) { assert( nucl_coord[nucl_num*k+i] == nucl_coord2[3*i+k] ); } } rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2, nucl_num*3); assert(rc == QMCKL_SUCCESS); for (int64_t i=0 ; i<3*nucl_num ; ++i) { assert( nucl_coord[i] == nucl_coord2[i] ); } double nucl_charge2[nucl_num]; rc = qmckl_get_nucleus_charge(context, nucl_charge2, nucl_num); assert(rc == QMCKL_NOT_PROVIDED); rc = qmckl_set_nucleus_charge(context, nucl_charge, nucl_num); assert(rc == QMCKL_SUCCESS); rc = qmckl_get_nucleus_charge(context, nucl_charge2, nucl_num); assert(rc == QMCKL_SUCCESS); for (int64_t i=0 ; i<nucl_num ; ++i) { assert( nucl_charge[i] == nucl_charge2[i] ); } assert(qmckl_nucleus_provided(context));
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 Asymptotic component for \(J_{ee}\)
Calculate the asymptotic component asymp_jasb to be substracted from the final
electron-electron jastrow factor \(J_{\text{ee}}\). The asymptotic component is calculated
via the bord_vector and the electron-electron rescale factor rescale_factor_kappa.
\[ J_{\text{ee}}^{\infty} = \frac{b_1 \kappa^{-1}}{1 + b_2 \kappa^{-1}} \]
3.1.1 Get
qmckl_exit_code qmckl_get_jastrow_asymp_jasb(qmckl_context context, double* const asymp_jasb, const int64_t size_max);
3.1.2 Compute
| Variable | Type | In/Out | Description |
|---|---|---|---|
context |
qmckl_context |
in | Global state |
bord_num |
int64_t |
in | Order of the polynomial |
bord_vector |
double[bord_num+1] |
in | Values of b |
rescale_factor_kappa_ee |
double |
in | Electron coordinates |
asymp_jasb |
double[2] |
out | Asymptotic value |
integer function qmckl_compute_asymp_jasb_f(context, bord_num, bord_vector, rescale_factor_kappa_ee, asymp_jasb) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: bord_num double precision , intent(in) :: bord_vector(bord_num + 1) double precision , intent(in) :: rescale_factor_kappa_ee double precision , intent(out) :: asymp_jasb(2) integer*8 :: i, p double precision :: kappa_inv, x, asym_one kappa_inv = 1.0d0 / rescale_factor_kappa_ee info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif if (bord_num <= 0) then info = QMCKL_INVALID_ARG_2 return endif asym_one = bord_vector(1) * kappa_inv / (1.0d0 + bord_vector(2) * kappa_inv) asymp_jasb(:) = (/asym_one, 0.5d0 * asym_one/) do i = 1, 2 x = kappa_inv do p = 2, bord_num x = x * kappa_inv asymp_jasb(i) = asymp_jasb(i) + bord_vector(p + 1) * x end do end do end function qmckl_compute_asymp_jasb_f
qmckl_exit_code qmckl_compute_asymp_jasb ( const qmckl_context context, const int64_t bord_num, const double* bord_vector, const double rescale_factor_kappa_ee, double* const asymp_jasb ) { double kappa_inv, x, asym_one; kappa_inv = 1.0 / rescale_factor_kappa_ee; if (context == QMCKL_NULL_CONTEXT){ return QMCKL_INVALID_CONTEXT; } if (bord_num <= 0) { return QMCKL_INVALID_ARG_2; } asym_one = bord_vector[0] * kappa_inv / (1.0 + bord_vector[1] * kappa_inv); asymp_jasb[0] = asym_one; asymp_jasb[1] = 0.5 * asym_one; for (int i = 0 ; i <= 1; ++i) { x = kappa_inv; for (int p = 1; p < bord_num; ++p){ x = x * kappa_inv; asymp_jasb[i] = asymp_jasb[i] + bord_vector[p + 1] * x; } } return QMCKL_SUCCESS; }
qmckl_exit_code qmckl_compute_asymp_jasb ( const qmckl_context context, const int64_t bord_num, const double* bord_vector, const double rescale_factor_kappa_ee, double* const asymp_jasb );
3.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* aord_vector = &(n2_aord_vector[0][0]); double* bord_vector = &(n2_bord_vector[0]); double* cord_vector = &(n2_cord_vector[0][0]); int64_t dim_cord_vect=0; /* Initialize the Jastrow data */ rc = qmckl_init_jastrow(context); assert(!qmckl_jastrow_provided(context)); /* Set the data */ rc = qmckl_set_jastrow_ord_num(context, aord_num, bord_num, cord_num); assert(rc == QMCKL_SUCCESS); rc = qmckl_set_jastrow_type_nucl_num(context, type_nucl_num); assert(rc == QMCKL_SUCCESS); rc = qmckl_set_jastrow_type_nucl_vector(context, type_nucl_vector, nucl_num); assert(rc == QMCKL_SUCCESS); rc = qmckl_set_jastrow_aord_vector(context, aord_vector,(aord_num+1)*type_nucl_num); assert(rc == QMCKL_SUCCESS); rc = qmckl_set_jastrow_bord_vector(context, bord_vector,(bord_num+1)); assert(rc == QMCKL_SUCCESS); rc = qmckl_get_jastrow_dim_cord_vect(context, &dim_cord_vect); assert(rc == QMCKL_SUCCESS); rc = qmckl_set_jastrow_cord_vector(context, cord_vector,dim_cord_vect*type_nucl_num); assert(rc == QMCKL_SUCCESS); /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_provided(context)); double asymp_jasb[2]; rc = qmckl_get_jastrow_asymp_jasb(context, asymp_jasb,2); // calculate asymp_jasb assert(fabs(asymp_jasb[0]-0.5323750557252571) < 1.e-12); assert(fabs(asymp_jasb[1]-0.31567342786262853) < 1.e-12);
3.2 Electron-electron component \(f_{ee}\)
Calculate the electron-electron jastrow component factor_ee using the asymp_jasb
componenet and the electron-electron rescaled distances ee_distance_rescaled.
\[
f_{ee} = \sum_{i,j
3.2.1 Get
qmckl_exit_code qmckl_get_jastrow_factor_ee(qmckl_context context, double* const factor_ee, const int64_t size_max);
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 |
up_num |
int64_t |
in | Number of alpha electrons |
bord_num |
int64_t |
in | Number of coefficients |
bord_vector |
double[bord_num+1] |
in | List of coefficients |
ee_distance_rescaled |
double[walk_num][elec_num][elec_num] |
in | Electron-electron distances |
asymp_jasb |
double[2] |
in | Electron-electron distances |
factor_ee |
double[walk_num] |
out | Electron-electron distances |
integer function qmckl_compute_factor_ee_f(context, walk_num, elec_num, up_num, bord_num, & bord_vector, ee_distance_rescaled, asymp_jasb, factor_ee) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: walk_num, elec_num, bord_num, up_num double precision , intent(in) :: bord_vector(bord_num + 1) double precision , intent(in) :: ee_distance_rescaled(elec_num, elec_num, walk_num) double precision , intent(in) :: asymp_jasb(2) double precision , intent(out) :: factor_ee(walk_num) integer*8 :: i, j, p, ipar, nw double precision :: pow_ser, x, spin_fact, power_ser info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif if (walk_num <= 0) then info = QMCKL_INVALID_ARG_2 return endif if (elec_num <= 0) then info = QMCKL_INVALID_ARG_3 return endif if (bord_num <= 0) then info = QMCKL_INVALID_ARG_4 return endif factor_ee = 0.0d0 do nw =1, walk_num do j = 1, elec_num do i = 1, j - 1 x = ee_distance_rescaled(i,j,nw) power_ser = 0.0d0 spin_fact = 1.0d0 ipar = 1 do p = 2, bord_num x = x * ee_distance_rescaled(i,j,nw) power_ser = power_ser + bord_vector(p + 1) * x end do if(j .LE. up_num .OR. i .GT. up_num) then spin_fact = 0.5d0 ipar = 2 endif factor_ee(nw) = factor_ee(nw) + spin_fact * bord_vector(1) * & ee_distance_rescaled(i,j,nw) / & (1.0d0 + bord_vector(2) * & ee_distance_rescaled(i,j,nw)) & -asymp_jasb(ipar) + power_ser end do end do end do end function qmckl_compute_factor_ee_f
qmckl_exit_code qmckl_compute_factor_ee ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t up_num, const int64_t bord_num, const double* bord_vector, const double* ee_distance_rescaled, const double* asymp_jasb, double* const factor_ee );
3.2.3 Test
/* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_provided(context)); double factor_ee[walk_num]; rc = qmckl_get_jastrow_factor_ee(context, factor_ee, walk_num); // calculate factor_ee assert(fabs(factor_ee[0]+4.282760865958113) < 1.e-12);
3.3 Electron-electron component derivative \(f'_{ee}\)
Calculate the derivative of the factor_ee using the ee_distance_rescaled and
the electron-electron rescaled distances derivatives ee_distance_rescaled_deriv_e.
There are four components, the gradient which has 3 components in the \(x, y, z\)
directions and the laplacian as the last component.
TODO: Add equation
3.3.1 Get
qmckl_exit_code qmckl_get_jastrow_factor_ee_deriv_e(qmckl_context context, double* const factor_ee_deriv_e, const int64_t size_max);
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 |
up_num |
int64_t |
in | Number of alpha electrons |
bord_num |
int64_t |
in | Number of coefficients |
bord_vector |
double[bord_num+1] |
in | List of coefficients |
ee_distance_rescaled |
double[walk_num][elec_num][elec_num] |
in | Electron-electron distances |
ee_distance_rescaled_deriv_e |
double[walk_num][4][elec_num][elec_num] |
in | Electron-electron distances |
asymp_jasb |
double[2] |
in | Electron-electron distances |
factor_ee_deriv_e |
double[walk_num][4][elec_num] |
out | Electron-electron distances |
integer function qmckl_compute_factor_ee_deriv_e_f(context, walk_num, elec_num, up_num, bord_num, & bord_vector, ee_distance_rescaled, ee_distance_rescaled_deriv_e, & asymp_jasb, factor_ee_deriv_e) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: walk_num, elec_num, bord_num, up_num double precision , intent(in) :: bord_vector(bord_num + 1) double precision , intent(in) :: ee_distance_rescaled(elec_num, elec_num,walk_num) double precision , intent(in) :: ee_distance_rescaled_deriv_e(4,elec_num, elec_num,walk_num) double precision , intent(in) :: asymp_jasb(2) double precision , intent(out) :: factor_ee_deriv_e(elec_num,4,walk_num) integer*8 :: i, j, p, ipar, nw, ii double precision :: x, spin_fact, y double precision :: den, invden, invden2, invden3, xinv double precision :: lap1, lap2, lap3, third double precision, dimension(3) :: pow_ser_g double precision, dimension(4) :: dx info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif if (walk_num <= 0) then info = QMCKL_INVALID_ARG_2 return endif if (elec_num <= 0) then info = QMCKL_INVALID_ARG_3 return endif if (bord_num <= 0) then info = QMCKL_INVALID_ARG_4 return endif factor_ee_deriv_e = 0.0d0 third = 1.0d0 / 3.0d0 do nw =1, walk_num do j = 1, elec_num do i = 1, elec_num x = ee_distance_rescaled(i,j,nw) if(abs(x) < 1.0d-18) cycle pow_ser_g = 0.0d0 spin_fact = 1.0d0 den = 1.0d0 + bord_vector(2) * x invden = 1.0d0 / den invden2 = invden * invden invden3 = invden2 * invden xinv = 1.0d0 / (x + 1.0d-18) ipar = 1 dx(1) = ee_distance_rescaled_deriv_e(1, i, j, nw) dx(2) = ee_distance_rescaled_deriv_e(2, i, j, nw) dx(3) = ee_distance_rescaled_deriv_e(3, i, j, nw) dx(4) = ee_distance_rescaled_deriv_e(4, i, j, nw) if((i .LE. up_num .AND. j .LE. up_num ) .OR. & (i .GT. up_num .AND. j .GT. up_num)) then spin_fact = 0.5d0 endif lap1 = 0.0d0 lap2 = 0.0d0 lap3 = 0.0d0 do ii = 1, 3 x = ee_distance_rescaled(i, j, nw) if(abs(x) < 1.0d-18) cycle do p = 2, bord_num y = p * bord_vector(p + 1) * x pow_ser_g(ii) = pow_ser_g(ii) + y * dx(ii) lap1 = lap1 + (p - 1) * y * xinv * dx(ii) * dx(ii) lap2 = lap2 + y x = x * ee_distance_rescaled(i, j, nw) end do lap3 = lap3 - 2.0d0 * bord_vector(2) * dx(ii) * dx(ii) factor_ee_deriv_e( j, ii, nw) = factor_ee_deriv_e( j, ii, nw) + spin_fact * bord_vector(1) * & dx(ii) * invden2 + pow_ser_g(ii) end do ii = 4 lap2 = lap2 * dx(ii) * third lap3 = lap3 + den * dx(ii) lap3 = lap3 * (spin_fact * bord_vector(1) * invden3) factor_ee_deriv_e( j, ii, nw) = factor_ee_deriv_e( j, ii, nw) + lap1 + lap2 + lap3 end do end do end do end function qmckl_compute_factor_ee_deriv_e_f
qmckl_exit_code qmckl_compute_factor_ee_deriv_e ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t up_num, const int64_t bord_num, const double* bord_vector, const double* ee_distance_rescaled, const double* ee_distance_rescaled_deriv_e, const double* asymp_jasb, double* const factor_ee_deriv_e );
3.3.3 Test
/* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_provided(context)); // calculate factor_ee_deriv_e double factor_ee_deriv_e[walk_num][4][elec_num]; rc = qmckl_get_jastrow_factor_ee_deriv_e(context, &(factor_ee_deriv_e[0][0][0]),walk_num*4*elec_num); // check factor_ee_deriv_e assert(fabs(factor_ee_deriv_e[0][0][0]-0.16364894652107934) < 1.e-12); assert(fabs(factor_ee_deriv_e[0][1][0]+0.6927548119830084 ) < 1.e-12); assert(fabs(factor_ee_deriv_e[0][2][0]-0.073267755223968 ) < 1.e-12); assert(fabs(factor_ee_deriv_e[0][3][0]-1.5111672803213185 ) < 1.e-12);
3.4 Electron-nucleus component \(f_{en}\)
Calculate the electron-electron jastrow component factor_en using the aord_vector
coeffecients and the electron-nucleus rescaled distances en_distance_rescaled.
\[
f_{en} = \sum_{i,j
3.4.1 Get
qmckl_exit_code qmckl_get_jastrow_factor_en(qmckl_context context, double* const factor_en, const int64_t size_max);
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 nucleii |
type_nucl_num |
int64_t |
in | Number of unique nuclei |
type_nucl_vector |
int64_t[nucl_num] |
in | IDs of unique nucleii |
aord_num |
int64_t |
in | Number of coefficients |
aord_vector |
double[aord_num+1][type_nucl_num] |
in | List of coefficients |
en_distance_rescaled |
double[walk_num][nucl_num][elec_num] |
in | Electron-nucleus distances |
factor_en |
double[walk_num] |
out | Electron-nucleus jastrow |
integer function qmckl_compute_factor_en_f(context, walk_num, elec_num, nucl_num, type_nucl_num, & type_nucl_vector, aord_num, aord_vector, & en_distance_rescaled, factor_en) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: walk_num, elec_num, aord_num, nucl_num, type_nucl_num integer*8 , intent(in) :: type_nucl_vector(nucl_num) double precision , intent(in) :: aord_vector(aord_num + 1, type_nucl_num) double precision , intent(in) :: en_distance_rescaled(elec_num, nucl_num, walk_num) double precision , intent(out) :: factor_en(walk_num) integer*8 :: i, a, p, ipar, nw double precision :: x, spin_fact, power_ser info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif if (walk_num <= 0) then info = QMCKL_INVALID_ARG_2 return endif if (elec_num <= 0) then info = QMCKL_INVALID_ARG_3 return endif if (nucl_num <= 0) then info = QMCKL_INVALID_ARG_4 return endif if (aord_num <= 0) then info = QMCKL_INVALID_ARG_7 return endif factor_en = 0.0d0 do nw =1, walk_num do a = 1, nucl_num do i = 1, elec_num x = en_distance_rescaled(i, a, nw) power_ser = 0.0d0 do p = 2, aord_num x = x * en_distance_rescaled(i, a, nw) power_ser = power_ser + aord_vector(p + 1, type_nucl_vector(a)) * x end do factor_en(nw) = factor_en(nw) + aord_vector(1, type_nucl_vector(a)) * & en_distance_rescaled(i, a, nw) / & (1.0d0 + aord_vector(2, type_nucl_vector(a)) * & en_distance_rescaled(i, a, nw)) & + power_ser end do end do end do end function qmckl_compute_factor_en_f
qmckl_exit_code qmckl_compute_factor_en ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t type_nucl_num, const int64_t* type_nucl_vector, const int64_t aord_num, const double* aord_vector, const double* en_distance_rescaled, double* const factor_en );
3.4.3 Test
/* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_provided(context)); double factor_en[walk_num]; rc = qmckl_get_jastrow_factor_en(context, factor_en,walk_num); // calculate factor_en assert(fabs(factor_en[0]+5.865822569188727) < 1.e-12);
3.5 Electron-nucleus component derivative \(f'_{en}\)
Calculate the electron-electron jastrow component factor_en_deriv_e derivative
with respect to the electron coordinates using the en_distance_rescaled and
en_distance_rescaled_deriv_e which are already calculated previously.
TODO: write equations.
3.5.1 Get
qmckl_exit_code qmckl_get_jastrow_factor_en_deriv_e(qmckl_context context, double* const factor_en_deriv_e, const int64_t size_max);
3.5.2 Compute
| Variable | Type | In/Out | Description |
|---|---|---|---|
context |
qmckl_context |
in | Global state |
walk_num |
int64_t |
in | Number of walkers |
elec_num |
int64_t |
in | Number of electrons |
nucl_num |
int64_t |
in | Number of nucleii |
type_nucl_num |
int64_t |
in | Number of unique nuclei |
type_nucl_vector |
int64_t[nucl_num] |
in | IDs of unique nucleii |
aord_num |
int64_t |
in | Number of coefficients |
aord_vector |
double[aord_num+1][type_nucl_num] |
in | List of coefficients |
en_distance_rescaled |
double[walk_num][nucl_num][elec_num] |
in | Electron-nucleus distances |
en_distance_rescaled_deriv_e |
double[walk_num][4][nucl_num][elec_num] |
in | Electron-nucleus distance derivatives |
factor_en_deriv_e |
double[walk_num][4][elec_num] |
out | Electron-nucleus jastrow |
integer function qmckl_compute_factor_en_deriv_e_f(context, walk_num, elec_num, nucl_num, type_nucl_num, & type_nucl_vector, aord_num, aord_vector, & en_distance_rescaled, en_distance_rescaled_deriv_e, factor_en_deriv_e) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: walk_num, elec_num, aord_num, nucl_num, type_nucl_num integer*8 , intent(in) :: type_nucl_vector(nucl_num) double precision , intent(in) :: aord_vector(aord_num + 1, type_nucl_num) double precision , intent(in) :: en_distance_rescaled(elec_num, nucl_num, walk_num) double precision , intent(in) :: en_distance_rescaled_deriv_e(4, elec_num, nucl_num, walk_num) double precision , intent(out) :: factor_en_deriv_e(elec_num,4,walk_num) integer*8 :: i, a, p, ipar, nw, ii double precision :: x, spin_fact, den, invden, invden2, invden3, xinv double precision :: y, lap1, lap2, lap3, third double precision, dimension(3) :: power_ser_g double precision, dimension(4) :: dx info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif if (walk_num <= 0) then info = QMCKL_INVALID_ARG_2 return endif if (elec_num <= 0) then info = QMCKL_INVALID_ARG_3 return endif if (nucl_num <= 0) then info = QMCKL_INVALID_ARG_4 return endif if (aord_num <= 0) then info = QMCKL_INVALID_ARG_7 return endif factor_en_deriv_e = 0.0d0 third = 1.0d0 / 3.0d0 do nw =1, walk_num do a = 1, nucl_num do i = 1, elec_num x = en_distance_rescaled(i,a,nw) if(abs(x) < 1.0d-18) continue power_ser_g = 0.0d0 den = 1.0d0 + aord_vector(2, type_nucl_vector(a)) * x invden = 1.0d0 / den invden2 = invden * invden invden3 = invden2 * invden xinv = 1.0d0 / x do ii = 1, 4 dx(ii) = en_distance_rescaled_deriv_e(ii,i,a,nw) end do lap1 = 0.0d0 lap2 = 0.0d0 lap3 = 0.0d0 do ii = 1, 3 x = en_distance_rescaled(i, a, nw) do p = 2, aord_num y = p * aord_vector(p + 1, type_nucl_vector(a)) * x power_ser_g(ii) = power_ser_g(ii) + y * dx(ii) lap1 = lap1 + (p - 1) * y * xinv * dx(ii) * dx(ii) lap2 = lap2 + y x = x * en_distance_rescaled(i, a, nw) end do lap3 = lap3 - 2.0d0 * aord_vector(2, type_nucl_vector(a)) * dx(ii) * dx(ii) factor_en_deriv_e(i, ii, nw) = factor_en_deriv_e(i, ii, nw) + aord_vector(1, type_nucl_vector(a)) & * dx(ii) * invden2 & + power_ser_g(ii) end do ii = 4 lap2 = lap2 * dx(ii) * third lap3 = lap3 + den * dx(ii) lap3 = lap3 * aord_vector(1, type_nucl_vector(a)) * invden3 factor_en_deriv_e(i, ii, nw) = factor_en_deriv_e(i, ii, nw) + lap1 + lap2 + lap3 end do end do end do end function qmckl_compute_factor_en_deriv_e_f
qmckl_exit_code qmckl_compute_factor_en_deriv_e ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t type_nucl_num, const int64_t* type_nucl_vector, const int64_t aord_num, const double* aord_vector, const double* en_distance_rescaled, const double* en_distance_rescaled_deriv_e, double* const factor_en_deriv_e );
3.5.3 Test
/* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_provided(context)); // calculate factor_en_deriv_e double factor_en_deriv_e[walk_num][4][elec_num]; rc = qmckl_get_jastrow_factor_en_deriv_e(context, &(factor_en_deriv_e[0][0][0]),walk_num*4*elec_num); // check factor_en_deriv_e assert(fabs(factor_en_deriv_e[0][0][0]-0.11609919541763383) < 1.e-12); assert(fabs(factor_en_deriv_e[0][1][0]+0.23301394780804574) < 1.e-12); assert(fabs(factor_en_deriv_e[0][2][0]-0.17548337641865783) < 1.e-12); assert(fabs(factor_en_deriv_e[0][3][0]+0.9667363412285741 ) < 1.e-12);
3.6 Electron-electron rescaled distances for each order
een_rescaled_e stores the table of the rescaled distances between all
pairs of electrons and raised to the power \(p\) defined by cord_num:
\[ C_{ij,p} = \left( 1 - \exp{-\kappa C_{ij}} \right)^p \]
where \(C_{ij}\) is the matrix of electron-electron distances.
3.6.1 Get
qmckl_exit_code qmckl_get_jastrow_een_rescaled_e(qmckl_context context, double* const distance_rescaled, const int64_t size_max);
3.6.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_kappa_ee |
double |
in | Factor to rescale ee distances |
ee_distance |
double[walk_num][elec_num][elec_num] |
in | Electron-electron distances |
een_rescaled_e |
double[walk_num][0:cord_num][elec_num][elec_num] |
out | Electron-electron rescaled distances |
integer function qmckl_compute_een_rescaled_e_f(context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, & ee_distance, een_rescaled_e) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: walk_num integer*8 , intent(in) :: elec_num integer*8 , intent(in) :: cord_num double precision , intent(in) :: rescale_factor_kappa_ee double precision , intent(in) :: ee_distance(elec_num,elec_num,walk_num) double precision , intent(out) :: een_rescaled_e(elec_num,elec_num,0:cord_num,walk_num) double precision,dimension(:,:),allocatable :: een_rescaled_e_ij double precision :: x integer*8 :: i, j, k, l, nw allocate(een_rescaled_e_ij(elec_num * (elec_num - 1) / 2, cord_num + 1)) info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif if (walk_num <= 0) then info = QMCKL_INVALID_ARG_2 return endif if (elec_num <= 0) then info = QMCKL_INVALID_ARG_3 return endif if (cord_num <= 0) then info = QMCKL_INVALID_ARG_4 return endif ! Prepare table of exponentiated distances raised to appropriate power een_rescaled_e = 0.0d0 do nw = 1, walk_num een_rescaled_e_ij = 0.0d0 een_rescaled_e_ij(:, 1) = 1.0d0 k = 0 do j = 1, elec_num do i = 1, j - 1 k = k + 1 een_rescaled_e_ij(k, 2) = dexp(-rescale_factor_kappa_ee * ee_distance(i, j, nw)) end do end do do l = 2, cord_num do k = 1, elec_num * (elec_num - 1)/2 een_rescaled_e_ij(k, l + 1) = een_rescaled_e_ij(k, l + 1 - 1) * een_rescaled_e_ij(k, 2) end do end do ! prepare the actual een table een_rescaled_e(:, :, 0, nw) = 1.0d0 do l = 1, cord_num k = 0 do j = 1, elec_num do i = 1, j - 1 k = k + 1 x = een_rescaled_e_ij(k, l + 1) een_rescaled_e(i, j, l, nw) = x een_rescaled_e(j, i, l, nw) = x end do end do end do do l = 0, cord_num do j = 1, elec_num een_rescaled_e(j, j, l, nw) = 0.0d0 end do end do end do end function qmckl_compute_een_rescaled_e_f
qmckl_exit_code qmckl_compute_een_rescaled_e ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t cord_num, const double rescale_factor_kappa_ee, const double* ee_distance, double* const een_rescaled_e );
3.6.3 Test
assert(qmckl_electron_provided(context)); double een_rescaled_e[walk_num][(cord_num + 1)][elec_num][elec_num]; rc = qmckl_get_jastrow_een_rescaled_e(context, &(een_rescaled_e[0][0][0][0]),elec_num*elec_num*(cord_num+1)*walk_num); // value of (0,2,1) assert(fabs(een_rescaled_e[0][1][0][2]-0.08084493981483197) < 1.e-12); assert(fabs(een_rescaled_e[0][1][0][3]-0.1066745707571846) < 1.e-12); assert(fabs(een_rescaled_e[0][1][0][4]-0.01754273169464735) < 1.e-12); assert(fabs(een_rescaled_e[0][2][1][3]-0.02214680362033448) < 1.e-12); assert(fabs(een_rescaled_e[0][2][1][4]-0.0005700154999202759) < 1.e-12); assert(fabs(een_rescaled_e[0][2][1][5]-0.3424402276009091) < 1.e-12);
3.7 Electron-electron rescaled distances for each order and derivatives
een_rescaled_e_deriv_e stores the table of the derivatives of the
rescaled distances between all pairs of electrons and raised to the
power \(p\) defined by cord_num. Here we take its derivatives
required for the een jastrow.
TODO: write formulae
3.7.1 Get
qmckl_exit_code qmckl_get_jastrow_een_rescaled_e_deriv_e(qmckl_context context, double* const distance_rescaled, const int64_t size_max);
3.7.2 Compute
| Variable | Type | In/Out | Description |
|---|---|---|---|
context |
qmckl_context |
in | Global state |
walk_num |
int64_t |
in | Number of walkers |
elec_num |
int64_t |
in | Number of electrons |
cord_num |
int64_t |
in | Order of polynomials |
rescale_factor_kappa_ee |
double |
in | Factor to rescale ee distances |
coord_new |
double[walk_num][3][elec_num] |
in | Electron coordinates |
ee_distance |
double[walk_num][elec_num][elec_num] |
in | Electron-electron distances |
een_rescaled_e |
double[walk_num][0:cord_num][elec_num][elec_num] |
in | Electron-electron distances |
een_rescaled_e_deriv_e |
double[walk_num][0:cord_num][elec_num][4][elec_num] |
out | Electron-electron rescaled distances |
integer function qmckl_compute_factor_een_rescaled_e_deriv_e_f(context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, & coord_new, ee_distance, een_rescaled_e, een_rescaled_e_deriv_e) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: walk_num integer*8 , intent(in) :: elec_num integer*8 , intent(in) :: cord_num double precision , intent(in) :: rescale_factor_kappa_ee double precision , intent(in) :: coord_new(elec_num,3,walk_num) double precision , intent(in) :: ee_distance(elec_num,elec_num,walk_num) double precision , intent(in) :: een_rescaled_e(elec_num,elec_num,0:cord_num,walk_num) double precision , intent(out) :: een_rescaled_e_deriv_e(elec_num,4,elec_num,0:cord_num,walk_num) double precision,dimension(:,:,:),allocatable :: elec_dist_deriv_e double precision :: x, rij_inv, kappa_l integer*8 :: i, j, k, l, nw, ii allocate(elec_dist_deriv_e(4,elec_num,elec_num)) info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif if (walk_num <= 0) then info = QMCKL_INVALID_ARG_2 return endif if (elec_num <= 0) then info = QMCKL_INVALID_ARG_3 return endif if (cord_num <= 0) then info = QMCKL_INVALID_ARG_4 return endif ! Prepare table of exponentiated distances raised to appropriate power een_rescaled_e_deriv_e = 0.0d0 do nw = 1, walk_num do j = 1, elec_num do i = 1, elec_num rij_inv = 1.0d0 / ee_distance(i, j, nw) do ii = 1, 3 elec_dist_deriv_e(ii, i, j) = (coord_new(i, ii, nw) - coord_new(j, ii, nw)) * rij_inv end do elec_dist_deriv_e(4, i, j) = 2.0d0 * rij_inv end do elec_dist_deriv_e(:, j, j) = 0.0d0 end do ! prepare the actual een table do l = 1, cord_num kappa_l = - dble(l) * rescale_factor_kappa_ee do j = 1, elec_num do i = 1, elec_num een_rescaled_e_deriv_e(i, 1, j, l, nw) = kappa_l * elec_dist_deriv_e(1, i, j) een_rescaled_e_deriv_e(i, 2, j, l, nw) = kappa_l * elec_dist_deriv_e(2, i, j) een_rescaled_e_deriv_e(i, 3, j, l, nw) = kappa_l * elec_dist_deriv_e(3, i, j) een_rescaled_e_deriv_e(i, 4, j, l, nw) = kappa_l * elec_dist_deriv_e(4, i, j) een_rescaled_e_deriv_e(i, 4, j, l, nw) = een_rescaled_e_deriv_e(i, 4, j, l, nw) & + een_rescaled_e_deriv_e(i, 1, j, l, nw) * een_rescaled_e_deriv_e(i, 1, j, l, nw) & + een_rescaled_e_deriv_e(i, 2, j, l, nw) * een_rescaled_e_deriv_e(i, 2, j, l, nw) & + een_rescaled_e_deriv_e(i, 3, j, l, nw) * een_rescaled_e_deriv_e(i, 3, j, l, nw) een_rescaled_e_deriv_e(i, 1, j, l, nw) = een_rescaled_e_deriv_e(i, 1, j, l, nw) * & een_rescaled_e(i, j, l, nw) een_rescaled_e_deriv_e(i, 3, j, l, nw) = een_rescaled_e_deriv_e(i, 2, j, l, nw) * & een_rescaled_e(i, j, l, nw) een_rescaled_e_deriv_e(i, 3, j, l, nw) = een_rescaled_e_deriv_e(i, 3, j, l, nw) * & een_rescaled_e(i, j, l, nw) een_rescaled_e_deriv_e(i, 4, j, l, nw) = een_rescaled_e_deriv_e(i, 4, j, l, nw) * & een_rescaled_e(i, j, l, nw) end do end do end do end do end function qmckl_compute_factor_een_rescaled_e_deriv_e_f
qmckl_exit_code qmckl_compute_een_rescaled_e_deriv_e ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t cord_num, const double rescale_factor_kappa_ee, const double* coord_new, const double* ee_distance, const double* een_rescaled_e, double* const een_rescaled_e_deriv_e );
3.7.3 Test
//assert(qmckl_electron_provided(context)); double een_rescaled_e_deriv_e[walk_num][(cord_num + 1)][elec_num][4][elec_num]; size_max=walk_num*(cord_num + 1)*elec_num*4*elec_num; rc = qmckl_get_jastrow_een_rescaled_e_deriv_e(context, &(een_rescaled_e_deriv_e[0][0][0][0][0]),size_max); // value of (0,0,0,2,1) assert(fabs(een_rescaled_e_deriv_e[0][1][0][0][2] + 0.05991352796887283 ) < 1.e-12); assert(fabs(een_rescaled_e_deriv_e[0][1][0][0][3] + 0.011714035071545248 ) < 1.e-12); assert(fabs(een_rescaled_e_deriv_e[0][1][0][0][4] + 0.00441398875758468 ) < 1.e-12); assert(fabs(een_rescaled_e_deriv_e[0][2][1][0][3] + 0.013553180060167595 ) < 1.e-12); assert(fabs(een_rescaled_e_deriv_e[0][2][1][0][4] + 0.00041342909359870457) < 1.e-12); assert(fabs(een_rescaled_e_deriv_e[0][2][1][0][5] + 0.5880599146214673 ) < 1.e-12);
3.8 Electron-nucleus rescaled distances for each order
een_rescaled_n stores the table of the rescaled distances between
electrons and nucleii raised to the power \(p\) defined by cord_num:
\[ C_{ia,p} = \left( 1 - \exp{-\kappa C_{ia}} \right)^p \]
where \(C_{ia}\) is the matrix of electron-nucleus distances.
3.8.1 Get
qmckl_exit_code qmckl_get_jastrow_een_rescaled_n(qmckl_context context, double* const distance_rescaled, const int64_t size_max);
3.8.2 Compute
| Variable | Type | In/Out | Description |
|---|---|---|---|
context |
qmckl_context |
in | Global state |
walk_num |
int64_t |
in | Number of walkers |
elec_num |
int64_t |
in | Number of electrons |
nucl_num |
int64_t |
in | Number of atoms |
cord_num |
int64_t |
in | Order of polynomials |
rescale_factor_kappa_en |
double |
in | Factor to rescale ee distances |
en_distance |
double[walk_num][elec_num][nucl_num] |
in | Electron-nucleus distances |
een_rescaled_n |
double[walk_num][0:cord_num][nucl_num][elec_num] |
out | Electron-nucleus rescaled distances |
integer function qmckl_compute_een_rescaled_n_f(context, walk_num, elec_num, nucl_num, cord_num, rescale_factor_kappa_en, & en_distance, een_rescaled_n) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: walk_num integer*8 , intent(in) :: elec_num integer*8 , intent(in) :: nucl_num integer*8 , intent(in) :: cord_num double precision , intent(in) :: rescale_factor_kappa_en double precision , intent(in) :: en_distance(elec_num,nucl_num,walk_num) double precision , intent(out) :: een_rescaled_n(elec_num,nucl_num,0:cord_num,walk_num) double precision :: x integer*8 :: i, a, k, l, nw info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif if (walk_num <= 0) then info = QMCKL_INVALID_ARG_2 return endif if (elec_num <= 0) then info = QMCKL_INVALID_ARG_3 return endif if (nucl_num <= 0) then info = QMCKL_INVALID_ARG_4 return endif if (cord_num <= 0) then info = QMCKL_INVALID_ARG_5 return endif ! Prepare table of exponentiated distances raised to appropriate power een_rescaled_n = 0.0d0 do nw = 1, walk_num ! prepare the actual een table een_rescaled_n(:, :, 0, nw) = 1.0d0 do a = 1, nucl_num do i = 1, elec_num een_rescaled_n(i, a, 1, nw) = dexp(-rescale_factor_kappa_en * en_distance(i, a, nw)) end do end do do l = 2, cord_num do a = 1, nucl_num do i = 1, elec_num een_rescaled_n(i, a, l, nw) = een_rescaled_n(i, a, l - 1, nw) * een_rescaled_n(i, a, 1, nw) end do end do end do end do end function qmckl_compute_een_rescaled_n_f
qmckl_exit_code qmckl_compute_een_rescaled_n ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, const double rescale_factor_kappa_en, const double* en_distance, double* const een_rescaled_n );
3.8.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_een_rescaled_n(context, &(een_rescaled_n[0][0][0][0]),size_max); // value of (0,2,1) assert(fabs(een_rescaled_n[0][1][0][2]-0.10612983920006765) < 1.e-12); assert(fabs(een_rescaled_n[0][1][0][3]-0.135652809635553) < 1.e-12); assert(fabs(een_rescaled_n[0][1][0][4]-0.023391817607642338) < 1.e-12); assert(fabs(een_rescaled_n[0][2][1][3]-0.880957224822116) < 1.e-12); assert(fabs(een_rescaled_n[0][2][1][4]-0.027185942659395074) < 1.e-12); assert(fabs(een_rescaled_n[0][2][1][5]-0.01343938025140174) < 1.e-12);
3.9 Electron-nucleus rescaled distances for each order and derivatives
een_rescaled_n_deriv_e stores the table of the rescaled distances between
electrons and nucleii raised to the power \(p\) defined by cord_num:
3.9.1 Get
qmckl_exit_code qmckl_get_jastrow_een_rescaled_n_deriv_e(qmckl_context context, double* const distance_rescaled, const int64_t size_max);
3.9.2 Compute
| Variable | Type | In/Out | Description |
|---|---|---|---|
context |
qmckl_context |
in | Global state |
walk_num |
int64_t |
in | Number of walkers |
elec_num |
int64_t |
in | Number of electrons |
nucl_num |
int64_t |
in | Number of atoms |
cord_num |
int64_t |
in | Order of polynomials |
rescale_factor_kappa_en |
double |
in | Factor to rescale ee distances |
coord_new |
double[walk_num][3][elec_num] |
in | Electron coordinates |
coord |
double[3][nucl_num] |
in | Nuclear coordinates |
en_distance |
double[walk_num][elec_num][nucl_num] |
in | Electron-nucleus distances |
een_rescaled_n |
double[walk_num][0:cord_num][nucl_num][elec_num] |
in | Electron-nucleus distances |
een_rescaled_n_deriv_e |
double[walk_num][0:cord_num][nucl_num][4][elec_num] |
out | Electron-nucleus rescaled distances |
integer function qmckl_compute_factor_een_rescaled_n_deriv_e_f(context, walk_num, elec_num, nucl_num, & cord_num, rescale_factor_kappa_en, & coord_new, coord, en_distance, een_rescaled_n, een_rescaled_n_deriv_e) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: walk_num integer*8 , intent(in) :: elec_num integer*8 , intent(in) :: nucl_num integer*8 , intent(in) :: cord_num double precision , intent(in) :: rescale_factor_kappa_en double precision , intent(in) :: coord_new(elec_num,3,walk_num) double precision , intent(in) :: coord(nucl_num,3) double precision , intent(in) :: en_distance(elec_num,nucl_num,walk_num) double precision , intent(in) :: een_rescaled_n(elec_num,nucl_num,0:cord_num,walk_num) double precision , intent(out) :: een_rescaled_n_deriv_e(elec_num,4,nucl_num,0:cord_num,walk_num) double precision,dimension(:,:,:),allocatable :: elnuc_dist_deriv_e double precision :: x, ria_inv, kappa_l integer*8 :: i, a, k, l, nw, ii allocate(elnuc_dist_deriv_e(4, elec_num, nucl_num)) info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif if (walk_num <= 0) then info = QMCKL_INVALID_ARG_2 return endif if (elec_num <= 0) then info = QMCKL_INVALID_ARG_3 return endif if (nucl_num <= 0) then info = QMCKL_INVALID_ARG_4 return endif if (cord_num <= 0) then info = QMCKL_INVALID_ARG_5 return endif ! Prepare table of exponentiated distances raised to appropriate power een_rescaled_n_deriv_e = 0.0d0 do nw = 1, walk_num ! prepare the actual een table do a = 1, nucl_num do i = 1, elec_num ria_inv = 1.0d0 / en_distance(i, a, nw) do ii = 1, 3 elnuc_dist_deriv_e(ii, i, a) = (coord_new(i, ii, nw) - coord(a, ii)) * ria_inv end do elnuc_dist_deriv_e(4, i, a) = 2.0d0 * ria_inv end do end do do l = 0, cord_num kappa_l = - dble(l) * rescale_factor_kappa_en do a = 1, nucl_num do i = 1, elec_num een_rescaled_n_deriv_e(i, 1, a, l, nw) = kappa_l * elnuc_dist_deriv_e(1, i, a) een_rescaled_n_deriv_e(i, 2, a, l, nw) = kappa_l * elnuc_dist_deriv_e(2, i, a) een_rescaled_n_deriv_e(i, 3, a, l, nw) = kappa_l * elnuc_dist_deriv_e(3, i, a) een_rescaled_n_deriv_e(i, 4, a, l, nw) = kappa_l * elnuc_dist_deriv_e(4, i, a) een_rescaled_n_deriv_e(i, 4, a, l, nw) = een_rescaled_n_deriv_e(i, 4, a, l, nw) & + een_rescaled_n_deriv_e(i, 1, a, l, nw) * een_rescaled_n_deriv_e(i, 1, a, l, nw) & + een_rescaled_n_deriv_e(i, 2, a, l, nw) * een_rescaled_n_deriv_e(i, 2, a, l, nw) & + een_rescaled_n_deriv_e(i, 3, a, l, nw) * een_rescaled_n_deriv_e(i, 3, a, l, nw) een_rescaled_n_deriv_e(i, 1, a, l, nw) = een_rescaled_n_deriv_e(i, 1, a, l, nw) * & een_rescaled_n(i, a, l, nw) een_rescaled_n_deriv_e(i, 2, a, l, nw) = een_rescaled_n_deriv_e(i, 2, a, l, nw) * & een_rescaled_n(i, a, l, nw) een_rescaled_n_deriv_e(i, 3, a, l, nw) = een_rescaled_n_deriv_e(i, 3, a, l, nw) * & een_rescaled_n(i, a, l, nw) een_rescaled_n_deriv_e(i, 4, a, l, nw) = een_rescaled_n_deriv_e(i, 4, a, l, nw) * & een_rescaled_n(i, a, l, nw) end do end do end do end do end function qmckl_compute_factor_een_rescaled_n_deriv_e_f
qmckl_exit_code qmckl_compute_factor_een_rescaled_n_deriv_e ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, const double rescale_factor_kappa_en, const double* coord_new, const double* coord, const double* en_distance, const double* een_rescaled_n, double* const een_rescaled_n_deriv_e );
3.9.3 Test
assert(qmckl_electron_provided(context)); double een_rescaled_n_deriv_e[walk_num][(cord_num + 1)][nucl_num][4][elec_num]; size_max=walk_num*(cord_num + 1)*nucl_num*4*elec_num; rc = qmckl_get_jastrow_een_rescaled_n_deriv_e(context, &(een_rescaled_n_deriv_e[0][0][0][0][0]),size_max); // value of (0,2,1) assert(fabs(een_rescaled_n_deriv_e[0][1][0][0][2]+0.07633444246999128 ) < 1.e-12); assert(fabs(een_rescaled_n_deriv_e[0][1][0][0][3]-0.00033282346259738276) < 1.e-12); assert(fabs(een_rescaled_n_deriv_e[0][1][0][0][4]+0.004775370547333061 ) < 1.e-12); assert(fabs(een_rescaled_n_deriv_e[0][2][1][0][3]-0.1362654644223866 ) < 1.e-12); assert(fabs(een_rescaled_n_deriv_e[0][2][1][0][4]+0.0231253431662794 ) < 1.e-12); assert(fabs(een_rescaled_n_deriv_e[0][2][1][0][5]-0.001593334817691633 ) < 1.e-12);
3.10 Prepare for electron-electron-nucleus Jastrow \(f_{een}\)
Prepare cord_vect_full and lkpm_combined_index tables required for the
calculation of the three-body jastrow factor_een and its derivative
factor_een_deriv_e.
3.10.1 Get
qmckl_exit_code qmckl_get_jastrow_dim_cord_vect(qmckl_context context, int64_t* const dim_cord_vect); qmckl_exit_code qmckl_get_jastrow_cord_vect_full(qmckl_context context, double* const cord_vect_full); qmckl_exit_code qmckl_get_jastrow_lkpm_combined_index(qmckl_context context, int64_t* const lkpm_combined_index); qmckl_exit_code qmckl_get_jastrow_tmp_c(qmckl_context context, double* const tmp_c); qmckl_exit_code qmckl_get_jastrow_dtmp_c(qmckl_context context, double* const dtmp_c);
3.10.2 Compute dimcordvect
| Variable | Type | In/Out | Description |
|---|---|---|---|
context |
qmckl_context |
in | Global state |
cord_num |
int64_t |
in | Order of polynomials |
dim_cord_vect |
int64_t |
out | dimension of cordvectfull table |
integer function qmckl_compute_dim_cord_vect_f(context, cord_num, dim_cord_vect) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: cord_num integer*8 , intent(out) :: dim_cord_vect double precision :: x integer*8 :: i, a, k, l, p, lmax info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif if (cord_num <= 0) then info = QMCKL_INVALID_ARG_2 return endif dim_cord_vect = 0 do p = 2, cord_num do k = p - 1, 0, -1 if (k .ne. 0) then lmax = p - k else lmax = p - k - 2 endif do l = lmax, 0, -1 if (iand(p - k - l, 1_8) == 1) cycle dim_cord_vect = dim_cord_vect + 1 end do end do end do end function qmckl_compute_dim_cord_vect_f
qmckl_exit_code qmckl_compute_dim_cord_vect ( const qmckl_context context, const int64_t cord_num, int64_t* const dim_cord_vect );
3.10.3 Compute cordvectfull
| Variable | Type | In/Out | Description |
|---|---|---|---|
context |
qmckl_context |
in | Global state |
nucl_num |
int64_t |
in | Number of atoms |
dim_cord_vect |
int64_t |
in | dimension of cord full table |
type_nucl_num |
int64_t |
in | dimension of cord full table |
type_nucl_vector |
int64_t[nucl_num] |
in | dimension of cord full table |
cord_vector |
double[dim_cord_vect][type_nucl_num] |
in | dimension of cord full table |
cord_vect_full |
double[dim_cord_vect][nucl_num] |
out | Full list of coefficients |
integer function qmckl_compute_cord_vect_full_f(context, nucl_num, dim_cord_vect, type_nucl_num, & type_nucl_vector, cord_vector, cord_vect_full) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: nucl_num integer*8 , intent(in) :: dim_cord_vect integer*8 , intent(in) :: type_nucl_num integer*8 , intent(in) :: type_nucl_vector(nucl_num) double precision , intent(in) :: cord_vector(type_nucl_num, dim_cord_vect) double precision , intent(out) :: cord_vect_full(nucl_num,dim_cord_vect) double precision :: x integer*8 :: i, a, k, l, nw info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif if (nucl_num <= 0) then info = QMCKL_INVALID_ARG_2 return endif if (type_nucl_num <= 0) then info = QMCKL_INVALID_ARG_4 return endif if (dim_cord_vect <= 0) then info = QMCKL_INVALID_ARG_5 return endif do a = 1, nucl_num cord_vect_full(a,1:dim_cord_vect) = cord_vector(type_nucl_vector(a),1:dim_cord_vect) end do end function qmckl_compute_cord_vect_full_f
qmckl_exit_code qmckl_compute_cord_vect_full ( const qmckl_context context, const int64_t nucl_num, const int64_t dim_cord_vect, const int64_t type_nucl_num, const int64_t* type_nucl_vector, const double* cord_vector, double* const cord_vect_full );
3.10.4 Compute lkpmcombinedindex
| Variable | Type | In/Out | Description |
|---|---|---|---|
context |
qmckl_context |
in | Global state |
cord_num |
int64_t |
in | Order of polynomials |
dim_cord_vect |
int64_t |
in | dimension of cord full table |
lpkm_combined_index |
int64_t[4][dim_cord_vect] |
out | Full list of combined indices |
integer function qmckl_compute_lkpm_combined_index_f(context, cord_num, dim_cord_vect, & lkpm_combined_index) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: cord_num integer*8 , intent(in) :: dim_cord_vect integer*8 , intent(out) :: lkpm_combined_index(dim_cord_vect, 4) double precision :: x integer*8 :: i, a, k, l, kk, p, lmax, m info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif if (cord_num <= 0) then info = QMCKL_INVALID_ARG_2 return endif if (dim_cord_vect <= 0) then info = QMCKL_INVALID_ARG_3 return endif kk = 0 do p = 2, cord_num do k = p - 1, 0, -1 if (k .ne. 0) then lmax = p - k else lmax = p - k - 2 end if do l = lmax, 0, -1 if (iand(p - k - l, 1_8) .eq. 1) cycle m = (p - k - l)/2 kk = kk + 1 lkpm_combined_index(kk, 1) = l lkpm_combined_index(kk, 2) = k lkpm_combined_index(kk, 3) = p lkpm_combined_index(kk, 4) = m end do end do end do end function qmckl_compute_lkpm_combined_index_f
qmckl_exit_code qmckl_compute_lkpm_combined_index ( const qmckl_context context, const int64_t cord_num, const int64_t dim_cord_vect, int64_t* const lpkm_combined_index );
3.10.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 nucleii |
walk_num |
int64_t |
in | Number of walkers |
een_rescaled_e |
double[walk_num][0:cord_num][elec_num][elec_num] |
in | Electron-electron rescaled factor |
een_rescaled_n |
double[walk_num][0:cord_num][nucl_num][elec_num] |
in | Electron-nucleus rescaled factor |
tmp_c |
double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num] |
out | vector of non-zero coefficients |
integer function qmckl_compute_tmp_c_f(context, cord_num, elec_num, nucl_num, & walk_num, een_rescaled_e, een_rescaled_n, tmp_c) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: cord_num integer*8 , intent(in) :: elec_num integer*8 , intent(in) :: nucl_num integer*8 , intent(in) :: walk_num double precision , intent(in) :: een_rescaled_e(elec_num, elec_num, 0:cord_num, walk_num) double precision , intent(in) :: een_rescaled_n(elec_num, nucl_num, 0:cord_num, walk_num) double precision , intent(out) :: tmp_c(elec_num, nucl_num,0:cord_num, 0:cord_num-1, walk_num) double precision :: x integer*8 :: i, j, a, l, kk, p, lmax, nw character :: TransA, TransB double precision :: alpha, beta integer*8 :: M, N, K, LDA, LDB, LDC TransA = 'N' TransB = 'N' alpha = 1.0d0 beta = 0.0d0 info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif if (cord_num <= 0) then info = QMCKL_INVALID_ARG_2 return endif if (elec_num <= 0) then info = QMCKL_INVALID_ARG_3 return endif if (nucl_num <= 0) then info = QMCKL_INVALID_ARG_4 return endif M = elec_num N = nucl_num*(cord_num + 1) K = elec_num LDA = size(een_rescaled_e,1) LDB = size(een_rescaled_n,1) LDC = size(tmp_c,1) do nw=1, walk_num do i=0, cord_num-1 info = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & een_rescaled_e(1,1,i,nw),LDA*1_8, & een_rescaled_n(1,1,0,nw),LDB*1_8, & beta, & tmp_c(1,1,0,i,nw),LDC) end do end do end function qmckl_compute_tmp_c_f
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 );
3.10.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 nucleii |
walk_num |
int64_t |
in | Number of walkers |
een_rescaled_e_deriv_e |
double[walk_num][0:cord_num][elec_num][4][elec_num] |
in | Electron-electron rescaled factor derivatives |
een_rescaled_n |
double[walk_num][0:cord_num][nucl_num][elec_num] |
in | Electron-nucleus rescaled factor |
dtmp_c |
double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num] |
out | vector of non-zero coefficients |
integer function qmckl_compute_dtmp_c_f(context, cord_num, elec_num, nucl_num, & walk_num, een_rescaled_e_deriv_e, een_rescaled_n, dtmp_c) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: cord_num integer*8 , intent(in) :: elec_num integer*8 , intent(in) :: nucl_num integer*8 , intent(in) :: walk_num double precision , intent(in) :: een_rescaled_e_deriv_e(elec_num, 4, elec_num, 0:cord_num, walk_num) double precision , intent(in) :: een_rescaled_n(elec_num, nucl_num, 0:cord_num, walk_num) double precision , intent(out) :: dtmp_c(elec_num, 4, nucl_num,0:cord_num, 0:cord_num-1, walk_num) double precision :: x integer*8 :: i, j, a, l, kk, p, lmax, nw, ii character :: TransA, TransB double precision :: alpha, beta integer*8 :: M, N, K, LDA, LDB, LDC TransA = 'N' TransB = 'N' alpha = 1.0d0 beta = 0.0d0 info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif if (cord_num <= 0) then info = QMCKL_INVALID_ARG_2 return endif if (elec_num <= 0) then info = QMCKL_INVALID_ARG_3 return endif if (nucl_num <= 0) then info = QMCKL_INVALID_ARG_4 return endif M = 4*elec_num N = nucl_num*(cord_num + 1) K = elec_num LDA = 4*size(een_rescaled_e_deriv_e,1) LDB = size(een_rescaled_n,1) LDC = 4*size(dtmp_c,1) do nw=1, walk_num do i=0, cord_num-1 info = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & een_rescaled_e_deriv_e(1,1,1,i,nw),LDA*1_8, & een_rescaled_n(1,1,0,nw),LDB*1_8, & beta, & dtmp_c(1,1,1,0,i,nw),LDC) end do end do end function qmckl_compute_dtmp_c_f
qmckl_exit_code qmckl_compute_dtmp_c ( const qmckl_context context, const int64_t cord_num, const int64_t elec_num, const int64_t nucl_num, const int64_t walk_num, const double* een_rescaled_e_deriv_e, const double* een_rescaled_n, double* const dtmp_c );
3.10.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_tmp_c(context, &(tmp_c[0][0][0][0][0])); double dtmp_c[walk_num][cord_num][cord_num+1][nucl_num][4][elec_num]; rc = qmckl_get_jastrow_dtmp_c(context, &(dtmp_c[0][0][0][0][0][0])); assert(fabs(tmp_c[0][0][1][0][0] - 2.7083473948352403) < 1e-12); assert(fabs(dtmp_c[0][1][0][0][0][0] - 0.237440520852232) < 1e-12);
3.11 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.11.1 Get
qmckl_exit_code qmckl_get_jastrow_factor_een(qmckl_context context, double* const factor_een, const int64_t size_max);
3.11.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 nucleii |
cord_num |
int64_t |
in | order of polynomials |
dim_cord_vect |
int64_t |
in | dimension of full coefficient vector |
cord_vect_full |
double[dim_cord_vect][nucl_num] |
in | full coefficient vector |
lkpm_combined_index |
int64_t[4][dim_cord_vect] |
in | combined indices |
een_rescaled_e |
double[walk_num][elec_num][elec_num][0:cord_num] |
in | Electron-nucleus rescaled |
een_rescaled_n |
double[walk_num][elec_num][nucl_num][0:cord_num] |
in | Electron-nucleus rescaled factor |
factor_een |
double[walk_num] |
out | Electron-nucleus jastrow |
integer function qmckl_compute_factor_een_naive_f(context, walk_num, elec_num, nucl_num, cord_num,& dim_cord_vect, cord_vect_full, lkpm_combined_index, & een_rescaled_e, een_rescaled_n, factor_een) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: walk_num, elec_num, cord_num, nucl_num, dim_cord_vect integer*8 , intent(in) :: lkpm_combined_index(dim_cord_vect,4) double precision , intent(in) :: cord_vect_full(nucl_num, dim_cord_vect) double precision , intent(in) :: een_rescaled_e(0:cord_num, elec_num, elec_num, walk_num) double precision , intent(in) :: een_rescaled_n(0:cord_num, nucl_num, elec_num, walk_num) double precision , intent(out) :: factor_een(walk_num) integer*8 :: i, a, j, l, k, p, m, n, nw double precision :: accu, accu2, cn info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif if (walk_num <= 0) then info = QMCKL_INVALID_ARG_2 return endif if (elec_num <= 0) then info = QMCKL_INVALID_ARG_3 return endif if (nucl_num <= 0) then info = QMCKL_INVALID_ARG_4 return endif if (cord_num <= 0) then info = QMCKL_INVALID_ARG_5 return endif factor_een = 0.0d0 do nw =1, walk_num do n = 1, dim_cord_vect l = lkpm_combined_index(n, 1) k = lkpm_combined_index(n, 2) p = lkpm_combined_index(n, 3) m = lkpm_combined_index(n, 4) do a = 1, nucl_num accu2 = 0.0d0 cn = cord_vect_full(a, n) do j = 1, elec_num accu = 0.0d0 do i = 1, elec_num accu = accu + een_rescaled_e(k,i,j,nw) * & een_rescaled_n(m,a,i,nw) !if(nw .eq. 1) then ! print *,l,k,p,m,j,i,een_rescaled_e(k,i,j,nw), een_rescaled_n(m,a,i,nw), accu !endif end do accu2 = accu2 + accu * een_rescaled_n(m + l,a,j,nw) !print *, l,m,nw,accu, accu2, een_rescaled_n(m + l, a, j, nw), cn, factor_een(nw) end do factor_een(nw) = factor_een(nw) + accu2 * cn end do end do end do end function qmckl_compute_factor_een_naive_f
qmckl_exit_code qmckl_compute_factor_een_naive ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, const int64_t dim_cord_vect, const double* cord_vect_full, const int64_t* lkpm_combined_index, const double* een_rescaled_e, const double* een_rescaled_n, double* const factor_een );
3.11.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 nucleii |
cord_num |
int64_t |
in | order of polynomials |
dim_cord_vect |
int64_t |
in | dimension of full coefficient vector |
cord_vect_full |
double[dim_cord_vect][nucl_num] |
in | full coefficient vector |
lkpm_combined_index |
int64_t[4][dim_cord_vect] |
in | combined indices |
tmp_c |
double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num] |
vector of non-zero coefficients | |
een_rescaled_n |
double[walk_num][0:cord_num][nucl_num][elec_num] |
in | Electron-nucleus rescaled factor |
factor_een |
double[walk_num] |
out | Electron-nucleus jastrow |
integer function qmckl_compute_factor_een_f(context, walk_num, elec_num, nucl_num, cord_num, & dim_cord_vect, cord_vect_full, lkpm_combined_index, & tmp_c, een_rescaled_n, factor_een) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: walk_num, elec_num, cord_num, nucl_num, dim_cord_vect integer*8 , intent(in) :: lkpm_combined_index(dim_cord_vect,4) double precision , intent(in) :: cord_vect_full(nucl_num, dim_cord_vect) double precision , intent(in) :: tmp_c(elec_num, nucl_num,0:cord_num, 0:cord_num-1, walk_num) double precision , intent(in) :: een_rescaled_n(elec_num, nucl_num, 0:cord_num, walk_num) double precision , intent(out) :: factor_een(walk_num) integer*8 :: i, a, j, l, k, p, m, n, nw double precision :: accu, accu2, cn info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif if (walk_num <= 0) then info = QMCKL_INVALID_ARG_2 return endif if (elec_num <= 0) then info = QMCKL_INVALID_ARG_3 return endif if (nucl_num <= 0) then info = QMCKL_INVALID_ARG_4 return endif if (cord_num <= 0) then info = QMCKL_INVALID_ARG_5 return endif factor_een = 0.0d0 do nw =1, walk_num do n = 1, dim_cord_vect l = lkpm_combined_index(n, 1) k = lkpm_combined_index(n, 2) p = lkpm_combined_index(n, 3) m = lkpm_combined_index(n, 4) do a = 1, nucl_num cn = cord_vect_full(a, n) if(cn == 0.d0) cycle accu = 0.0d0 do j = 1, elec_num accu = accu + een_rescaled_n(j,a,m,nw) * tmp_c(j,a,m+l,k,nw) end do factor_een(nw) = factor_een(nw) + accu * cn end do end do end do end function qmckl_compute_factor_een_f
qmckl_exit_code qmckl_compute_factor_een ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, const int64_t dim_cord_vect, const double* cord_vect_full, const int64_t* lkpm_combined_index, double* tmp_c, const double* een_rescaled_n, double* const factor_een );
3.11.4 Test
/* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_provided(context)); double factor_een[walk_num]; rc = qmckl_get_jastrow_factor_een(context, &(factor_een[0]),walk_num); assert(fabs(factor_een[0] + 0.37407972141304213) < 1e-12);
3.12 Electron-electron-nucleus Jastrow \(f_{een}\) derivative
Calculate the electron-electron-nuclear three-body jastrow component factor_een_deriv_e
using the above prepared tables.
TODO: write equations.
3.12.1 Get
qmckl_exit_code qmckl_get_jastrow_factor_een_deriv_e(qmckl_context context, double* const factor_een_deriv_e, const int64_t size_max);
3.12.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 nucleii |
cord_num |
int64_t |
in | order of polynomials |
dim_cord_vect |
int64_t |
in | dimension of full coefficient vector |
cord_vect_full |
double[dim_cord_vect][nucl_num] |
in | full coefficient vector |
lkpm_combined_index |
int64_t[4][dim_cord_vect] |
in | combined indices |
een_rescaled_e |
double[walk_num][elec_num][elec_num][0:cord_num] |
in | Electron-nucleus rescaled |
een_rescaled_n |
double[walk_num][elec_num][nucl_num][0:cord_num] |
in | Electron-nucleus rescaled factor |
een_rescaled_e_deriv_e |
double[walk_num][elec_num][4][elec_num][0:cord_num] |
in | Electron-nucleus rescaled |
een_rescaled_n_deriv_e |
double[walk_num][elec_num][4][nucl_num][0:cord_num] |
in | Electron-nucleus rescaled factor |
factor_een_deriv_e |
double[walk_num][4][elec_num] |
out | Electron-nucleus jastrow |
integer function qmckl_compute_factor_een_deriv_e_naive_f(context, walk_num, elec_num, nucl_num, cord_num, & dim_cord_vect, cord_vect_full, lkpm_combined_index, & een_rescaled_e, een_rescaled_n, & een_rescaled_e_deriv_e, een_rescaled_n_deriv_e, factor_een_deriv_e)& result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: walk_num, elec_num, cord_num, nucl_num, dim_cord_vect integer*8 , intent(in) :: lkpm_combined_index(dim_cord_vect, 4) double precision , intent(in) :: cord_vect_full(nucl_num, dim_cord_vect) double precision , intent(in) :: een_rescaled_e(0:cord_num, elec_num, elec_num, walk_num) double precision , intent(in) :: een_rescaled_n(0:cord_num, nucl_num, elec_num, walk_num) double precision , intent(in) :: een_rescaled_e_deriv_e(0:cord_num, elec_num, 4, elec_num, walk_num) double precision , intent(in) :: een_rescaled_n_deriv_e(0:cord_num, nucl_num, 4, elec_num, walk_num) double precision , intent(out) :: factor_een_deriv_e(elec_num, 4, walk_num) integer*8 :: i, a, j, l, k, p, m, n, nw double precision :: accu, accu2, cn double precision :: daccu(1:4), daccu2(1:4) info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif if (walk_num <= 0) then info = QMCKL_INVALID_ARG_2 return endif if (elec_num <= 0) then info = QMCKL_INVALID_ARG_3 return endif if (nucl_num <= 0) then info = QMCKL_INVALID_ARG_4 return endif if (cord_num <= 0) then info = QMCKL_INVALID_ARG_5 return endif factor_een_deriv_e = 0.0d0 do nw =1, walk_num do n = 1, dim_cord_vect l = lkpm_combined_index(n, 1) k = lkpm_combined_index(n, 2) p = lkpm_combined_index(n, 3) m = lkpm_combined_index(n, 4) do a = 1, nucl_num cn = cord_vect_full(a, n) do j = 1, elec_num accu = 0.0d0 accu2 = 0.0d0 daccu = 0.0d0 daccu2 = 0.0d0 do i = 1, elec_num accu = accu + een_rescaled_e(k, i, j, nw) * & een_rescaled_n(m, a, i, nw) accu2 = accu2 + een_rescaled_e(k, i, j, nw) * & een_rescaled_n(m + l, a, i, nw) daccu(1:4) = daccu(1:4) + een_rescaled_e_deriv_e(k, j, 1:4, i, nw) * & een_rescaled_n(m, a, i, nw) daccu2(1:4) = daccu2(1:4) + een_rescaled_e_deriv_e(k, j, 1:4, i, nw) * & een_rescaled_n(m + l, a, i, nw) end do factor_een_deriv_e(j, 1:4, nw) = factor_een_deriv_e(j, 1:4, nw) + & (accu * een_rescaled_n_deriv_e(m + l, a, 1:4, j, nw) & + daccu(1:4) * een_rescaled_n(m + l, a, j, nw) & + daccu2(1:4) * een_rescaled_n(m, a, j, nw) & + accu2 * een_rescaled_n_deriv_e(m, a, 1:4, j, nw)) * cn factor_een_deriv_e(j, 4, nw) = factor_een_deriv_e(j, 4, nw) + 2.0d0 * ( & daccu (1) * een_rescaled_n_deriv_e(m + l, a, 1, j, nw) + & daccu (2) * een_rescaled_n_deriv_e(m + l, a, 2, j, nw) + & daccu (3) * een_rescaled_n_deriv_e(m + l, a, 3, j, nw) + & daccu2(1) * een_rescaled_n_deriv_e(m, a, 1, j, nw ) + & daccu2(2) * een_rescaled_n_deriv_e(m, a, 2, j, nw ) + & daccu2(3) * een_rescaled_n_deriv_e(m, a, 3, j, nw ) ) * cn end do end do end do end do end function qmckl_compute_factor_een_deriv_e_naive_f
qmckl_exit_code qmckl_compute_factor_een_deriv_e_naive ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, const int64_t dim_cord_vect, const double* cord_vect_full, const int64_t* lkpm_combined_index, const double* een_rescaled_e, const double* een_rescaled_n, const double* een_rescaled_e_deriv_e, const double* een_rescaled_n_deriv_e, double* const factor_een_deriv_e );
3.12.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 nucleii |
cord_num |
int64_t |
in | order of polynomials |
dim_cord_vect |
int64_t |
in | dimension of full coefficient vector |
cord_vect_full |
double[dim_cord_vect][nucl_num] |
in | full coefficient vector |
lkpm_combined_index |
int64_t[4][dim_cord_vect] |
in | combined indices |
tmp_c |
double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num] |
in | Temporary intermediate tensor |
dtmp_c |
double[walk_num][0:cord_num-1][0:cord_num][nucl_num][4][elec_num] |
in | vector of non-zero coefficients |
een_rescaled_n |
double[walk_num][0:cord_num][nucl_num][elec_num] |
in | Electron-nucleus rescaled factor |
een_rescaled_n_deriv_e |
double[walk_num][0:cord_num][nucl_num][4][elec_num] |
in | Derivative of Electron-nucleus rescaled factor |
factor_een_deriv_e |
double[walk_num][4][elec_num] |
out | Derivative of Electron-nucleus jastrow |
integer function qmckl_compute_factor_een_deriv_e_f(context, walk_num, elec_num, nucl_num, cord_num, dim_cord_vect, & cord_vect_full, lkpm_combined_index, & tmp_c, dtmp_c, een_rescaled_n, een_rescaled_n_deriv_e, factor_een_deriv_e) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: walk_num, elec_num, cord_num, nucl_num, dim_cord_vect integer*8 , intent(in) :: lkpm_combined_index(dim_cord_vect,4) double precision , intent(in) :: cord_vect_full(nucl_num, dim_cord_vect) double precision , intent(in) :: tmp_c(elec_num, nucl_num,0:cord_num, 0:cord_num-1, walk_num) double precision , intent(in) :: dtmp_c(elec_num, 4, nucl_num,0:cord_num, 0:cord_num-1, walk_num) double precision , intent(in) :: een_rescaled_n(elec_num, nucl_num, 0:cord_num, walk_num) double precision , intent(in) :: een_rescaled_n_deriv_e(elec_num, 4, nucl_num, 0:cord_num, walk_num) double precision , intent(out) :: factor_een_deriv_e(elec_num,4,walk_num) integer*8 :: i, a, j, l, k, p, m, n, nw, ii double precision :: accu, accu2, cn info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif if (walk_num <= 0) then info = QMCKL_INVALID_ARG_2 return endif if (elec_num <= 0) then info = QMCKL_INVALID_ARG_3 return endif if (nucl_num <= 0) then info = QMCKL_INVALID_ARG_4 return endif if (cord_num <= 0) then info = QMCKL_INVALID_ARG_5 return endif factor_een_deriv_e = 0.0d0 do nw =1, walk_num do n = 1, dim_cord_vect l = lkpm_combined_index(n, 1) k = lkpm_combined_index(n, 2) p = lkpm_combined_index(n, 3) m = lkpm_combined_index(n, 4) do a = 1, nucl_num cn = cord_vect_full(a, n) if(cn == 0.d0) cycle do ii = 1, 4 do j = 1, elec_num factor_een_deriv_e(j,ii,nw) = factor_een_deriv_e(j,ii,nw) + (& tmp_c(j,a,m,k,nw) * een_rescaled_n_deriv_e(j,ii,a,m+l,nw) + & (dtmp_c(j,ii,a,m,k,nw)) * een_rescaled_n(j,a,m+l,nw) + & (dtmp_c(j,ii,a,m+l,k,nw)) * een_rescaled_n(j,a,m ,nw) + & tmp_c(j,a,m+l,k,nw) * een_rescaled_n_deriv_e(j,ii,a,m,nw) & ) * cn end do end do cn = cn + cn do j = 1, elec_num factor_een_deriv_e(j,4,nw) = factor_een_deriv_e(j,4,nw) + (& (dtmp_c(j,1,a,m ,k,nw)) * een_rescaled_n_deriv_e(j,1,a,m+l,nw) + & (dtmp_c(j,2,a,m ,k,nw)) * een_rescaled_n_deriv_e(j,2,a,m+l,nw) + & (dtmp_c(j,3,a,m ,k,nw)) * een_rescaled_n_deriv_e(j,3,a,m+l,nw) + & (dtmp_c(j,1,a,m+l,k,nw)) * een_rescaled_n_deriv_e(j,1,a,m ,nw) + & (dtmp_c(j,2,a,m+l,k,nw)) * een_rescaled_n_deriv_e(j,2,a,m ,nw) + & (dtmp_c(j,3,a,m+l,k,nw)) * een_rescaled_n_deriv_e(j,3,a,m ,nw) & ) * cn end do end do end do end do end function qmckl_compute_factor_een_deriv_e_f
qmckl_exit_code qmckl_compute_factor_een_deriv_e ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, const int64_t dim_cord_vect, const double* cord_vect_full, const int64_t* lkpm_combined_index, const double* tmp_c, const double* dtmp_c, const double* een_rescaled_n, const double* een_rescaled_n_deriv_e, double* const factor_een_deriv_e );
3.12.4 Test
/* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_provided(context)); double factor_een_deriv_e[4][walk_num][elec_num]; rc = qmckl_get_jastrow_factor_een_deriv_e(context, &(factor_een_deriv_e[0][0][0]),4*walk_num*elec_num); assert(fabs(factor_een_deriv_e[0][0][0] + 0.0005481671107226865) < 1e-12);