From 368604633d126b6771e62c0a2009c76588c1070a Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 9 Sep 2022 10:36:38 +0200 Subject: [PATCH] Moved rescale factors in Jastrow, removed kappa from names --- org/qmckl_context.org | 3 + org/qmckl_determinant.org | 70 +- org/qmckl_electron.org | 1315 +------------------ org/qmckl_error.org | 63 +- org/qmckl_jastrow.org | 2611 ++++++++++++++++++++++++++++--------- org/qmckl_nucleus.org | 320 +---- org/qmckl_trexio.org | 73 +- 7 files changed, 2134 insertions(+), 2321 deletions(-) diff --git a/org/qmckl_context.org b/org/qmckl_context.org index dfcd906..cb00776 100644 --- a/org/qmckl_context.org +++ b/org/qmckl_context.org @@ -292,6 +292,9 @@ qmckl_context qmckl_context_create() { rc = qmckl_init_determinant(context); assert (rc == QMCKL_SUCCESS); + + rc = qmckl_init_jastrow(context); + assert (rc == QMCKL_SUCCESS); } /* Allocate qmckl_memory_struct */ diff --git a/org/qmckl_determinant.org b/org/qmckl_determinant.org index ba46707..78d54f1 100644 --- a/org/qmckl_determinant.org +++ b/org/qmckl_determinant.org @@ -1111,21 +1111,21 @@ const double* nucl_charge = chbrclf_charge; const double* nucl_coord = &(chbrclf_nucl_coord[0][0]); rc = qmckl_set_electron_num (context, chbrclf_elec_up_num, chbrclf_elec_dn_num); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); + +rc = qmckl_set_electron_coord (context, 'N', chbrclf_walk_num, elec_coord, chbrclf_walk_num*chbrclf_elec_num*3); +qmckl_check(context, rc); assert(qmckl_electron_provided(context)); -rc = qmckl_set_electron_coord (context, 'N', chbrclf_walk_num, elec_coord, chbrclf_walk_num*chbrclf_elec_num*3); -assert(rc == QMCKL_SUCCESS); - rc = qmckl_set_nucleus_num (context, chbrclf_nucl_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), chbrclf_nucl_num*3); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); rc = qmckl_set_nucleus_charge(context, nucl_charge, chbrclf_nucl_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(qmckl_nucleus_provided(context)); @@ -1145,57 +1145,57 @@ const char typ = 'G'; assert(!qmckl_ao_basis_provided(context)); rc = qmckl_set_ao_basis_type (context, typ); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(!qmckl_ao_basis_provided(context)); rc = qmckl_set_ao_basis_shell_num (context, chbrclf_shell_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(!qmckl_ao_basis_provided(context)); rc = qmckl_set_ao_basis_prim_num (context, chbrclf_prim_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(!qmckl_ao_basis_provided(context)); rc = qmckl_set_ao_basis_nucleus_index (context, nucleus_index, chbrclf_nucl_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(!qmckl_ao_basis_provided(context)); rc = qmckl_set_ao_basis_nucleus_shell_num (context, nucleus_shell_num, chbrclf_nucl_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(!qmckl_ao_basis_provided(context)); rc = qmckl_set_ao_basis_shell_ang_mom (context, shell_ang_mom, chbrclf_shell_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(!qmckl_ao_basis_provided(context)); rc = qmckl_set_ao_basis_shell_factor (context, shell_factor, chbrclf_shell_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(!qmckl_ao_basis_provided(context)); rc = qmckl_set_ao_basis_shell_prim_num (context, shell_prim_num, chbrclf_shell_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(!qmckl_ao_basis_provided(context)); rc = qmckl_set_ao_basis_shell_prim_index (context, shell_prim_index, chbrclf_shell_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(!qmckl_ao_basis_provided(context)); rc = qmckl_set_ao_basis_exponent (context, exponent, chbrclf_prim_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(!qmckl_ao_basis_provided(context)); rc = qmckl_set_ao_basis_coefficient (context, coefficient, chbrclf_prim_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(!qmckl_ao_basis_provided(context)); rc = qmckl_set_ao_basis_prim_factor (context, prim_factor, chbrclf_prim_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); rc = qmckl_set_ao_basis_ao_num(context, chbrclf_ao_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); rc = qmckl_set_ao_basis_ao_factor (context, ao_factor, chbrclf_ao_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(qmckl_ao_basis_provided(context)); @@ -1203,22 +1203,22 @@ assert(qmckl_ao_basis_provided(context)); double ao_vgl[chbrclf_walk_num*chbrclf_elec_num][5][chbrclf_ao_num]; rc = qmckl_get_ao_basis_ao_vgl(context, &(ao_vgl[0][0][0]), (int64_t) 5*chbrclf_walk_num*chbrclf_elec_num*chbrclf_ao_num); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); /* Set up MO data */ rc = qmckl_set_mo_basis_mo_num(context, chbrclf_mo_num); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); const double * mo_coefficient = &(chbrclf_mo_coef[0]); rc = qmckl_set_mo_basis_coefficient(context, mo_coefficient); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(qmckl_mo_basis_provided(context)); double mo_vgl[chbrclf_walk_num*chbrclf_elec_num][5][chbrclf_mo_num]; rc = qmckl_get_mo_basis_mo_vgl(context, &(mo_vgl[0][0][0]), 5*chbrclf_walk_num*chbrclf_elec_num*chbrclf_mo_num); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); /* Set up determinant data */ @@ -1238,19 +1238,19 @@ for(k = 0; k < det_num_beta; ++k) mo_index_beta[k][i][j] = j + 1; rc = qmckl_set_determinant_type (context, typ); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); rc = qmckl_set_determinant_det_num_alpha (context, det_num_alpha); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); rc = qmckl_set_determinant_det_num_beta (context, det_num_beta); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); rc = qmckl_set_determinant_mo_index_alpha (context, &(mo_index_alpha[0][0][0])); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); rc = qmckl_set_determinant_mo_index_beta (context, &(mo_index_beta[0][0][0])); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); // Get slater-determinant @@ -1258,10 +1258,10 @@ double det_vgl_alpha[det_num_alpha][chbrclf_walk_num][5][chbrclf_elec_up_num][ch double det_vgl_beta[det_num_beta][chbrclf_walk_num][5][chbrclf_elec_dn_num][chbrclf_elec_dn_num]; rc = qmckl_get_det_vgl_alpha(context, &(det_vgl_alpha[0][0][0][0][0])); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); rc = qmckl_get_det_vgl_beta(context, &(det_vgl_beta[0][0][0][0][0])); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); #+end_src @@ -2018,10 +2018,10 @@ double det_inv_matrix_alpha[det_num_alpha][chbrclf_walk_num][chbrclf_elec_up_num double det_inv_matrix_beta[det_num_beta][chbrclf_walk_num][chbrclf_elec_dn_num][chbrclf_elec_dn_num]; rc = qmckl_get_det_inv_matrix_alpha(context, &(det_inv_matrix_alpha[0][0][0][0])); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); rc = qmckl_get_det_inv_matrix_beta(context, &(det_inv_matrix_beta[0][0][0][0])); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); #+end_src @@ -2034,7 +2034,7 @@ assert (rc == QMCKL_SUCCESS); *** Test #+begin_src c :tangle (eval c_test) rc = qmckl_context_destroy(context); - assert (rc == QMCKL_SUCCESS); + qmckl_check(context, rc); return 0; } diff --git a/org/qmckl_electron.org b/org/qmckl_electron.org index f3144ce..a7164d7 100644 --- a/org/qmckl_electron.org +++ b/org/qmckl_electron.org @@ -85,8 +85,6 @@ int main() { | ~num~ | ~int64_t~ | Total number of electrons | | ~up_num~ | ~int64_t~ | Number of up-spin electrons | | ~down_num~ | ~int64_t~ | Number of down-spin electrons | - | ~rescale_factor_kappa_ee~ | ~double~ | The distance scaling factor | - | ~rescale_factor_kappa_en~ | ~double[nucl_num]~ | The distance scaling factor | | ~provided~ | ~bool~ | If true, ~electron~ is valid | | ~walker~ | ~qmckl_point~ | Current set of walkers | | ~walker_old~ | ~qmckl_point~ | Previous set of walkers | @@ -99,18 +97,10 @@ int main() { | ~ee_distance_date~ | ~uint64_t~ | Last modification date of the electron-electron distances | | ~en_distance~ | ~double[walker.num][nucl_num][num]~ | Electron-nucleus distances | | ~en_distance_date~ | ~uint64_t~ | Last modification date of the electron-electron distances | - | ~ee_distance_rescaled~ | ~double[walker.num][num][num]~ | Electron-electron rescaled distances | - | ~ee_distance_rescaled_date~ | ~uint64_t~ | Last modification date of the electron-electron distances | - | ~ee_distance_rescaled_deriv_e~ | ~double[walker.num][4][num][num]~ | Electron-electron rescaled distances derivatives | - | ~ee_distance_rescaled_deriv_e_date~ | ~uint64_t~ | Last modification date of the electron-electron distance derivatives | | ~ee_potential~ | ~double[walker.num]~ | Electron-electron potential energy | | ~ee_potential_date~ | ~uint64_t~ | Last modification date of the electron-electron potential | | ~en_potential~ | ~double[walker.num]~ | Electron-nucleus potential energy | | ~en_potential_date~ | ~int64_t~ | Date when the electron-nucleus potential energy was computed | - | ~en_distance_rescaled~ | ~double[walker.num][nucl_num][num]~ | Electron-nucleus distances | - | ~en_distance_rescaled_date~ | ~uint64_t~ | Last modification date of the electron-electron distances | - | ~en_distance_rescaled_deriv_e~ | ~double[walker.num][4][nucl_num][num]~ | Electron-electron rescaled distances derivatives | - | ~en_distance_rescaled_deriv_e_date~ | ~uint64_t~ | Last modification date of the electron-electron distance derivatives | ** Data structure @@ -127,24 +117,14 @@ typedef struct qmckl_electron_struct { int64_t down_num; qmckl_walker walker; qmckl_walker walker_old; - double rescale_factor_kappa_ee; - double* rescale_factor_kappa_en; - uint64_t ee_distance_date; - uint64_t en_distance_date; - uint64_t ee_potential_date; - uint64_t en_potential_date; - uint64_t ee_distance_rescaled_date; - uint64_t ee_distance_rescaled_deriv_e_date; - uint64_t en_distance_rescaled_date; - uint64_t en_distance_rescaled_deriv_e_date; + uint64_t ee_distance_date; + uint64_t en_distance_date; + uint64_t ee_potential_date; + uint64_t en_potential_date; double* ee_distance; double* en_distance; double* ee_potential; double* en_potential; - double* ee_distance_rescaled; - double* ee_distance_rescaled_deriv_e; - double* en_distance_rescaled; - double* en_distance_rescaled_deriv_e; int32_t uninitialized; bool provided; } qmckl_electron_struct; @@ -172,10 +152,7 @@ qmckl_exit_code qmckl_init_electron(qmckl_context context) { qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); - ctx->electron.uninitialized = (1 << 2) - 1; - - /* Default values */ - ctx->electron.rescale_factor_kappa_ee = 1.0; + ctx->electron.uninitialized = (1 << 1) - 1; return QMCKL_SUCCESS; } @@ -210,9 +187,6 @@ bool qmckl_electron_provided(const qmckl_context context) { #+begin_src c :comments org :tangle (eval h_func) qmckl_exit_code qmckl_set_electron_num (qmckl_context context, const int64_t up_num, const int64_t down_num); qmckl_exit_code qmckl_set_electron_coord (qmckl_context context, const char transp, const int64_t walk_num, const double* coord, const int64_t size_max); - -qmckl_exit_code qmckl_set_electron_rescale_factor_ee (qmckl_context context, const double kappa_ee); -qmckl_exit_code qmckl_set_electron_rescale_factor_en (qmckl_context context, const double* kappa_en, const int size_max); #+end_src #+NAME:pre2 @@ -274,80 +248,6 @@ qmckl_set_electron_num(qmckl_context context, #+end_src - Next we set the rescale parameter for the rescaled distance metric. - - #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -qmckl_exit_code -qmckl_set_electron_rescale_factor_ee(qmckl_context context, - const double rescale_factor_kappa_ee) { - - int32_t mask = 0; // can be changed - - <> - - if (rescale_factor_kappa_ee <= 0.0) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_2, - "qmckl_set_electron_rescale_factor_ee", - "rescale_factor_kappa_ee <= 0.0"); - } - - ctx->electron.rescale_factor_kappa_ee = rescale_factor_kappa_ee; - - <> -} - - #+end_src - - #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -qmckl_exit_code -qmckl_set_electron_rescale_factor_en(qmckl_context context, - const double* rescale_factor_kappa_en, const int64_t size_max) { - - int32_t mask = 1 << 1; - - <> - - if (rescale_factor_kappa_en == NULL) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_2, - "qmckl_set_electron_rescale_factor_en", - "Null pointer"); - } - - if (size_max < ctx->nucleus.num) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_3, - "qmckl_set_electron_rescale_factor_en", - "Array too small"); - } - - - if (ctx->electron.rescale_factor_kappa_en != NULL) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_3, - "qmckl_set_electron_rescale_factor_en", - "Already set"); - } - - qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; - mem_info.size = ctx->nucleus.num * sizeof(double); - ctx->electron.rescale_factor_kappa_en = (double*) qmckl_malloc(context, mem_info); - - for (int64_t i=0 ; inucleus.num ; ++i) { - if (rescale_factor_kappa_en[i] <= 0.0) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_2, - "qmckl_set_electron_rescale_factor_en", - "rescale_factor_kappa_en <= 0.0"); - } - ctx->electron.rescale_factor_kappa_en[i] = rescale_factor_kappa_en[i]; - } - - <> -} - #+end_src - #+begin_src f90 :comments org :tangle (eval fh_func) :noweb yes interface integer(c_int32_t) function qmckl_set_electron_num(context, alpha, beta) bind(C) @@ -384,7 +284,7 @@ qmckl_set_electron_coord(qmckl_context context, const int64_t size_max) { - int32_t mask = 0; // coord can be changed + int32_t mask = 0; <> @@ -588,74 +488,6 @@ qmckl_get_electron_walk_num (const qmckl_context context, int64_t* const walk_nu } #+end_src -*** Scaling factors Kappa - - #+begin_src c :comments org :tangle (eval h_func) :exports none -qmckl_exit_code qmckl_get_electron_rescale_factor_ee (const qmckl_context context, double* const rescale_factor_kappa_ee); -qmckl_exit_code qmckl_get_electron_rescale_factor_en (const qmckl_context context, double* const rescale_factor_kappa_en, const int64_t size_max); - #+end_src - - #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -qmckl_exit_code -qmckl_get_electron_rescale_factor_ee (const qmckl_context context, double* const rescale_factor_kappa_ee) { - - if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { - return QMCKL_INVALID_CONTEXT; - } - - if (rescale_factor_kappa_ee == NULL) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_2, - "qmckl_get_electron_rescale_factor_ee", - "rescale_factor_kappa_ee is a null pointer"); - } - - qmckl_context_struct* const ctx = (qmckl_context_struct*) context; - assert (ctx != NULL); - - assert (ctx->electron.rescale_factor_kappa_ee > 0.0); - - ,*rescale_factor_kappa_ee = ctx->electron.rescale_factor_kappa_ee; - return QMCKL_SUCCESS; -} - - -qmckl_exit_code -qmckl_get_electron_rescale_factor_en (const qmckl_context context, double* const rescale_factor_kappa_en, const int64_t size_max) { - if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { - return QMCKL_INVALID_CONTEXT; - } - - if (rescale_factor_kappa_en == NULL) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_2, - "qmckl_get_electron_rescale_factor_en", - "rescale_factor_kappa_en is a null pointer"); - } - - if (size_max < nucl_num) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_3, - "qmckl_get_electron_rescale_factor_en", - "Array to small"); - } - - qmckl_context_struct* const ctx = (qmckl_context_struct*) context; - assert (ctx != NULL); - - if (ctx->electron.rescale_factor_kappa_en != NULL) { - for (int64_t i=0 ; inucleus.num ; ++i) { - rescale_factor_kappa_en[i] = ctx->electron.rescale_factor_kappa_en[i]; - } - } else { - for (int64_t i=0 ; inucleus.num ; ++i) { - rescale_factor_kappa_en[i] = 1.0; - } - } - - return QMCKL_SUCCESS; -} - #+end_src *** Electron coordinates @@ -747,9 +579,6 @@ int64_t nucl_num = chbrclf_nucl_num; double* charge = chbrclf_charge; double* nucl_coord = &(chbrclf_nucl_coord[0][0]); -double rescale_factor_kappa_ee = 1.0; -double rescale_factor_kappa_en[nucl_num] = { 1.0, 1.0, 1.0, 1.0, 1.0} ; -double nucl_rescale_factor_kappa = 1.0; double* elec_coord = &(chbrclf_elec_coord[0][0][0]); @@ -771,69 +600,41 @@ assert(rc == QMCKL_NOT_PROVIDED); rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(qmckl_electron_provided(context)); rc = qmckl_get_electron_up_num (context, &n); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(n == elec_up_num); rc = qmckl_get_electron_down_num (context, &n); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(n == elec_dn_num); rc = qmckl_get_electron_num (context, &n); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(n == elec_num); -double k_ee = 0.; -double k_en[nucl_num] = { 0., 0., 0., 0., 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[0]), nucl_num); -assert(rc == QMCKL_SUCCESS); -for (int i=0 ; ielectron.num * ctx->electron.num * ctx->electron.walker.num; - memcpy(distance_rescaled, ctx->electron.ee_distance_rescaled, sze * sizeof(double)); - - return QMCKL_SUCCESS; -} - #+end_src - -*** Provide :noexport: - - #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none -qmckl_exit_code qmckl_provide_ee_distance_rescaled(qmckl_context context); - #+end_src - - #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -qmckl_exit_code qmckl_provide_ee_distance_rescaled(qmckl_context context) -{ - - if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { - return QMCKL_NULL_CONTEXT; - } - - qmckl_context_struct* const ctx = (qmckl_context_struct*) context; - assert (ctx != NULL); - - - /* Compute if necessary */ - if (ctx->electron.walker.point.date > ctx->electron.ee_distance_rescaled_date) { - - if (ctx->electron.walker.num > ctx->electron.walker_old.num) { - free(ctx->electron.ee_distance_rescaled); - ctx->electron.ee_distance_rescaled = NULL; - } - - /* Allocate array */ - if (ctx->electron.ee_distance_rescaled == NULL) { - - qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; - mem_info.size = ctx->electron.num * ctx->electron.num * - ctx->electron.walker.num * sizeof(double); - double* ee_distance_rescaled = (double*) qmckl_malloc(context, mem_info); - - if (ee_distance_rescaled == NULL) { - return qmckl_failwith( context, - QMCKL_ALLOCATION_FAILED, - "qmckl_ee_distance_rescaled", - NULL); - } - ctx->electron.ee_distance_rescaled = ee_distance_rescaled; - } - - qmckl_exit_code rc = - qmckl_compute_ee_distance_rescaled(context, - ctx->electron.num, - ctx->electron.rescale_factor_kappa_ee, - ctx->electron.walker.num, - ctx->electron.walker.point.coord.data, - ctx->electron.ee_distance_rescaled); - if (rc != QMCKL_SUCCESS) { - return rc; - } - - ctx->electron.ee_distance_rescaled_date = ctx->date; - } - - return QMCKL_SUCCESS; -} - #+end_src - -*** Compute - :PROPERTIES: - :Name: qmckl_compute_ee_distance_rescaled - :CRetType: qmckl_exit_code - :FRetType: qmckl_exit_code - :END: - - #+NAME: qmckl_ee_distance_rescaled_args - | Variable | Type | In/Out | Description | - |---------------------------+----------------------------------------+--------+--------------------------------------| - | ~context~ | ~qmckl_context~ | in | Global state | - | ~elec_num~ | ~int64_t~ | in | Number of electrons | - | ~rescale_factor_kappa_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 | - - #+begin_src f90 :comments org :tangle (eval f) :noweb yes -integer function qmckl_compute_ee_distance_rescaled_f(context, elec_num, rescale_factor_kappa_ee, walk_num, & - coord, ee_distance_rescaled) & - result(info) - use qmckl - implicit none - integer(qmckl_context), intent(in) :: context - integer*8 , intent(in) :: elec_num - double precision , intent(in) :: rescale_factor_kappa_ee - integer*8 , intent(in) :: walk_num - double precision , intent(in) :: coord(elec_num,walk_num,3) - double precision , intent(out) :: ee_distance_rescaled(elec_num,elec_num,walk_num) - - integer*8 :: k - - info = QMCKL_SUCCESS - - if (context == QMCKL_NULL_CONTEXT) then - info = QMCKL_INVALID_CONTEXT - return - endif - - if (elec_num <= 0) then - info = QMCKL_INVALID_ARG_2 - return - endif - - if (walk_num <= 0) then - info = QMCKL_INVALID_ARG_3 - return - endif - - do k=1,walk_num - info = qmckl_distance_rescaled(context, 'T', 'T', elec_num, elec_num, & - coord(1,k,1), elec_num * walk_num, & - coord(1,k,1), elec_num * walk_num, & - ee_distance_rescaled(1,1,k), elec_num, rescale_factor_kappa_ee) - if (info /= QMCKL_SUCCESS) then - exit - endif - end do - -end function qmckl_compute_ee_distance_rescaled_f - #+end_src - - #+begin_src c :tangle (eval h_private_func) :comments org :exports none -qmckl_exit_code qmckl_compute_ee_distance_rescaled ( - const qmckl_context context, - const int64_t elec_num, - const double rescale_factor_kappa_ee, - const int64_t walk_num, - const double* coord, - double* const ee_distance_rescaled ); - #+end_src - - #+CALL: generate_c_interface(table=qmckl_ee_distance_rescaled_args,rettyp=get_value("CRetType"),fname=get_value("Name")) - - #+RESULTS: - #+begin_src f90 :tangle (eval f) :comments org :exports none - integer(c_int32_t) function qmckl_compute_ee_distance_rescaled & - (context, elec_num, rescale_factor_kappa_ee, walk_num, coord, ee_distance_rescaled) & - bind(C) result(info) - - use, intrinsic :: iso_c_binding - implicit none - - integer (c_int64_t) , intent(in) , value :: context - integer (c_int64_t) , intent(in) , value :: elec_num - real (c_double ) , intent(in) , value :: rescale_factor_kappa_ee - integer (c_int64_t) , intent(in) , value :: walk_num - real (c_double ) , intent(in) :: coord(elec_num,3,walk_num) - real (c_double ) , intent(out) :: ee_distance_rescaled(elec_num,elec_num,walk_num) - - integer(c_int32_t), external :: qmckl_compute_ee_distance_rescaled_f - info = qmckl_compute_ee_distance_rescaled_f & - (context, elec_num, rescale_factor_kappa_ee, walk_num, coord, ee_distance_rescaled) - - end function qmckl_compute_ee_distance_rescaled - #+end_src - -*** Test - - #+begin_src python :results output :exports none -import numpy as np - -kappa = 1.0 - -elec_1_w1 = np.array( [ -2.26995253563, -5.15737533569, -2.22940072417 ]) -elec_2_w1 = np.array( [ 3.51983380318, -1.08717381954, -1.19617708027 ]) -elec_1_w2 = np.array( [ -2.34410619736, -3.20016115904, -1.53496759012 ]) -elec_2_w2 = np.array( [ 3.17996025085, -1.40260577202, 1.49473607540 ]) - -print ( "[0][0][0] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_1_w1-elec_1_w1)) )/kappa ) -print ( "[0][1][0] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_1_w1-elec_2_w1)) )/kappa ) -print ( "[1][0][0] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_2_w1-elec_1_w1)) )/kappa ) -print ( "[0][0][1] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_1_w2-elec_1_w2)) )/kappa ) -print ( "[0][1][1] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_1_w2-elec_2_w2)) )/kappa ) -print ( "[1][0][1] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_2_w2-elec_1_w2)) )/kappa ) - #+end_src - - #+RESULTS: - : [0][0][0] : 0.0 - : [0][1][0] : 0.9992169566605263 - : [1][0][0] : 0.9992169566605263 - : [0][0][1] : 0.0 - : [0][1][1] : 0.9985724058042633 - : [1][0][1] : 0.9985724058042633 - - #+begin_src c :tangle (eval c_test) -assert(qmckl_electron_provided(context)); - - -double ee_distance_rescaled[walk_num * elec_num * elec_num]; -rc = qmckl_get_electron_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.9992169566605263) < 1.e-12); - -// (0,0,1) == 0. -assert(ee_distance_rescaled[elec_num*elec_num] == 0.); - -// (1,0,1) == (0,1,1) -assert(ee_distance_rescaled[elec_num*elec_num+1] == ee_distance_rescaled[elec_num*elec_num+elec_num]); - -// value of (1,0,1) -assert(fabs(ee_distance_rescaled[elec_num*elec_num+1]-0.9985724058042633) < 1.e-12); - - #+end_src - -** Electron-electron rescaled distance gradients and Laplacian with respect to electron coords - - The rescaled distances which is given as $R = (1 - \exp{-\kappa r})/\kappa$ - needs to be perturbed with respect to the electorn coordinates. - This data is stored in the ~ee_distance_rescaled_deriv_e~ tensor. The - The first three elements of this three index tensor ~[4][num][num]~ gives the - derivatives in the x, y, and z directions $dx, dy, dz$ and the last index - gives the Laplacian $\partial x^2 + \partial y^2 + \partial z^2$. - -*** Get - - #+begin_src c :comments org :tangle (eval h_func) :noweb yes -qmckl_exit_code qmckl_get_electron_ee_distance_rescaled_deriv_e(qmckl_context context, double* const distance_rescaled_deriv_e); - #+end_src - - #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -qmckl_exit_code qmckl_get_electron_ee_distance_rescaled_deriv_e(qmckl_context context, double* const distance_rescaled_deriv_e) -{ - if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { - return QMCKL_NULL_CONTEXT; - } - - qmckl_exit_code rc; - - rc = qmckl_provide_ee_distance_rescaled_deriv_e(context); - if (rc != QMCKL_SUCCESS) return rc; - - qmckl_context_struct* const ctx = (qmckl_context_struct*) context; - assert (ctx != NULL); - - size_t sze = 4 * ctx->electron.num * ctx->electron.num * ctx->electron.walker.num; - memcpy(distance_rescaled_deriv_e, ctx->electron.ee_distance_rescaled_deriv_e, sze * sizeof(double)); - - return QMCKL_SUCCESS; -} - #+end_src - -*** Provide :noexport: - - #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none -qmckl_exit_code qmckl_provide_ee_distance_rescaled_deriv_e(qmckl_context context); - #+end_src - - #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -qmckl_exit_code qmckl_provide_ee_distance_rescaled_deriv_e(qmckl_context context) -{ - - if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { - return QMCKL_NULL_CONTEXT; - } - - qmckl_context_struct* const ctx = (qmckl_context_struct*) context; - assert (ctx != NULL); - - - /* Compute if necessary */ - if (ctx->electron.walker.point.date > ctx->electron.ee_distance_rescaled_deriv_e_date) { - - if (ctx->electron.walker.num > ctx->electron.walker_old.num) { - free(ctx->electron.ee_distance_rescaled_deriv_e); - ctx->electron.ee_distance_rescaled_deriv_e = NULL; - } - - /* Allocate array */ - if (ctx->electron.ee_distance_rescaled_deriv_e == NULL) { - - qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; - mem_info.size = 4 * ctx->electron.num * ctx->electron.num * - ctx->electron.walker.num * sizeof(double); - double* ee_distance_rescaled_deriv_e = (double*) qmckl_malloc(context, mem_info); - - if (ee_distance_rescaled_deriv_e == NULL) { - return qmckl_failwith( context, - QMCKL_ALLOCATION_FAILED, - "qmckl_ee_distance_rescaled_deriv_e", - NULL); - } - ctx->electron.ee_distance_rescaled_deriv_e = ee_distance_rescaled_deriv_e; - } - - qmckl_exit_code rc = - qmckl_compute_ee_distance_rescaled_deriv_e(context, - ctx->electron.num, - ctx->electron.rescale_factor_kappa_ee, - ctx->electron.walker.num, - ctx->electron.walker.point.coord.data, - ctx->electron.ee_distance_rescaled_deriv_e); - if (rc != QMCKL_SUCCESS) { - return rc; - } - - ctx->electron.ee_distance_rescaled_date = ctx->date; - } - - return QMCKL_SUCCESS; -} - #+end_src - -*** Compute - :PROPERTIES: - :Name: qmckl_compute_ee_distance_rescaled_deriv_e - :CRetType: qmckl_exit_code - :FRetType: qmckl_exit_code - :END: - - #+NAME: qmckl_ee_distance_rescaled_deriv_e_args - | Variable | Type | In/Out | Description | - |---------------------------+-------------------------------------------+--------+-------------------------------------------------| - | ~context~ | ~qmckl_context~ | in | Global state | - | ~elec_num~ | ~int64_t~ | in | Number of electrons | - | ~rescale_factor_kappa_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_deriv_e~ | ~double[walk_num][4][elec_num][elec_num]~ | out | Electron-electron rescaled distance derivatives | - - #+begin_src f90 :comments org :tangle (eval f) :noweb yes -integer function qmckl_compute_ee_distance_rescaled_deriv_e_f(context, elec_num, rescale_factor_kappa_ee, walk_num, & - coord, ee_distance_rescaled_deriv_e) & - result(info) - use qmckl - implicit none - integer(qmckl_context), intent(in) :: context - integer*8 , intent(in) :: elec_num - double precision , intent(in) :: rescale_factor_kappa_ee - integer*8 , intent(in) :: walk_num - double precision , intent(in) :: coord(elec_num,walk_num,3) - double precision , intent(out) :: ee_distance_rescaled_deriv_e(4,elec_num,elec_num,walk_num) - - integer*8 :: k - - info = QMCKL_SUCCESS - - if (context == QMCKL_NULL_CONTEXT) then - info = QMCKL_INVALID_CONTEXT - return - endif - - if (elec_num <= 0) then - info = QMCKL_INVALID_ARG_2 - return - endif - - if (walk_num <= 0) then - info = QMCKL_INVALID_ARG_3 - return - endif - - do k=1,walk_num - info = qmckl_distance_rescaled_deriv_e(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_deriv_e(1,1,1,k), elec_num, rescale_factor_kappa_ee) - if (info /= QMCKL_SUCCESS) then - exit - endif - end do - -end function qmckl_compute_ee_distance_rescaled_deriv_e_f - #+end_src - - #+begin_src c :tangle (eval h_private_func) :comments org :exports none -qmckl_exit_code qmckl_compute_ee_distance_rescaled_deriv_e ( - const qmckl_context context, - const int64_t elec_num, - const double rescale_factor_kappa_ee, - const int64_t walk_num, - const double* coord, - double* const ee_distance_rescaled_deriv_e ); - #+end_src - - #+CALL: generate_c_interface(table=qmckl_ee_distance_rescaled_deriv_e_args,rettyp=get_value("CRetType"),fname=get_value("Name")) - - #+RESULTS: - #+begin_src f90 :tangle (eval f) :comments org :exports none - integer(c_int32_t) function qmckl_compute_ee_distance_rescaled_deriv_e & - (context, elec_num, rescale_factor_kappa_ee, walk_num, coord, ee_distance_rescaled_deriv_e) & - bind(C) result(info) - - use, intrinsic :: iso_c_binding - implicit none - - integer (c_int64_t) , intent(in) , value :: context - integer (c_int64_t) , intent(in) , value :: elec_num - real (c_double ) , intent(in) , value :: rescale_factor_kappa_ee - integer (c_int64_t) , intent(in) , value :: walk_num - real (c_double ) , intent(in) :: coord(elec_num,3,walk_num) - real (c_double ) , intent(out) :: ee_distance_rescaled_deriv_e(4,elec_num,elec_num,walk_num) - - integer(c_int32_t), external :: qmckl_compute_ee_distance_rescaled_deriv_e_f - info = qmckl_compute_ee_distance_rescaled_deriv_e_f & - (context, elec_num, rescale_factor_kappa_ee, walk_num, coord, ee_distance_rescaled_deriv_e) - - end function qmckl_compute_ee_distance_rescaled_deriv_e - #+end_src - -*** Test - - #+begin_src python :results output :exports none -import numpy as np - -# TODO - #+end_src - - #+begin_src c :tangle (eval c_test) -assert(qmckl_electron_provided(context)); - - -double ee_distance_rescaled_deriv_e[4 * walk_num * elec_num * elec_num]; -rc = qmckl_get_electron_ee_distance_rescaled_deriv_e(context, ee_distance_rescaled_deriv_e); - -// 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); - - #+end_src - ** Electron-electron potential ~ee_potential~ is given by @@ -1653,7 +969,6 @@ qmckl_exit_code qmckl_provide_ee_potential(qmckl_context context) if (ctx->electron.walker.num > ctx->electron.walker_old.num) { free(ctx->electron.ee_potential); - ctx->electron.ee_distance_rescaled_deriv_e = NULL; } /* Allocate array */ @@ -1696,13 +1011,13 @@ qmckl_exit_code qmckl_provide_ee_potential(qmckl_context context) :END: #+NAME: qmckl_ee_potential_args - | Variable | Type | In/Out | Description | - |----------------+----------------------------------------+--------+--------------------------------------| - | ~context~ | ~qmckl_context~ | in | Global state | - | ~elec_num~ | ~int64_t~ | in | Number of electrons | - | ~walk_num~ | ~int64_t~ | in | Number of walkers | - | ~ee_distance~ | ~double[walk_num][elec_num][elec_num]~ | in | Electron-electron rescaled distances | - | ~ee_potential~ | ~double[walk_num]~ | out | Electron-electron potential | + | Variable | Type | In/Out | Description | + |----------------+----------------------------------------+--------+-----------------------------| + | ~context~ | ~qmckl_context~ | in | Global state | + | ~elec_num~ | ~int64_t~ | in | Number of electrons | + | ~walk_num~ | ~int64_t~ | in | Number of walkers | + | ~ee_distance~ | ~double[walk_num][elec_num][elec_num]~ | in | Electron-electron distances | + | ~ee_potential~ | ~double[walk_num]~ | out | Electron-electron potential | #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_ee_potential_f(context, elec_num, walk_num, & @@ -1790,7 +1105,7 @@ end function qmckl_compute_ee_potential_f double ee_potential[walk_num]; rc = qmckl_get_electron_ee_potential(context, &(ee_potential[0])); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); #+end_src ** Electron-nucleus distances @@ -2043,20 +1358,20 @@ assert(!qmckl_nucleus_provided(context)); assert(qmckl_electron_provided(context)); rc = qmckl_set_nucleus_num (context, nucl_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); rc = qmckl_set_nucleus_charge (context, charge, nucl_num); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); rc = qmckl_set_nucleus_coord (context, 'T', nucl_coord, 3*nucl_num); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(qmckl_nucleus_provided(context)); double en_distance[walk_num][nucl_num][elec_num]; rc = qmckl_get_electron_en_distance(context, &(en_distance[0][0][0])); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); // (e,n,w) in Fortran notation // (1,1,1) @@ -2079,582 +1394,6 @@ assert(fabs(en_distance[1][0][1] - 3.1804527583077356) < 1.e-12); #+end_src -** Electron-nucleus rescaled distances - - ~en_distance_rescaled~ stores the matrix of the rescaled distances between - electrons and nuclei. - - \[ - C_{ij} = \left( 1 - \exp{-\kappa C_{ij}}\right)/\kappa - \] - - where \(C_{ij}\) is the matrix of electron-nucleus distances. - -*** Get - - #+begin_src c :comments org :tangle (eval h_func) :noweb yes -qmckl_exit_code qmckl_get_electron_en_distance_rescaled(qmckl_context context, double* distance_rescaled); - #+end_src - - - #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -qmckl_exit_code qmckl_get_electron_en_distance_rescaled(qmckl_context context, double* distance_rescaled) -{ - - if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { - return QMCKL_NULL_CONTEXT; - } - - qmckl_exit_code rc; - - rc = qmckl_provide_en_distance_rescaled(context); - if (rc != QMCKL_SUCCESS) return rc; - - qmckl_context_struct* const ctx = (qmckl_context_struct*) context; - assert (ctx != NULL); - - size_t sze = ctx->electron.num * ctx->nucleus.num * ctx->electron.walker.num; - memcpy(distance_rescaled, ctx->electron.en_distance_rescaled, sze * sizeof(double)); - - return QMCKL_SUCCESS; -} - #+end_src - -*** Provide :noexport: - - #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none -qmckl_exit_code qmckl_provide_en_distance_rescaled(qmckl_context context); - #+end_src - - #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -qmckl_exit_code qmckl_provide_en_distance_rescaled(qmckl_context context) -{ - - if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { - return QMCKL_NULL_CONTEXT; - } - - qmckl_context_struct* const ctx = (qmckl_context_struct*) context; - assert (ctx != NULL); - - if (!(ctx->nucleus.provided)) { - return QMCKL_NOT_PROVIDED; - } - - /* Compute if necessary */ - if (ctx->electron.walker.point.date > ctx->electron.en_distance_rescaled_date) { - - if (ctx->electron.walker.num > ctx->electron.walker_old.num) { - free(ctx->electron.en_distance_rescaled); - ctx->electron.en_distance_rescaled = NULL; - } - - /* Allocate array */ - if (ctx->electron.en_distance_rescaled == NULL) { - - qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; - mem_info.size = ctx->electron.num * ctx->nucleus.num * - ctx->electron.walker.num * sizeof(double); - double* en_distance_rescaled = (double*) qmckl_malloc(context, mem_info); - - if (en_distance_rescaled == NULL) { - return qmckl_failwith( context, - QMCKL_ALLOCATION_FAILED, - "qmckl_en_distance_rescaled", - NULL); - } - ctx->electron.en_distance_rescaled = en_distance_rescaled; - } - - qmckl_exit_code rc = - qmckl_compute_en_distance_rescaled(context, - ctx->electron.num, - ctx->nucleus.num, - ctx->electron.rescale_factor_kappa_en, - ctx->electron.walker.num, - ctx->electron.walker.point.coord.data, - ctx->nucleus.coord.data, - ctx->electron.en_distance_rescaled); - if (rc != QMCKL_SUCCESS) { - return rc; - } - - ctx->electron.en_distance_rescaled_date = ctx->date; - } - - return QMCKL_SUCCESS; -} - #+end_src - -*** Compute - :PROPERTIES: - :Name: qmckl_compute_en_distance_rescaled - :CRetType: qmckl_exit_code - :FRetType: qmckl_exit_code - :END: - - #+NAME: qmckl_en_distance_rescaled_args - | 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 | - | ~rescale_factor_kappa_en~ | ~double[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 | - - #+begin_src f90 :comments org :tangle (eval f) :noweb yes -integer function qmckl_compute_en_distance_rescaled_f(context, elec_num, nucl_num, rescale_factor_kappa_en, walk_num, elec_coord, & - nucl_coord, en_distance_rescaled) & - result(info) - use qmckl - implicit none - integer(qmckl_context), intent(in) :: context - integer*8 , intent(in) :: elec_num - integer*8 , intent(in) :: nucl_num - double precision , intent(in) :: rescale_factor_kappa_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(elec_num,nucl_num,walk_num) - - integer*8 :: 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', 'T', elec_num, 1, & - elec_coord(1,k,1), elec_num*walk_num, coord, 1, & - en_distance_rescaled(1,i,k), elec_num, rescale_factor_kappa_en(i)) - if (info /= QMCKL_SUCCESS) then - return - endif - end do - end do - -end function qmckl_compute_en_distance_rescaled_f - #+end_src - - #+begin_src c :tangle (eval h_private_func) :comments org :exports none -qmckl_exit_code qmckl_compute_en_distance_rescaled ( - const qmckl_context context, - const int64_t elec_num, - const int64_t nucl_num, - const double* rescale_factor_kappa_en, - const int64_t walk_num, - const double* elec_coord, - const double* nucl_coord, - double* const en_distance_rescaled ); - #+end_src - - #+CALL: generate_c_interface(table=qmckl_en_distance_rescaled_args,rettyp=get_value("CRetType"),fname=get_value("Name")) - - #+RESULTS: - #+begin_src f90 :tangle (eval f) :comments org :exports none - integer(c_int32_t) function qmckl_compute_en_distance_rescaled & - (context, & - elec_num, & - nucl_num, & - rescale_factor_kappa_en, & - walk_num, & - elec_coord, & - nucl_coord, & - en_distance_rescaled) & - bind(C) result(info) - - use, intrinsic :: iso_c_binding - implicit none - - integer (c_int64_t) , intent(in) , value :: context - integer (c_int64_t) , intent(in) , value :: elec_num - integer (c_int64_t) , intent(in) , value :: nucl_num - real (c_double ) , intent(in) :: rescale_factor_kappa_en(nucl_num) - integer (c_int64_t) , intent(in) , value :: walk_num - real (c_double ) , intent(in) :: elec_coord(elec_num,walk_num,3) - real (c_double ) , intent(in) :: nucl_coord(elec_num,3) - real (c_double ) , intent(out) :: en_distance_rescaled(elec_num,nucl_num,walk_num) - - integer(c_int32_t), external :: qmckl_compute_en_distance_rescaled_f - info = qmckl_compute_en_distance_rescaled_f & - (context, & - elec_num, & - nucl_num, & - rescale_factor_kappa_en, & - walk_num, & - elec_coord, & - nucl_coord, & - en_distance_rescaled) - - end function qmckl_compute_en_distance_rescaled - #+end_src - -*** Test - - #+begin_src python :results output :exports none -import numpy as np - -kappa = 1.0 - -elec_1_w1 = np.array( [ -2.26995253563, -5.15737533569, -2.22940072417 ]) -elec_2_w1 = np.array( [ 3.51983380318, -1.08717381954, -1.19617708027 ]) -elec_1_w2 = np.array( [ -2.34410619736, -3.20016115904, -1.53496759012 ]) -elec_2_w2 = np.array( [ 3.17996025085, -1.40260577202, 1.49473607540 ]) -nucl_1 = np.array( [ 1.096243353458458e+00, 8.907054016973815e-01, 7.777092280258892e-01 ] ) -nucl_2 = np.array( [ 1.168459237342663e+00, 1.125660720053393e+00, 2.833370314829343e+00 ] ) - -print ( "[0][0][0] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_1_w1-nucl_1)) )/kappa ) -print ( "[0][1][0] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_1_w1-nucl_2)) )/kappa ) -print ( "[0][0][1] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_2_w1-nucl_1)) )/kappa ) -print ( "[1][0][0] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_1_w2-nucl_1)) )/kappa ) -print ( "[1][1][0] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_1_w2-nucl_2)) )/kappa ) -print ( "[1][0][1] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_2_w2-nucl_1)) )/kappa ) - - #+end_src - - #+RESULTS: - : [0][0][0] : 0.9994721712909764 - : [0][1][0] : 0.9998448354439821 - : [0][0][1] : 0.9752498074577688 - : [1][0][0] : 0.9970444172399963 - : [1][1][0] : 0.9991586325813303 - : [1][0][1] : 0.9584331688679852 - - #+begin_src c :tangle (eval c_test) - -assert(qmckl_electron_provided(context)); -assert(qmckl_nucleus_provided(context)); - -double en_distance_rescaled[walk_num][nucl_num][elec_num]; - -rc = 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.9994721712909764) < 1.e-12); - -// (1,2,1) -assert(fabs(en_distance_rescaled[0][1][0] - 0.9998448354439821) < 1.e-12); - -// (2,1,1) -assert(fabs(en_distance_rescaled[0][0][1] - 0.9752498074577688) < 1.e-12); - -// (1,1,2) -assert(fabs(en_distance_rescaled[1][0][0] - 0.9970444172399963) < 1.e-12); - -// (1,2,2) -assert(fabs(en_distance_rescaled[1][1][0] - 0.9991586325813303) < 1.e-12); - -// (2,1,2) -assert(fabs(en_distance_rescaled[1][0][1] - 0.9584331688679852) < 1.e-12); - - #+end_src - -** Electron-nucleus rescaled distance gradients and laplacian with respect to electron coords - - The rescaled distances which is given as $R = (1 - \exp{-\kappa r})/\kappa$ - needs to be perturbed with respect to the nuclear coordinates. - This data is stored in the ~en_distance_rescaled_deriv_e~ tensor. The - The first three elements of this three index tensor ~[4][nucl_num][elec_num]~ gives the - derivatives in the x, y, and z directions $dx, dy, dz$ and the last index - gives the Laplacian $\partial x^2 + \partial y^2 + \partial z^2$. - -*** Get - - #+begin_src c :comments org :tangle (eval h_func) :noweb yes -qmckl_exit_code qmckl_get_electron_en_distance_rescaled_deriv_e(qmckl_context context, double* distance_rescaled_deriv_e); - #+end_src - - #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -qmckl_exit_code qmckl_get_electron_en_distance_rescaled_deriv_e(qmckl_context context, double* distance_rescaled_deriv_e) -{ - - if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { - return QMCKL_NULL_CONTEXT; - } - - qmckl_exit_code rc; - - rc = qmckl_provide_en_distance_rescaled_deriv_e(context); - if (rc != QMCKL_SUCCESS) return rc; - - qmckl_context_struct* const ctx = (qmckl_context_struct*) context; - assert (ctx != NULL); - - size_t sze = 4 * ctx->electron.num * ctx->nucleus.num * ctx->electron.walker.num; - memcpy(distance_rescaled_deriv_e, ctx->electron.en_distance_rescaled_deriv_e, sze * sizeof(double)); - - return QMCKL_SUCCESS; -} - #+end_src - -*** Provide :noexport: - - #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none -qmckl_exit_code qmckl_provide_en_distance_rescaled_deriv_e(qmckl_context context); - #+end_src - - #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -qmckl_exit_code qmckl_provide_en_distance_rescaled_deriv_e(qmckl_context context) -{ - - if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { - return QMCKL_NULL_CONTEXT; - } - - qmckl_context_struct* const ctx = (qmckl_context_struct*) context; - assert (ctx != NULL); - - if (!(ctx->nucleus.provided)) { - return QMCKL_NOT_PROVIDED; - } - - /* Compute if necessary */ - if (ctx->electron.walker.point.date > ctx->electron.en_distance_rescaled_deriv_e_date) { - - if (ctx->electron.walker.num > ctx->electron.walker_old.num) { - free(ctx->electron.en_distance_rescaled_deriv_e); - ctx->electron.en_distance_rescaled_deriv_e = NULL; - } - - /* Allocate array */ - if (ctx->electron.en_distance_rescaled_deriv_e == NULL) { - - qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; - mem_info.size = 4 * ctx->electron.num * ctx->nucleus.num * - ctx->electron.walker.num * sizeof(double); - double* en_distance_rescaled_deriv_e = (double*) qmckl_malloc(context, mem_info); - - if (en_distance_rescaled_deriv_e == NULL) { - return qmckl_failwith( context, - QMCKL_ALLOCATION_FAILED, - "qmckl_en_distance_rescaled_deriv_e", - NULL); - } - ctx->electron.en_distance_rescaled_deriv_e = en_distance_rescaled_deriv_e; - } - - qmckl_exit_code rc = - qmckl_compute_en_distance_rescaled_deriv_e(context, - ctx->electron.num, - ctx->nucleus.num, - ctx->electron.rescale_factor_kappa_en, - ctx->electron.walker.num, - ctx->electron.walker.point.coord.data, - ctx->nucleus.coord.data, - ctx->electron.en_distance_rescaled_deriv_e); - if (rc != QMCKL_SUCCESS) { - return rc; - } - - ctx->electron.en_distance_rescaled_deriv_e_date = ctx->date; - } - - return QMCKL_SUCCESS; -} - #+end_src - -*** Compute - :PROPERTIES: - :Name: qmckl_compute_en_distance_rescaled_deriv_e - :CRetType: qmckl_exit_code - :FRetType: qmckl_exit_code - :END: - - #+NAME: qmckl_en_distance_rescaled_deriv_e_args - | 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 | - | ~rescale_factor_kappa_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_deriv_e~ | ~double[walk_num][nucl_num][elec_num][4]~ | out | Electron-nucleus distance derivatives | - - #+begin_src f90 :comments org :tangle (eval f) :noweb yes -integer function qmckl_compute_en_distance_rescaled_deriv_e_f(context, elec_num, nucl_num, & - rescale_factor_kappa_en, walk_num, elec_coord, & - nucl_coord, en_distance_rescaled_deriv_e) & - result(info) - use qmckl - implicit none - integer(qmckl_context), intent(in) :: context - integer*8 , intent(in) :: elec_num - integer*8 , intent(in) :: nucl_num - double precision , intent(in) :: rescale_factor_kappa_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_deriv_e(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_deriv_e(context, 'T', 'T', elec_num, 1, & - elec_coord(1,k,1), elec_num*walk_num, & - coord, 1, & - en_distance_rescaled_deriv_e(1,1,i,k), elec_num, rescale_factor_kappa_en(i)) - if (info /= QMCKL_SUCCESS) then - return - endif - end do - end do - -end function qmckl_compute_en_distance_rescaled_deriv_e_f - #+end_src - - #+begin_src c :tangle (eval h_private_func) :comments org :exports none -qmckl_exit_code qmckl_compute_en_distance_rescaled_deriv_e ( - const qmckl_context context, - const int64_t elec_num, - const int64_t nucl_num, - const double* rescale_factor_kappa_en, - const int64_t walk_num, - const double* elec_coord, - const double* nucl_coord, - double* const en_distance_rescaled_deriv_e ); - #+end_src - - #+CALL: generate_c_interface(table=qmckl_en_distance_rescaled_deriv_e_args,rettyp=get_value("CRetType"),fname=get_value("Name")) - - #+RESULTS: - #+begin_src f90 :tangle (eval f) :comments org :exports none - integer(c_int32_t) function qmckl_compute_en_distance_rescaled_deriv_e & - (context, & - elec_num, & - nucl_num, & - rescale_factor_kappa_en, & - walk_num, & - elec_coord, & - nucl_coord, & - en_distance_rescaled_deriv_e) & - bind(C) result(info) - - use, intrinsic :: iso_c_binding - implicit none - - integer (c_int64_t) , intent(in) , value :: context - integer (c_int64_t) , intent(in) , value :: elec_num - integer (c_int64_t) , intent(in) , value :: nucl_num - real (c_double ) , intent(in) :: rescale_factor_kappa_en(nucl_num) - integer (c_int64_t) , intent(in) , value :: walk_num - real (c_double ) , intent(in) :: elec_coord(elec_num,walk_num,3) - real (c_double ) , intent(in) :: nucl_coord(elec_num,3) - real (c_double ) , intent(out) :: en_distance_rescaled_deriv_e(4,elec_num,nucl_num,walk_num) - - integer(c_int32_t), external :: qmckl_compute_en_distance_rescaled_deriv_e_f - info = qmckl_compute_en_distance_rescaled_deriv_e_f & - (context, & - elec_num, & - nucl_num, & - rescale_factor_kappa_en, & - walk_num, & - elec_coord, & - nucl_coord, & - en_distance_rescaled_deriv_e) - - end function qmckl_compute_en_distance_rescaled_deriv_e - #+end_src - -*** Test - - #+begin_src python :results output :exports none -import numpy as np - -# TODO - #+end_src - - - #+begin_src c :tangle (eval c_test) - -assert(qmckl_electron_provided(context)); - -rc = qmckl_set_nucleus_rescale_factor (context, nucl_rescale_factor_kappa); -assert(rc == QMCKL_SUCCESS); - -assert(qmckl_nucleus_provided(context)); - -double en_distance_rescaled_deriv_e[walk_num][4][nucl_num][elec_num]; - -rc = qmckl_get_electron_en_distance_rescaled_deriv_e(context, &(en_distance_rescaled_deriv_e[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); - - #+end_src - ** Electron-nucleus potential ~en_potential~ stores the ~en~ potential energy @@ -2868,7 +1607,7 @@ end function qmckl_compute_en_potential_f double en_potential[walk_num]; rc = qmckl_get_electron_en_potential(context, &(en_potential[0])); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); #+end_src ** Generate initial coordinates diff --git a/org/qmckl_error.org b/org/qmckl_error.org index 33829e4..a80db17 100644 --- a/org/qmckl_error.org +++ b/org/qmckl_error.org @@ -608,6 +608,62 @@ qmckl_last_error(qmckl_context context, char* buffer) { end interface #+end_src +* Helper functions for debugging + + The following function prints to ~stderr~ an error message is the return code is + not ~QMCKL_SUCCESS~. + + # Header + #+begin_src c :comments org :tangle (eval h_func) +qmckl_exit_code +qmckl_check(qmckl_context context, qmckl_exit_code rc); + #+end_src + + # Source + #+begin_src c :tangle (eval c) :exports none +#include + +qmckl_exit_code +qmckl_check(qmckl_context context, qmckl_exit_code rc) +{ + + char fname[QMCKL_MAX_FUN_LEN]; + char message[QMCKL_MAX_MSG_LEN]; + + if (rc != QMCKL_SUCCESS) { + fprintf(stderr, "===========\nQMCKL ERROR\n%s\n", qmckl_string_of_error(rc)); + qmckl_get_error(context, &rc, fname, message); + fprintf(stderr, "Function: %s\nMessage: %s\n===========\n", fname, message); + } + + return rc; +} + #+end_src + + It should be used as: + + #+begin_src c +rc = qmckl_check(context, + qmckl_...(context, ...) + ); +assert (rc == QMCKL_SUCCESS); + #+end_src + +** Fortran inteface + + #+begin_src f90 :tangle (eval fh_func) :exports none :noweb yes + interface + function qmckl_check (context, rc) bind(C, name='qmckl_check') + use, intrinsic :: iso_c_binding + import + implicit none + integer(qmckl_exit_code) :: qmckl_check + integer (c_int64_t) , intent(in), value :: context + integer(qmckl_exit_code), intent(in) :: rc + end function qmckl_check + end interface + #+end_src + * End of files :noexport: #+begin_src c :comments link :tangle (eval h_private_type) @@ -623,12 +679,13 @@ qmckl_last_error(qmckl_context context, char* buffer) { qmckl_exit_code exit_code; exit_code = 1; - assert (qmckl_set_error(context, exit_code, "qmckl_transpose", "Success") == QMCKL_SUCCESS); + assert (qmckl_set_error(context, exit_code, "my_function", "Message") == QMCKL_SUCCESS); assert (qmckl_get_error(context, &exit_code, function_name, message) == QMCKL_SUCCESS); assert (exit_code == 1); - assert (strcmp(function_name,"qmckl_transpose") == 0); - assert (strcmp(message,"Success") == 0); + + assert (strcmp(function_name,"my_function") == 0); + assert (strcmp(message,"Message") == 0); exit_code = qmckl_context_destroy(context); assert(exit_code == QMCKL_SUCCESS); diff --git a/org/qmckl_jastrow.org b/org/qmckl_jastrow.org index f06d00b..f078f66 100644 --- a/org/qmckl_jastrow.org +++ b/org/qmckl_jastrow.org @@ -137,6 +137,8 @@ int main() { | Variable | Type | In/Out | Description | |---------------------------+---------------------------------------+--------+-------------------------------------------------------------------| | ~uninitialized~ | ~int32_t~ | in | Keeps bits set for uninitialized data | + | ~rescale_factor_ee~ | ~double~ | in | The distance scaling factor | + | ~rescale_factor_en~ | ~double[nucl_num]~ | in | The distance scaling factor | | ~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 | @@ -160,24 +162,32 @@ int main() { computed data: - | Variable | Type | In/Out | - |-------------------------------+-------------------------------------------------------------------+---------------------------------------------------------------------------------------------------------| - | ~dim_cord_vect~ | ~int64_t~ | Number of unique C coefficients | - | ~dim_cord_vect_date~ | ~uint64_t~ | Number of unique C coefficients | - | ~asymp_jasb~ | ~double[2]~ | Asymptotic component | - | ~asymp_jasb_date~ | ~uint64_t~ | Asymptotic component | - | ~cord_vect_full~ | ~double[dim_cord_vect][nucl_num]~ | vector of non-zero coefficients | - | ~cord_vect_full_date~ | ~uint64_t~ | Keep track of changes here | - | ~lkpm_combined_index~ | ~int64_t[4][dim_cord_vect]~ | Transform l,k,p, and m into consecutive indices | - | ~lkpm_combined_index_date~ | ~uint64_t~ | Transform l,k,p, and m into consecutive indices | - | ~tmp_c~ | ~double[walker.num][cord_num][cord_num+1][nucl_num][elec_num]~ | vector of non-zero coefficients | - | ~dtmp_c~ | ~double[walker.num][elec_num][4][nucl_num][cord_num+1][cord_num]~ | vector of non-zero coefficients | - | ~een_rescaled_n~ | ~double[walker.num][cord_num+1][nucl_num][elec_num]~ | The electron-electron rescaled distances raised to the powers defined by cord | - | ~een_rescaled_n_date~ | ~uint64_t~ | Keep track of the date of creation | - | ~een_rescaled_e_deriv_e~ | ~double[walker.num][cord_num+1][elec_num][4][elec_num]~ | The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons | - | ~een_rescaled_e_deriv_e_date~ | ~uint64_t~ | Keep track of the date of creation | - | ~een_rescaled_n_deriv_e~ | ~double[walker.num][cord_num+1][nucl_num][4][elec_num]~ | The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons | - | ~een_rescaled_n_deriv_e_date~ | ~uint64_t~ | Keep track of the date of creation | + | Variable | Type | In/Out | + |-------------------------------------+-------------------------------------------------------------------+---------------------------------------------------------------------------------------------------------| + | ~dim_cord_vect~ | ~int64_t~ | Number of unique C coefficients | + | ~dim_cord_vect_date~ | ~uint64_t~ | Number of unique C coefficients | + | ~asymp_jasb~ | ~double[2]~ | Asymptotic component | + | ~asymp_jasb_date~ | ~uint64_t~ | Asymptotic component | + | ~cord_vect_full~ | ~double[dim_cord_vect][nucl_num]~ | vector of non-zero coefficients | + | ~cord_vect_full_date~ | ~uint64_t~ | Keep track of changes here | + | ~lkpm_combined_index~ | ~int64_t[4][dim_cord_vect]~ | Transform l,k,p, and m into consecutive indices | + | ~lkpm_combined_index_date~ | ~uint64_t~ | Transform l,k,p, and m into consecutive indices | + | ~tmp_c~ | ~double[walker.num][cord_num][cord_num+1][nucl_num][elec_num]~ | vector of non-zero coefficients | + | ~dtmp_c~ | ~double[walker.num][elec_num][4][nucl_num][cord_num+1][cord_num]~ | vector of non-zero coefficients | + | ~ee_distance_rescaled~ | ~double[walker.num][num][num]~ | Electron-electron rescaled distances | + | ~ee_distance_rescaled_date~ | ~uint64_t~ | Last modification date of the electron-electron distances | + | ~ee_distance_rescaled_deriv_e~ | ~double[walker.num][4][num][num]~ | Electron-electron rescaled distances derivatives | + | ~ee_distance_rescaled_deriv_e_date~ | ~uint64_t~ | Last modification date of the electron-electron distance derivatives | + | ~en_distance_rescaled~ | ~double[walker.num][nucl_num][num]~ | Electron-nucleus distances | + | ~en_distance_rescaled_date~ | ~uint64_t~ | Last modification date of the electron-electron distances | + | ~en_distance_rescaled_deriv_e~ | ~double[walker.num][4][nucl_num][num]~ | Electron-electron rescaled distances derivatives | + | ~en_distance_rescaled_deriv_e_date~ | ~uint64_t~ | Last modification date of the electron-electron distance derivatives | + | ~een_rescaled_n~ | ~double[walker.num][cord_num+1][nucl_num][elec_num]~ | The electron-electron rescaled distances raised to the powers defined by cord | + | ~een_rescaled_n_date~ | ~uint64_t~ | Keep track of the date of creation | + | ~een_rescaled_e_deriv_e~ | ~double[walker.num][cord_num+1][elec_num][4][elec_num]~ | The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons | + | ~een_rescaled_e_deriv_e_date~ | ~uint64_t~ | Keep track of the date of creation | + | ~een_rescaled_n_deriv_e~ | ~double[walker.num][cord_num+1][nucl_num][4][elec_num]~ | The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons | + | ~een_rescaled_n_deriv_e_date~ | ~uint64_t~ | Keep track of the date of creation | #+NAME: jastrow_data #+BEGIN_SRC python :results none :exports none @@ -348,6 +358,8 @@ typedef struct qmckl_jastrow_struct{ uint64_t factor_ee_deriv_e_date; uint64_t factor_en_deriv_e_date; uint64_t factor_een_deriv_e_date; + double rescale_factor_ee; + double* rescale_factor_en; int64_t* type_nucl_vector; double * aord_vector; double * bord_vector; @@ -367,6 +379,14 @@ typedef struct qmckl_jastrow_struct{ uint64_t lkpm_combined_index_date; double * tmp_c; double * dtmp_c; + uint64_t ee_distance_rescaled_date; + uint64_t ee_distance_rescaled_deriv_e_date; + uint64_t en_distance_rescaled_date; + uint64_t en_distance_rescaled_deriv_e_date; + double* ee_distance_rescaled; + double* ee_distance_rescaled_deriv_e; + double* en_distance_rescaled; + double* en_distance_rescaled_deriv_e; double * een_rescaled_e; double * een_rescaled_n; uint64_t een_rescaled_e_date; @@ -406,7 +426,7 @@ qmckl_exit_code qmckl_init_jastrow(qmckl_context context) { qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); - ctx->jastrow.uninitialized = (1 << 6) - 1; + ctx->jastrow.uninitialized = (1 << 8) - 1; /* Default values */ @@ -414,6 +434,494 @@ qmckl_exit_code qmckl_init_jastrow(qmckl_context context) { } #+end_src +** Initialization functions + + To prepare for the Jastrow and its derivative, all the following functions need to be + called. + + #+begin_src c :comments org :tangle (eval h_func) +qmckl_exit_code qmckl_set_jastrow_rescale_factor_ee (qmckl_context context, const double kappa_ee); +qmckl_exit_code qmckl_set_jastrow_rescale_factor_en (qmckl_context context, const double* kappa_en, const int64_t size_max); +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); + #+end_src + + #+NAME:pre2 + #+begin_src c :exports none +if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + +qmckl_context_struct* const ctx = (qmckl_context_struct*) context; + +if (mask != 0 && !(ctx->jastrow.uninitialized & mask)) { + printf("%d %d\n", mask, ctx->jastrow.uninitialized ); + return qmckl_failwith( context, + QMCKL_ALREADY_SET, + "qmckl_set_jastrow_*", + NULL); + } + #+end_src + + #+NAME:post2 + #+begin_src c :exports none +ctx->jastrow.uninitialized &= ~mask; +ctx->jastrow.provided = (ctx->jastrow.uninitialized == 0); +if (ctx->jastrow.provided) { + qmckl_exit_code rc_ = qmckl_finalize_jastrow(context); + if (rc_ != QMCKL_SUCCESS) return rc_; + } + +return QMCKL_SUCCESS; + #+end_src + + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +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) +{ + + int32_t mask = 1 << 0; + +<> + + if (aord_num <= 0) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_set_jastrow_ord_num", + "aord_num <= 0"); + } + + if (bord_num <= 0) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_set_jastrow_ord_num", + "bord_num <= 0"); + } + + if (cord_num <= 0) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_set_jastrow_ord_num", + "cord_num <= 0"); + } + + ctx->jastrow.aord_num = aord_num; + ctx->jastrow.bord_num = bord_num; + ctx->jastrow.cord_num = cord_num; + + <> +} + + +qmckl_exit_code +qmckl_set_jastrow_type_nucl_num(qmckl_context context, const int64_t type_nucl_num) +{ + int32_t mask = 1 << 1; + +<> + + if (type_nucl_num <= 0) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_set_jastrow_type_nucl_num", + "type_nucl_num < 0"); + } + + ctx->jastrow.type_nucl_num = type_nucl_num; + + <> +} + + +qmckl_exit_code +qmckl_set_jastrow_type_nucl_vector(qmckl_context context, + int64_t const * type_nucl_vector, + const int64_t nucl_num) +{ + + int32_t mask = 1 << 2; + +<> + + int64_t type_nucl_num; + qmckl_exit_code rc = qmckl_get_jastrow_type_nucl_num(context, &type_nucl_num); + if (rc != QMCKL_SUCCESS) return rc; + + if (type_nucl_num == 0) { + return qmckl_failwith( context, + QMCKL_FAILURE, + "qmckl_set_jastrow_type_nucl_vector", + "type_nucl_num is not set"); + } + + if (type_nucl_vector == NULL) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_set_jastrow_type_nucl_vector", + "type_nucl_vector = NULL"); + } + + if (ctx->jastrow.type_nucl_vector != NULL) { + rc = qmckl_free(context, ctx->jastrow.type_nucl_vector); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_set_type_nucl_vector", + "Unable to free ctx->jastrow.type_nucl_vector"); + } + } + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = nucl_num * sizeof(int64_t); + int64_t* new_array = (int64_t*) qmckl_malloc(context, mem_info); + + if(new_array == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_set_jastrow_type_nucl_vector", + NULL); + } + + memcpy(new_array, type_nucl_vector, mem_info.size); + + ctx->jastrow.type_nucl_vector = new_array; + + <> +} + + +qmckl_exit_code +qmckl_set_jastrow_aord_vector(qmckl_context context, + double const * aord_vector, + const int64_t size_max) +{ + int32_t mask = 1 << 3; + +<> + + int64_t aord_num; + qmckl_exit_code rc = qmckl_get_jastrow_aord_num(context, &aord_num); + if (rc != QMCKL_SUCCESS) return rc; + + int64_t type_nucl_num; + rc = qmckl_get_jastrow_type_nucl_num(context, &type_nucl_num); + if (rc != QMCKL_SUCCESS) return rc; + + if (aord_num == 0) { + return qmckl_failwith( context, + QMCKL_FAILURE, + "qmckl_set_jastrow_coefficient", + "aord_num is not set"); + } + + if (aord_vector == NULL) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_set_jastrow_aord_vector", + "aord_vector = NULL"); + } + + if (ctx->jastrow.aord_vector != NULL) { + rc = qmckl_free(context, ctx->jastrow.aord_vector); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_set_jastrow_aord_vector", + "Unable to free ctx->jastrow.aord_vector"); + } + } + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = (aord_num + 1) * type_nucl_num * sizeof(double); + + if ((size_t) size_max < mem_info.size/sizeof(double)) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_set_jastrow_aord_vector", + "Array too small. Expected (aord_num+1)*type_nucl_num"); + } + + double* new_array = (double*) qmckl_malloc(context, mem_info); + + if(new_array == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_set_jastrow_coefficient", + NULL); + } + + memcpy(new_array, aord_vector, mem_info.size); + + ctx->jastrow.aord_vector = new_array; + + <> +} + + +qmckl_exit_code +qmckl_set_jastrow_bord_vector(qmckl_context context, + double const * bord_vector, + const int64_t size_max) +{ + int32_t mask = 1 << 4; + +<> + + int64_t bord_num; + qmckl_exit_code rc = qmckl_get_jastrow_bord_num(context, &bord_num); + if (rc != QMCKL_SUCCESS) return rc; + + if (bord_num == 0) { + return qmckl_failwith( context, + QMCKL_FAILURE, + "qmckl_set_jastrow_coefficient", + "bord_num is not set"); + } + + if (bord_vector == NULL) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_set_jastrow_bord_vector", + "bord_vector = NULL"); + } + + if (ctx->jastrow.bord_vector != NULL) { + rc = qmckl_free(context, ctx->jastrow.bord_vector); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_set_jastrow_bord_vector", + "Unable to free ctx->jastrow.bord_vector"); + } + } + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = (bord_num + 1) * sizeof(double); + + if ((size_t) size_max < mem_info.size/sizeof(double)) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_set_jastrow_bord_vector", + "Array too small. Expected (bord_num+1)"); + } + + double* new_array = (double*) qmckl_malloc(context, mem_info); + + if(new_array == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_set_jastrow_coefficient", + NULL); + } + + memcpy(new_array, bord_vector, mem_info.size); + + ctx->jastrow.bord_vector = new_array; + + <> +} + + +qmckl_exit_code +qmckl_set_jastrow_cord_vector(qmckl_context context, + double const * cord_vector, + const int64_t size_max) +{ + int32_t mask = 1 << 5; + +<> + + qmckl_exit_code rc = qmckl_provide_dim_cord_vect(context); + if (rc != QMCKL_SUCCESS) return rc; + + int64_t dim_cord_vect; + rc = qmckl_get_jastrow_dim_cord_vect(context, &dim_cord_vect); + if (rc != QMCKL_SUCCESS) return rc; + + int64_t type_nucl_num; + rc = qmckl_get_jastrow_type_nucl_num(context, &type_nucl_num); + if (rc != QMCKL_SUCCESS) return rc; + + if (dim_cord_vect == 0) { + return qmckl_failwith( context, + QMCKL_FAILURE, + "qmckl_set_jastrow_coefficient", + "dim_cord_vect is not set"); + } + + if (cord_vector == NULL) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_set_jastrow_cord_vector", + "cord_vector = NULL"); + } + + if (ctx->jastrow.cord_vector != NULL) { + rc = qmckl_free(context, ctx->jastrow.cord_vector); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_set_jastrow_cord_vector", + "Unable to free ctx->jastrow.cord_vector"); + } + } + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = dim_cord_vect * type_nucl_num * sizeof(double); + + if ((size_t) size_max < mem_info.size/sizeof(double)) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_set_jastrow_cord_vector", + "Array too small. Expected dim_cord_vect * type_nucl_num"); + } + + double* new_array = (double*) qmckl_malloc(context, mem_info); + + if(new_array == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_set_jastrow_coefficient", + NULL); + } + + memcpy(new_array, cord_vector, mem_info.size); + + ctx->jastrow.cord_vector = new_array; + + <> +} + +qmckl_exit_code +qmckl_set_jastrow_rescale_factor_ee(qmckl_context context, + const double rescale_factor_ee) { + + int32_t mask = 1 << 6; + + <> + + if (rescale_factor_ee <= 0.0) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_set_jastrow_rescale_factor_ee", + "rescale_factor_ee <= 0.0"); + } + + ctx->jastrow.rescale_factor_ee = rescale_factor_ee; + + <> +} + + +qmckl_exit_code +qmckl_set_jastrow_rescale_factor_en(qmckl_context context, + const double* rescale_factor_en, + const int64_t size_max) { + + int32_t mask = 1 << 7; + + <> + + if (rescale_factor_en == NULL) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_set_jastrow_rescale_factor_en", + "Null pointer"); + } + + if (size_max < ctx->nucleus.num) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_set_jastrow_rescale_factor_en", + "Array too small"); + } + + + if (ctx->jastrow.rescale_factor_en != NULL) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_set_jastrow_rescale_factor_en", + "Already set"); + } + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->nucleus.num * sizeof(double); + ctx->jastrow.rescale_factor_en = (double*) qmckl_malloc(context, mem_info); + + for (int64_t i=0 ; inucleus.num ; ++i) { + if (rescale_factor_en[i] <= 0.0) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_set_jastrow_rescale_factor_en", + "rescale_factor_en <= 0.0"); + } + ctx->jastrow.rescale_factor_en[i] = rescale_factor_en[i]; + } + + <> +} + #+end_src + + 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. + + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none +qmckl_exit_code qmckl_finalize_jastrow(qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_finalize_jastrow(qmckl_context context) { + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_INVALID_CONTEXT; + } + + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; + assert (ctx != NULL); + + /* ----------------------------------- */ + /* Check for the necessary information */ + /* ----------------------------------- */ + + /* Check for the electron data + 1. elec_num + 2. ee_distances_rescaled + ,*/ + if (!(ctx->electron.provided)) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_electron", + NULL); + } + + /* Check for the nucleus data + 1. nucl_num + 2. en_distances_rescaled + ,*/ + if (!(ctx->nucleus.provided)) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_nucleus", + NULL); + } + + /* Decide if the Jastrow if offloaded on GPU or not */ +#if defined(HAVE_HPC) && (defined(HAVE_CUBLAS_OFFLOAD) || defined(HAVE_OPENACC_OFFLOAD) || defined(HAVE_OPENMP_OFFLOAD)) + ctx->jastrow.gpu_offload = true; // ctx->electron.num > 100; +#endif + + qmckl_exit_code rc = QMCKL_SUCCESS; + return rc; + + +} + #+end_src + ** Access functions #+begin_src c :comments org :tangle (eval h_func) :exports none @@ -425,7 +933,9 @@ qmckl_exit_code qmckl_get_jastrow_type_nucl_vector (qmckl_context context, int qmckl_exit_code qmckl_get_jastrow_aord_vector (qmckl_context context, double * const aord_vector, const int64_t size_max); qmckl_exit_code qmckl_get_jastrow_bord_vector (qmckl_context context, double * const bord_vector, const int64_t size_max); qmckl_exit_code qmckl_get_jastrow_cord_vector (qmckl_context context, double * const cord_vector, const int64_t size_max); - #+end_src +qmckl_exit_code qmckl_get_jastrow_rescale_factor_ee (const qmckl_context context, double* const rescale_factor_ee); +qmckl_exit_code qmckl_get_jastrow_rescale_factor_en (const qmckl_context context, double* const rescale_factor_en, const int64_t size_max); + #+end_src Along with these core functions, calculation of the jastrow factor requires the following additional information to be set: @@ -452,13 +962,6 @@ bool qmckl_jastrow_provided(const qmckl_context context) { } #+end_src - #+NAME:post - #+begin_src c :exports none -if ( (ctx->jastrow.uninitialized & mask) != 0) { - return NULL; -} - #+end_src - #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_jastrow_aord_num (const qmckl_context context, int64_t* const aord_num) { @@ -722,425 +1225,75 @@ qmckl_get_jastrow_cord_vector (const qmckl_context context, return QMCKL_SUCCESS; } - #+end_src - -** Initialization functions - - To prepare for the Jastrow and its derivative, all the following functions need to be - called. - - #+begin_src c :comments org :tangle (eval h_func) -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); - #+end_src - - #+NAME:pre2 - #+begin_src c :exports none -if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { - return QMCKL_NULL_CONTEXT; - } - -qmckl_context_struct* const ctx = (qmckl_context_struct*) context; - -if (mask != 0 && !(ctx->jastrow.uninitialized & mask)) { - printf("%d %d\n", mask, ctx->jastrow.uninitialized ); - return qmckl_failwith( context, - QMCKL_ALREADY_SET, - "qmckl_set_jastrow_*", - NULL); - } - #+end_src - - #+NAME:post2 - #+begin_src c :exports none -ctx->jastrow.uninitialized &= ~mask; -ctx->jastrow.provided = (ctx->jastrow.uninitialized == 0); -if (ctx->jastrow.provided) { - qmckl_exit_code rc_ = qmckl_finalize_jastrow(context); - if (rc_ != QMCKL_SUCCESS) return rc_; - } - -return QMCKL_SUCCESS; - #+end_src - - - #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none 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) -{ - - int32_t mask = 1 << 0; - -<> - - if (aord_num <= 0) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_2, - "qmckl_set_jastrow_ord_num", - "aord_num <= 0"); - } - - if (bord_num <= 0) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_2, - "qmckl_set_jastrow_ord_num", - "bord_num <= 0"); - } - - if (cord_num <= 0) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_2, - "qmckl_set_jastrow_ord_num", - "cord_num <= 0"); - } - - ctx->jastrow.aord_num = aord_num; - ctx->jastrow.bord_num = bord_num; - ctx->jastrow.cord_num = cord_num; - - <> -} - - -qmckl_exit_code -qmckl_set_jastrow_type_nucl_num(qmckl_context context, const int64_t type_nucl_num) -{ - int32_t mask = 1 << 1; - -<> - - if (type_nucl_num <= 0) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_2, - "qmckl_set_jastrow_type_nucl_num", - "type_nucl_num < 0"); - } - - ctx->jastrow.type_nucl_num = type_nucl_num; - - <> -} - - -qmckl_exit_code -qmckl_set_jastrow_type_nucl_vector(qmckl_context context, - int64_t const * type_nucl_vector, - const int64_t nucl_num) -{ - - int32_t mask = 1 << 2; - -<> - - int64_t type_nucl_num; - qmckl_exit_code rc = qmckl_get_jastrow_type_nucl_num(context, &type_nucl_num); - if (rc != QMCKL_SUCCESS) return rc; - - if (type_nucl_num == 0) { - return qmckl_failwith( context, - QMCKL_FAILURE, - "qmckl_set_jastrow_type_nucl_vector", - "type_nucl_num is not set"); - } - - if (type_nucl_vector == NULL) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_2, - "qmckl_set_jastrow_type_nucl_vector", - "type_nucl_vector = NULL"); - } - - if (ctx->jastrow.type_nucl_vector != NULL) { - rc = qmckl_free(context, ctx->jastrow.type_nucl_vector); - if (rc != QMCKL_SUCCESS) { - return qmckl_failwith( context, rc, - "qmckl_set_type_nucl_vector", - NULL); - } - } - - qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; - mem_info.size = nucl_num * sizeof(int64_t); - int64_t* new_array = (int64_t*) qmckl_malloc(context, mem_info); - - if(new_array == NULL) { - return qmckl_failwith( context, - QMCKL_ALLOCATION_FAILED, - "qmckl_set_jastrow_type_nucl_vector", - NULL); - } - - memcpy(new_array, type_nucl_vector, mem_info.size); - - ctx->jastrow.type_nucl_vector = new_array; - - <> -} - - -qmckl_exit_code -qmckl_set_jastrow_aord_vector(qmckl_context context, - double const * aord_vector, - const int64_t size_max) -{ - int32_t mask = 1 << 3; - -<> - - int64_t aord_num; - qmckl_exit_code rc = qmckl_get_jastrow_aord_num(context, &aord_num); - if (rc != QMCKL_SUCCESS) return rc; - - int64_t type_nucl_num; - rc = qmckl_get_jastrow_type_nucl_num(context, &type_nucl_num); - if (rc != QMCKL_SUCCESS) return rc; - - if (aord_num == 0) { - return qmckl_failwith( context, - QMCKL_FAILURE, - "qmckl_set_jastrow_coefficient", - "aord_num is not set"); - } - - if (aord_vector == NULL) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_2, - "qmckl_set_jastrow_aord_vector", - "aord_vector = NULL"); - } - - if (ctx->jastrow.aord_vector != NULL) { - rc = qmckl_free(context, ctx->jastrow.aord_vector); - if (rc != QMCKL_SUCCESS) { - return qmckl_failwith( context, rc, - "qmckl_set_ord_vector", - NULL); - } - } - - qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; - mem_info.size = (aord_num + 1) * type_nucl_num * sizeof(double); - - if ((size_t) size_max < mem_info.size/sizeof(double)) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_3, - "qmckl_set_jastrow_aord_vector", - "Array too small. Expected (aord_num+1)*type_nucl_num"); - } - - double* new_array = (double*) qmckl_malloc(context, mem_info); - - if(new_array == NULL) { - return qmckl_failwith( context, - QMCKL_ALLOCATION_FAILED, - "qmckl_set_jastrow_coefficient", - NULL); - } - - memcpy(new_array, aord_vector, mem_info.size); - - ctx->jastrow.aord_vector = new_array; - - <> -} - - -qmckl_exit_code -qmckl_set_jastrow_bord_vector(qmckl_context context, - double const * bord_vector, - const int64_t size_max) -{ - int32_t mask = 1 << 4; - -<> - - int64_t bord_num; - qmckl_exit_code rc = qmckl_get_jastrow_bord_num(context, &bord_num); - if (rc != QMCKL_SUCCESS) return rc; - - if (bord_num == 0) { - return qmckl_failwith( context, - QMCKL_FAILURE, - "qmckl_set_jastrow_coefficient", - "bord_num is not set"); - } - - if (bord_vector == NULL) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_2, - "qmckl_set_jastrow_bord_vector", - "bord_vector = NULL"); - } - - if (ctx->jastrow.bord_vector != NULL) { - rc = qmckl_free(context, ctx->jastrow.bord_vector); - if (rc != QMCKL_SUCCESS) { - return qmckl_failwith( context, rc, - "qmckl_set_ord_vector", - NULL); - } - } - - qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; - mem_info.size = (bord_num + 1) * sizeof(double); - - if ((size_t) size_max < mem_info.size/sizeof(double)) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_3, - "qmckl_set_jastrow_bord_vector", - "Array too small. Expected (bord_num+1)"); - } - - double* new_array = (double*) qmckl_malloc(context, mem_info); - - if(new_array == NULL) { - return qmckl_failwith( context, - QMCKL_ALLOCATION_FAILED, - "qmckl_set_jastrow_coefficient", - NULL); - } - - memcpy(new_array, bord_vector, mem_info.size); - - ctx->jastrow.bord_vector = new_array; - - <> -} - - -qmckl_exit_code -qmckl_set_jastrow_cord_vector(qmckl_context context, - double const * cord_vector, - const int64_t size_max) -{ - int32_t mask = 1 << 5; - -<> - - qmckl_exit_code rc = qmckl_provide_dim_cord_vect(context); - if (rc != QMCKL_SUCCESS) return rc; - - int64_t dim_cord_vect; - rc = qmckl_get_jastrow_dim_cord_vect(context, &dim_cord_vect); - if (rc != QMCKL_SUCCESS) return rc; - - int64_t type_nucl_num; - rc = qmckl_get_jastrow_type_nucl_num(context, &type_nucl_num); - if (rc != QMCKL_SUCCESS) return rc; - - if (dim_cord_vect == 0) { - return qmckl_failwith( context, - QMCKL_FAILURE, - "qmckl_set_jastrow_coefficient", - "dim_cord_vect is not set"); - } - - if (cord_vector == NULL) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_2, - "qmckl_set_jastrow_cord_vector", - "cord_vector = NULL"); - } - - if (ctx->jastrow.cord_vector != NULL) { - rc = qmckl_free(context, ctx->jastrow.cord_vector); - if (rc != QMCKL_SUCCESS) { - return qmckl_failwith( context, rc, - "qmckl_set_cord_vector", - NULL); - } - } - - qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; - mem_info.size = dim_cord_vect * type_nucl_num * sizeof(double); - - if ((size_t) size_max < mem_info.size/sizeof(double)) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_3, - "qmckl_set_jastrow_cord_vector", - "Array too small. Expected dim_cord_vect * type_nucl_num"); - } - - double* new_array = (double*) qmckl_malloc(context, mem_info); - - if(new_array == NULL) { - return qmckl_failwith( context, - QMCKL_ALLOCATION_FAILED, - "qmckl_set_jastrow_coefficient", - NULL); - } - - memcpy(new_array, cord_vector, mem_info.size); - - ctx->jastrow.cord_vector = new_array; - - <> -} - - #+end_src - - 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. - - #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none -qmckl_exit_code qmckl_finalize_jastrow(qmckl_context context); - #+end_src - - #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -qmckl_exit_code qmckl_finalize_jastrow(qmckl_context context) { +qmckl_get_jastrow_rescale_factor_ee (const qmckl_context context, + double* const rescale_factor_ee) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_INVALID_CONTEXT; } + if (rescale_factor_ee == NULL) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_get_jastrow_rescale_factor_ee", + "rescale_factor_ee is a null pointer"); + } + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); - /* ----------------------------------- */ - /* Check for the necessary information */ - /* ----------------------------------- */ + int32_t mask = 1 << 6; - /* Check for the electron data - 1. elec_num - 2. ee_distances_rescaled - ,*/ - if (!(ctx->electron.provided)) { - return qmckl_failwith( context, - QMCKL_NOT_PROVIDED, - "qmckl_electron", - NULL); + if ( (ctx->jastrow.uninitialized & mask) != 0) { + return QMCKL_NOT_PROVIDED; } + assert (ctx->jastrow.rescale_factor_ee > 0.0); - /* Check for the nucleus data - 1. nucl_num - 2. en_distances_rescaled - ,*/ - if (!(ctx->nucleus.provided)) { - return qmckl_failwith( context, - QMCKL_NOT_PROVIDED, - "qmckl_nucleus", - NULL); - } - - /* Decide if the Jastrow if offloaded on GPU or not */ -#if defined(HAVE_HPC) && (defined(HAVE_CUBLAS_OFFLOAD) || defined(HAVE_OPENACC_OFFLOAD) || defined(HAVE_OPENMP_OFFLOAD)) - ctx->jastrow.gpu_offload = true; // ctx->electron.num > 100; -#endif - - qmckl_exit_code rc = QMCKL_SUCCESS; - return rc; - - + ,*rescale_factor_ee = ctx->jastrow.rescale_factor_ee; + return QMCKL_SUCCESS; } - #+end_src + + +qmckl_exit_code +qmckl_get_jastrow_rescale_factor_en (const qmckl_context context, + double* const rescale_factor_en, + const int64_t size_max) { + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_INVALID_CONTEXT; + } + + if (rescale_factor_en == NULL) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_get_jastrow_rescale_factor_en", + "rescale_factor_en is a null pointer"); + } + + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; + assert (ctx != NULL); + + int32_t mask = 1 << 7; + + if ( (ctx->jastrow.uninitialized & mask) != 0) { + return QMCKL_NOT_PROVIDED; + } + + if (size_max < ctx->nucleus.num) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_get_jastrow_rescale_factor_en", + "Array to small"); + } + + assert(ctx->jastrow.rescale_factor_en != NULL); + for (int64_t i=0 ; inucleus.num ; ++i) { + rescale_factor_en[i] = ctx->jastrow.rescale_factor_en[i]; + } + + return QMCKL_SUCCESS; +} + #+end_src ** Test #+begin_src c :tangle (eval c_test) @@ -1149,13 +1302,12 @@ 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; +int64_t nucl_num = n2_nucl_num; +double rescale_factor_ee = 1.0; +double rescale_factor_en[2] = { 1.0, 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; @@ -1165,42 +1317,23 @@ qmckl_exit_code rc; assert(!qmckl_electron_provided(context)); -rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num); +rc = qmckl_check(context, + qmckl_set_electron_num (context, elec_up_num, elec_dn_num) + ); assert(rc == QMCKL_SUCCESS); assert(qmckl_electron_provided(context)); -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); - - -rc = qmckl_set_electron_coord (context, 'N', walk_num, elec_coord, walk_num*3*elec_num); +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_get_electron_coord (context, 'N', 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] ); @@ -1211,34 +1344,27 @@ for (int64_t i=0 ; i<3*elec_num ; ++i) { assert(!qmckl_nucleus_provided(context)); -rc = qmckl_set_nucleus_num (context, nucl_num); +rc = qmckl_check(context, + qmckl_set_nucleus_num (context, nucl_num) + ); assert(rc == QMCKL_SUCCESS); assert(!qmckl_nucleus_provided(context)); -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); +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_get_nucleus_coord (context, 'N', nucl_coord2, nucl_num*3); +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 ; idate > ctx->jastrow.asymp_jasb_date) { @@ -1373,7 +1500,7 @@ qmckl_exit_code qmckl_provide_asymp_jasb(qmckl_context context) rc = qmckl_compute_asymp_jasb(context, ctx->jastrow.bord_num, ctx->jastrow.bord_vector, - rescale_factor_kappa_ee, + ctx->jastrow.rescale_factor_ee, ctx->jastrow.asymp_jasb); if (rc != QMCKL_SUCCESS) { return rc; @@ -1394,28 +1521,28 @@ qmckl_exit_code qmckl_provide_asymp_jasb(qmckl_context context) :END: #+NAME: qmckl_asymp_jasb_args - | 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 | + | 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_ee~ | ~double~ | in | Electron coordinates | + | ~asymp_jasb~ | ~double[2]~ | out | Asymptotic value | #+begin_src f90 :comments org :tangle (eval f) :noweb yes -integer function qmckl_compute_asymp_jasb_f(context, bord_num, bord_vector, rescale_factor_kappa_ee, asymp_jasb) & +integer function qmckl_compute_asymp_jasb_f(context, bord_num, bord_vector, rescale_factor_ee, asymp_jasb) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: bord_num double precision , intent(in) :: bord_vector(bord_num + 1) - double precision , intent(in) :: rescale_factor_kappa_ee + double precision , intent(in) :: rescale_factor_ee double precision , intent(out) :: asymp_jasb(2) integer*8 :: i, p double precision :: kappa_inv, x, asym_one - kappa_inv = 1.0d0 / rescale_factor_kappa_ee + kappa_inv = 1.0d0 / rescale_factor_ee info = QMCKL_SUCCESS @@ -1448,7 +1575,7 @@ 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, + const double rescale_factor_ee, double* const asymp_jasb ) { if (context == QMCKL_NULL_CONTEXT){ @@ -1459,7 +1586,7 @@ qmckl_exit_code qmckl_compute_asymp_jasb ( return QMCKL_INVALID_ARG_2; } - const double kappa_inv = 1.0 / rescale_factor_kappa_ee; + const double kappa_inv = 1.0 / rescale_factor_ee; const double asym_one = bord_vector[0] * kappa_inv / (1.0 + bord_vector[1] * kappa_inv); asymp_jasb[0] = asym_one; asymp_jasb[1] = 0.5 * asym_one; @@ -1483,7 +1610,7 @@ 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, + const double rescale_factor_ee, double* const asymp_jasb ); #+end_src @@ -1527,26 +1654,64 @@ 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); +rc = qmckl_check(context, + 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); +rc = qmckl_check(context, + 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); +rc = qmckl_check(context, + 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); +rc = qmckl_check(context, + 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)); +rc = qmckl_check(context, + 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); +rc = qmckl_check(context, + 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); +rc = qmckl_check(context, + qmckl_set_jastrow_cord_vector(context, cord_vector,dim_cord_vect*type_nucl_num) + ); assert(rc == QMCKL_SUCCESS); +double k_ee = 0.; +double k_en[2] = { 0., 0. }; +rc = qmckl_check(context, + qmckl_set_jastrow_rescale_factor_en(context, rescale_factor_en, nucl_num) + ); +assert(rc == QMCKL_SUCCESS); + +rc = qmckl_check(context, + qmckl_set_jastrow_rescale_factor_ee(context, rescale_factor_ee) + ); +assert(rc == QMCKL_SUCCESS); + +rc = qmckl_check(context, + qmckl_get_jastrow_rescale_factor_ee (context, &k_ee) + ); +assert(rc == QMCKL_SUCCESS); +assert(k_ee == rescale_factor_ee); + +rc = qmckl_check(context, + qmckl_get_jastrow_rescale_factor_en (context, &(k_en[0]), nucl_num) + ); +assert(rc == QMCKL_SUCCESS); +for (int i=0 ; idate > ctx->jastrow.factor_ee_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { - free(ctx->jastrow.factor_ee); - ctx->jastrow.factor_ee = NULL; + if (ctx->jastrow.factor_ee != NULL) { + rc = qmckl_free(context, ctx->jastrow.factor_ee); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_provide_factor_ee", + "Unable to free ctx->jastrow.factor_ee"); + } + ctx->jastrow.factor_ee = NULL; + } } /* Allocate array */ @@ -1660,7 +1832,7 @@ qmckl_exit_code qmckl_provide_factor_ee(qmckl_context context) ctx->electron.up_num, ctx->jastrow.bord_num, ctx->jastrow.bord_vector, - ctx->electron.ee_distance_rescaled, + ctx->jastrow.ee_distance_rescaled, ctx->jastrow.asymp_jasb, ctx->jastrow.factor_ee); if (rc != QMCKL_SUCCESS) { @@ -1892,9 +2064,12 @@ print("factor_ee :",factor_ee) assert(qmckl_jastrow_provided(context)); double factor_ee[walk_num]; -rc = qmckl_get_jastrow_factor_ee(context, factor_ee, walk_num); +rc = qmckl_check(context, + qmckl_get_jastrow_factor_ee(context, factor_ee, walk_num) + ); // calculate factor_ee +printf("%e\n%e\n\n",factor_ee[0],-4.282760865958113 ); assert(fabs(factor_ee[0]+4.282760865958113) < 1.e-12); #+end_src @@ -1979,8 +2154,15 @@ qmckl_exit_code qmckl_provide_factor_ee_deriv_e(qmckl_context context) if (ctx->date > ctx->jastrow.factor_ee_deriv_e_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { - free(ctx->jastrow.factor_ee_deriv_e); - ctx->jastrow.factor_ee_deriv_e = NULL; + if (ctx->jastrow.factor_ee_deriv_e != NULL) { + rc = qmckl_free(context, ctx->jastrow.factor_ee_deriv_e); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_provide_factor_ee_deriv_e", + "Unable to free ctx->jastrow.factor_ee_deriv_e"); + } + ctx->jastrow.factor_ee_deriv_e = NULL; + } } /* Allocate array */ @@ -2005,8 +2187,8 @@ qmckl_exit_code qmckl_provide_factor_ee_deriv_e(qmckl_context context) ctx->electron.up_num, ctx->jastrow.bord_num, ctx->jastrow.bord_vector, - ctx->electron.ee_distance_rescaled, - ctx->electron.ee_distance_rescaled_deriv_e, + ctx->jastrow.ee_distance_rescaled, + ctx->jastrow.ee_distance_rescaled_deriv_e, ctx->jastrow.factor_ee_deriv_e); if (rc != QMCKL_SUCCESS) { return rc; @@ -2562,8 +2744,15 @@ qmckl_exit_code qmckl_provide_factor_en(qmckl_context context) if (ctx->date > ctx->jastrow.factor_en_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { - free(ctx->jastrow.factor_en); - ctx->jastrow.factor_en = NULL; + if (ctx->jastrow.factor_en != NULL) { + rc = qmckl_free(context, ctx->jastrow.factor_en); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_provide_factor_en", + "Unable to free ctx->jastrow.factor_en"); + } + ctx->jastrow.factor_en = NULL; + } } /* Allocate array */ @@ -2590,7 +2779,7 @@ qmckl_exit_code qmckl_provide_factor_en(qmckl_context context) ctx->jastrow.type_nucl_vector, ctx->jastrow.aord_num, ctx->jastrow.aord_vector, - ctx->electron.en_distance_rescaled, + ctx->jastrow.en_distance_rescaled, ctx->jastrow.factor_en); if (rc != QMCKL_SUCCESS) { return rc; @@ -2913,8 +3102,15 @@ qmckl_exit_code qmckl_provide_factor_en_deriv_e(qmckl_context context) if (ctx->date > ctx->jastrow.factor_en_deriv_e_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { - free(ctx->jastrow.factor_en_deriv_e); - ctx->jastrow.factor_en_deriv_e = NULL; + if (ctx->jastrow.factor_en_deriv_e != NULL) { + rc = qmckl_free(context, ctx->jastrow.factor_en_deriv_e); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_provide_factor_en_deriv_e", + "Unable to free ctx->jastrow.factor_en_deriv_e"); + } + ctx->jastrow.factor_en_deriv_e = NULL; + } } /* Allocate array */ @@ -2941,8 +3137,8 @@ qmckl_exit_code qmckl_provide_factor_en_deriv_e(qmckl_context context) ctx->jastrow.type_nucl_vector, ctx->jastrow.aord_num, ctx->jastrow.aord_vector, - ctx->electron.en_distance_rescaled, - ctx->electron.en_distance_rescaled_deriv_e, + ctx->jastrow.en_distance_rescaled, + ctx->jastrow.en_distance_rescaled_deriv_e, ctx->jastrow.factor_en_deriv_e); if (rc != QMCKL_SUCCESS) { return rc; @@ -3258,7 +3454,506 @@ assert(fabs(factor_en_deriv_e[0][3][0]+0.9667363412285741 ) < 1.e-12); #+end_src -** Electron-electron rescaled distances for each order +** Electron-electron rescaled distances + + ~ee_distance_rescaled~ stores the matrix of the rescaled distances between all + pairs of electrons: + + \[ + C_{ij} = \left( 1 - \exp{-\kappa C_{ij}}\right)/\kappa + \] + + where \(C_{ij}\) is the matrix of electron-electron distances. + +*** Get + + #+begin_src c :comments org :tangle (eval h_func) :noweb yes +qmckl_exit_code qmckl_get_jastrow_ee_distance_rescaled(qmckl_context context, double* const distance_rescaled); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_get_jastrow_ee_distance_rescaled(qmckl_context context, double* const distance_rescaled) +{ + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + + qmckl_exit_code rc; + + rc = qmckl_provide_ee_distance_rescaled(context); + if (rc != QMCKL_SUCCESS) return rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; + assert (ctx != NULL); + + size_t sze = ctx->electron.num * ctx->electron.num * ctx->electron.walker.num; + memcpy(distance_rescaled, ctx->jastrow.ee_distance_rescaled, sze * sizeof(double)); + + return QMCKL_SUCCESS; +} + #+end_src + +*** Provide :noexport: + + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none +qmckl_exit_code qmckl_provide_ee_distance_rescaled(qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_provide_ee_distance_rescaled(qmckl_context context) +{ + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; + assert (ctx != NULL); + + + /* Compute if necessary */ + if (ctx->electron.walker.point.date > ctx->jastrow.ee_distance_rescaled_date) { + + if (ctx->electron.walker.num > ctx->electron.walker_old.num) { + if (ctx->jastrow.ee_distance_rescaled != NULL) { + qmckl_exit_code rc = qmckl_free(context, ctx->jastrow.ee_distance_rescaled); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_provide_ee_distance_rescaled", + "Unable to free ctx->jastrow.ee_distance_rescaled"); + } + ctx->jastrow.ee_distance_rescaled = NULL; + } + } + + /* Allocate array */ + if (ctx->jastrow.ee_distance_rescaled == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->electron.num * ctx->electron.num * + ctx->electron.walker.num * sizeof(double); + double* ee_distance_rescaled = (double*) qmckl_malloc(context, mem_info); + + if (ee_distance_rescaled == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_provide_ee_distance_rescaled", + NULL); + } + ctx->jastrow.ee_distance_rescaled = ee_distance_rescaled; + } + + qmckl_exit_code rc = + qmckl_compute_ee_distance_rescaled(context, + ctx->electron.num, + ctx->jastrow.rescale_factor_ee, + ctx->electron.walker.num, + ctx->electron.walker.point.coord.data, + ctx->jastrow.ee_distance_rescaled); + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->jastrow.ee_distance_rescaled_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + #+end_src + +*** Compute + :PROPERTIES: + :Name: qmckl_compute_ee_distance_rescaled + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_ee_distance_rescaled_args + | 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 | + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer function qmckl_compute_ee_distance_rescaled_f(context, elec_num, rescale_factor_ee, walk_num, & + coord, ee_distance_rescaled) & + result(info) + use qmckl + implicit none + integer(qmckl_context), intent(in) :: context + integer*8 , intent(in) :: elec_num + double precision , intent(in) :: rescale_factor_ee + integer*8 , intent(in) :: walk_num + double precision , intent(in) :: coord(elec_num,walk_num,3) + double precision , intent(out) :: ee_distance_rescaled(elec_num,elec_num,walk_num) + + integer*8 :: k + + info = QMCKL_SUCCESS + + if (context == QMCKL_NULL_CONTEXT) then + info = QMCKL_INVALID_CONTEXT + return + endif + + if (elec_num <= 0) then + info = QMCKL_INVALID_ARG_2 + return + endif + + if (walk_num <= 0) then + info = QMCKL_INVALID_ARG_3 + return + endif + + do k=1,walk_num + info = qmckl_distance_rescaled(context, 'T', 'T', elec_num, elec_num, & + coord(1,k,1), elec_num * walk_num, & + coord(1,k,1), elec_num * walk_num, & + ee_distance_rescaled(1,1,k), elec_num, rescale_factor_ee) + if (info /= QMCKL_SUCCESS) then + exit + endif + end do + +end function qmckl_compute_ee_distance_rescaled_f + #+end_src + + #+begin_src c :tangle (eval h_private_func) :comments org :exports none +qmckl_exit_code qmckl_compute_ee_distance_rescaled ( + const qmckl_context context, + const int64_t elec_num, + const double rescale_factor_ee, + const int64_t walk_num, + const double* coord, + double* const ee_distance_rescaled ); + #+end_src + + #+CALL: generate_c_interface(table=qmckl_ee_distance_rescaled_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + + #+RESULTS: + #+begin_src f90 :tangle (eval f) :comments org :exports none + integer(c_int32_t) function qmckl_compute_ee_distance_rescaled & + (context, elec_num, rescale_factor_ee, walk_num, coord, ee_distance_rescaled) & + bind(C) result(info) + + use, intrinsic :: iso_c_binding + implicit none + + integer (c_int64_t) , 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,3,walk_num) + real (c_double ) , intent(out) :: ee_distance_rescaled(elec_num,elec_num,walk_num) + + integer(c_int32_t), external :: qmckl_compute_ee_distance_rescaled_f + info = qmckl_compute_ee_distance_rescaled_f & + (context, elec_num, rescale_factor_ee, walk_num, coord, ee_distance_rescaled) + + end function qmckl_compute_ee_distance_rescaled + #+end_src + +*** Test + + #+begin_src python :results output :exports none +import numpy as np + +kappa = 1.0 + +elec_1_w1 = np.array( [-0.250655104764153, 0.503070975550133 , -0.166554344502303]) +elec_2_w1 = np.array( [-0.587812193472177, -0.128751981129274 , 0.187773606533075]) +elec_5_w1 = np.array( [-0.127732483187947, -0.138975497694196 , -8.669850480215846E-002]) +elec_6_w1 = np.array( [-0.232271834949124, -1.059321673434182E-002 , -0.504862241464867]) + +print ( "[0][0] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_1_w1-elec_1_w1)) )/kappa ) +print ( "[0][1] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_1_w1-elec_2_w1)) )/kappa ) +print ( "[1][0] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_2_w1-elec_1_w1)) )/kappa ) +print ( "[5][5] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_5_w1-elec_5_w1)) )/kappa ) +print ( "[5][6] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_5_w1-elec_6_w1)) )/kappa ) +print ( "[6][5] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_6_w1-elec_5_w1)) )/kappa ) + #+end_src + + #+RESULTS: + : [0][0] : 0.0 + : [0][1] : 0.5502278003524018 + : [1][0] : 0.5502278003524018 + : [5][5] : 0.0 + : [5][6] : 0.3622098222364193 + : [6][5] : 0.3622098222364193 + + #+begin_src c :tangle (eval c_test) +assert(qmckl_electron_provided(context)); + + +double ee_distance_rescaled[walk_num * elec_num * elec_num]; +rc = qmckl_get_jastrow_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.5502278003524018) < 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.3622098222364193) < 1.e-12); + + #+end_src + +** Electron-electron rescaled distance gradients and Laplacian with respect to electron coords + + The rescaled distances which is given as $R = (1 - \exp{-\kappa r})/\kappa$ + needs to be perturbed with respect to the electorn coordinates. + This data is stored in the ~ee_distance_rescaled_deriv_e~ tensor. The + The first three elements of this three index tensor ~[4][num][num]~ gives the + derivatives in the x, y, and z directions $dx, dy, dz$ and the last index + gives the Laplacian $\partial x^2 + \partial y^2 + \partial z^2$. + +*** Get + + #+begin_src c :comments org :tangle (eval h_func) :noweb yes +qmckl_exit_code qmckl_get_jastrow_ee_distance_rescaled_deriv_e(qmckl_context context, double* const distance_rescaled_deriv_e); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_get_jastrow_ee_distance_rescaled_deriv_e(qmckl_context context, double* const distance_rescaled_deriv_e) +{ + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + + qmckl_exit_code rc; + + rc = qmckl_provide_ee_distance_rescaled_deriv_e(context); + if (rc != QMCKL_SUCCESS) return rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; + assert (ctx != NULL); + + size_t sze = 4 * ctx->electron.num * ctx->electron.num * ctx->electron.walker.num; + memcpy(distance_rescaled_deriv_e, ctx->jastrow.ee_distance_rescaled_deriv_e, sze * sizeof(double)); + + return QMCKL_SUCCESS; +} + #+end_src + +*** Provide :noexport: + + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none +qmckl_exit_code qmckl_provide_ee_distance_rescaled_deriv_e(qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_provide_ee_distance_rescaled_deriv_e(qmckl_context context) +{ + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; + assert (ctx != NULL); + + + /* Compute if necessary */ + if (ctx->electron.walker.point.date > ctx->jastrow.ee_distance_rescaled_deriv_e_date) { + + if (ctx->electron.walker.num > ctx->electron.walker_old.num) { + if (ctx->jastrow.ee_distance_rescaled_deriv_e != NULL) { + qmckl_exit_code rc = qmckl_free(context, ctx->jastrow.ee_distance_rescaled_deriv_e); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_provide_ee_distance_rescaled_deriv_e", + "Unable to free ctx->jastrow.ee_distance_rescaled_deriv_e"); + } + ctx->jastrow.ee_distance_rescaled_deriv_e = NULL; + } + } + + /* Allocate array */ + if (ctx->jastrow.ee_distance_rescaled_deriv_e == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = 4 * ctx->electron.num * ctx->electron.num * + ctx->electron.walker.num * sizeof(double); + double* ee_distance_rescaled_deriv_e = (double*) qmckl_malloc(context, mem_info); + + if (ee_distance_rescaled_deriv_e == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_provide_ee_distance_rescaled_deriv_e", + NULL); + } + ctx->jastrow.ee_distance_rescaled_deriv_e = ee_distance_rescaled_deriv_e; + } + + qmckl_exit_code rc = + qmckl_compute_ee_distance_rescaled_deriv_e(context, + ctx->electron.num, + ctx->jastrow.rescale_factor_ee, + ctx->electron.walker.num, + ctx->electron.walker.point.coord.data, + ctx->jastrow.ee_distance_rescaled_deriv_e); + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->jastrow.ee_distance_rescaled_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + #+end_src + +*** Compute + :PROPERTIES: + :Name: qmckl_compute_ee_distance_rescaled_deriv_e + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_ee_distance_rescaled_deriv_e_args + | 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_deriv_e~ | ~double[walk_num][4][elec_num][elec_num]~ | out | Electron-electron rescaled distance derivatives | + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer function qmckl_compute_ee_distance_rescaled_deriv_e_f(context, elec_num, rescale_factor_ee, walk_num, & + coord, ee_distance_rescaled_deriv_e) & + result(info) + use qmckl + implicit none + integer(qmckl_context), intent(in) :: context + integer*8 , intent(in) :: elec_num + double precision , intent(in) :: rescale_factor_ee + integer*8 , intent(in) :: walk_num + double precision , intent(in) :: coord(elec_num,walk_num,3) + double precision , intent(out) :: ee_distance_rescaled_deriv_e(4,elec_num,elec_num,walk_num) + + integer*8 :: k + + info = QMCKL_SUCCESS + + if (context == QMCKL_NULL_CONTEXT) then + info = QMCKL_INVALID_CONTEXT + return + endif + + if (elec_num <= 0) then + info = QMCKL_INVALID_ARG_2 + return + endif + + if (walk_num <= 0) then + info = QMCKL_INVALID_ARG_3 + return + endif + + do k=1,walk_num + info = qmckl_distance_rescaled_deriv_e(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_deriv_e(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_deriv_e_f + #+end_src + + #+begin_src c :tangle (eval h_private_func) :comments org :exports none +qmckl_exit_code qmckl_compute_ee_distance_rescaled_deriv_e ( + const qmckl_context context, + const int64_t elec_num, + const double rescale_factor_ee, + const int64_t walk_num, + const double* coord, + double* const ee_distance_rescaled_deriv_e ); + #+end_src + + #+CALL: generate_c_interface(table=qmckl_ee_distance_rescaled_deriv_e_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + + #+RESULTS: + #+begin_src f90 :tangle (eval f) :comments org :exports none + integer(c_int32_t) function qmckl_compute_ee_distance_rescaled_deriv_e & + (context, elec_num, rescale_factor_ee, walk_num, coord, ee_distance_rescaled_deriv_e) & + bind(C) result(info) + + use, intrinsic :: iso_c_binding + implicit none + + integer (c_int64_t) , 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,3,walk_num) + real (c_double ) , intent(out) :: ee_distance_rescaled_deriv_e(4,elec_num,elec_num,walk_num) + + integer(c_int32_t), external :: qmckl_compute_ee_distance_rescaled_deriv_e_f + info = qmckl_compute_ee_distance_rescaled_deriv_e_f & + (context, elec_num, rescale_factor_ee, walk_num, coord, ee_distance_rescaled_deriv_e) + + end function qmckl_compute_ee_distance_rescaled_deriv_e + #+end_src + +*** Test + + #+begin_src python :results output :exports none +import numpy as np + +# TODO + #+end_src + + #+begin_src c :tangle (eval c_test) +assert(qmckl_electron_provided(context)); + + +double ee_distance_rescaled_deriv_e[4 * walk_num * elec_num * elec_num]; +rc = qmckl_get_jastrow_ee_distance_rescaled_deriv_e(context, ee_distance_rescaled_deriv_e); + +// 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); + + #+end_src + +** Electron-electron-nucleus 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~: @@ -3334,8 +4029,15 @@ qmckl_exit_code qmckl_provide_een_rescaled_e(qmckl_context context) if (ctx->date > ctx->jastrow.een_rescaled_e_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { - free(ctx->jastrow.een_rescaled_e); - ctx->jastrow.een_rescaled_e = NULL; + if (ctx->jastrow.een_rescaled_e != NULL) { + rc = qmckl_free(context, ctx->jastrow.een_rescaled_e); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_provide_een_rescaled_e", + "Unable to free ctx->jastrow.een_rescaled_e"); + } + ctx->jastrow.een_rescaled_e = NULL; + } } /* Allocate array */ @@ -3349,7 +4051,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_e(qmckl_context context) if (een_rescaled_e == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, - "qmckl_een_rescaled_e", + "qmckl_provide_een_rescaled_e", NULL); } ctx->jastrow.een_rescaled_e = een_rescaled_e; @@ -3359,7 +4061,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_e(qmckl_context context) ctx->electron.walker.num, ctx->electron.num, ctx->jastrow.cord_num, - ctx->electron.rescale_factor_kappa_ee, + ctx->jastrow.rescale_factor_ee, ctx->electron.ee_distance, ctx->jastrow.een_rescaled_e); if (rc != QMCKL_SUCCESS) { @@ -3381,19 +4083,19 @@ qmckl_exit_code qmckl_provide_een_rescaled_e(qmckl_context context) :END: #+NAME: qmckl_factor_een_rescaled_e_args - | 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 | + | 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 | + | ~een_rescaled_e~ | ~double[walk_num][0:cord_num][elec_num][elec_num]~ | out | Electron-electron rescaled distances | #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_een_rescaled_e_doc_f( & - context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, & + context, walk_num, elec_num, cord_num, rescale_factor_ee, & ee_distance, een_rescaled_e) & result(info) use qmckl @@ -3402,7 +4104,7 @@ integer function qmckl_compute_een_rescaled_e_doc_f( & 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) :: 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 @@ -3444,7 +4146,7 @@ integer function qmckl_compute_een_rescaled_e_doc_f( & 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)) + een_rescaled_e_ij(k, 2) = dexp(-rescale_factor_ee * ee_distance(i, j, nw)) end do end do @@ -3489,7 +4191,7 @@ end function qmckl_compute_een_rescaled_e_doc_f const int64_t walk_num, const int64_t elec_num, const int64_t cord_num, - const double rescale_factor_kappa_ee, + const double rescale_factor_ee, const double* ee_distance, double* const een_rescaled_e ); #+end_src @@ -3499,7 +4201,7 @@ end function qmckl_compute_een_rescaled_e_doc_f #+RESULTS: #+begin_src f90 :tangle (eval f) :comments org :exports none integer(c_int32_t) function qmckl_compute_een_rescaled_e_doc & - (context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, & + (context, walk_num, elec_num, cord_num, rescale_factor_ee, & ee_distance, een_rescaled_e) & bind(C) result(info) @@ -3510,13 +4212,13 @@ end function qmckl_compute_een_rescaled_e_doc_f integer (c_int64_t) , intent(in) , value :: walk_num integer (c_int64_t) , intent(in) , value :: elec_num integer (c_int64_t) , intent(in) , value :: cord_num - real (c_double ) , intent(in) , value :: rescale_factor_kappa_ee + real (c_double ) , intent(in) , value :: rescale_factor_ee real (c_double ) , intent(in) :: ee_distance(elec_num,elec_num,walk_num) real (c_double ) , intent(out) :: een_rescaled_e(elec_num,elec_num,0:cord_num,walk_num) integer(c_int32_t), external :: qmckl_compute_een_rescaled_e_doc_f info = qmckl_compute_een_rescaled_e_doc_f & - (context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, ee_distance, een_rescaled_e) + (context, walk_num, elec_num, cord_num, rescale_factor_ee, ee_distance, een_rescaled_e) end function qmckl_compute_een_rescaled_e_doc #+end_src @@ -3527,7 +4229,7 @@ qmckl_exit_code qmckl_compute_een_rescaled_e_hpc ( const int64_t walk_num, const int64_t elec_num, const int64_t cord_num, - const double rescale_factor_kappa_ee, + const double rescale_factor_ee, const double* ee_distance, double* const een_rescaled_e ) { @@ -3589,8 +4291,8 @@ qmckl_exit_code qmckl_compute_een_rescaled_e_hpc ( k = 0; for (int i = 0; i < elec_num; ++i) { for (int j = 0; j < i; ++j) { - // een_rescaled_e_ij(k, 2) = dexp(-rescale_factor_kappa_ee * ee_distance(i, j, nw)); - een_rescaled_e_ij[k + elec_pairs] = exp(-rescale_factor_kappa_ee * \ + // een_rescaled_e_ij(k, 2) = dexp(-rescale_factor_ee * ee_distance(i, j, nw)); + een_rescaled_e_ij[k + elec_pairs] = exp(-rescale_factor_ee * \ ee_distance[j + i*elec_num + nw*(elec_num*elec_num)]); k = k + 1; } @@ -3648,7 +4350,7 @@ qmckl_exit_code qmckl_compute_een_rescaled_e_hpc ( const int64_t walk_num, const int64_t elec_num, const int64_t cord_num, - const double rescale_factor_kappa_ee, + const double rescale_factor_ee, const double* ee_distance, double* const een_rescaled_e ); #+end_src @@ -3659,7 +4361,7 @@ qmckl_exit_code qmckl_compute_een_rescaled_e_hpc ( const int64_t walk_num, const int64_t elec_num, const int64_t cord_num, - const double rescale_factor_kappa_ee, + const double rescale_factor_ee, const double* ee_distance, double* const een_rescaled_e ); #+end_src @@ -3670,7 +4372,7 @@ qmckl_exit_code qmckl_compute_een_rescaled_e_hpc ( const int64_t walk_num, const int64_t elec_num, const int64_t cord_num, - const double rescale_factor_kappa_ee, + const double rescale_factor_ee, const double* ee_distance, double* const een_rescaled_e ); #+end_src @@ -3681,14 +4383,14 @@ qmckl_exit_code qmckl_compute_een_rescaled_e_hpc ( const int64_t walk_num, const int64_t elec_num, const int64_t cord_num, - const double rescale_factor_kappa_ee, + const double rescale_factor_ee, const double* ee_distance, double* const een_rescaled_e ) { #ifdef HAVE_HPC - return qmckl_compute_een_rescaled_e_hpc(context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, ee_distance, een_rescaled_e); + return qmckl_compute_een_rescaled_e_hpc(context, walk_num, elec_num, cord_num, rescale_factor_ee, ee_distance, een_rescaled_e); #else - return qmckl_compute_een_rescaled_e_doc(context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, ee_distance, een_rescaled_e); + return qmckl_compute_een_rescaled_e_doc(context, walk_num, elec_num, cord_num, rescale_factor_ee, ee_distance, een_rescaled_e); #endif } #+end_src @@ -3770,7 +4472,7 @@ 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); #+end_src -** Electron-electron rescaled distances for each order and derivatives +** Electron-electron-nucleus 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 @@ -3845,8 +4547,15 @@ qmckl_exit_code qmckl_provide_een_rescaled_e_deriv_e(qmckl_context context) if (ctx->date > ctx->jastrow.een_rescaled_e_deriv_e_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { - free(ctx->jastrow.een_rescaled_e_deriv_e); - ctx->jastrow.een_rescaled_e_deriv_e = NULL; + if (ctx->jastrow.een_rescaled_e_deriv_e != NULL) { + rc = qmckl_free(context, ctx->jastrow.een_rescaled_e_deriv_e); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_provide_een_rescaled_e_deriv_e", + "Unable to free ctx->jastrow.een_rescaled_e_deriv_e"); + } + ctx->jastrow.een_rescaled_e_deriv_e = NULL; + } } /* Allocate array */ @@ -3860,7 +4569,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_e_deriv_e(qmckl_context context) if (een_rescaled_e_deriv_e == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, - "qmckl_een_rescaled_e_deriv_e", + "qmckl_provide_een_rescaled_e_deriv_e", NULL); } ctx->jastrow.een_rescaled_e_deriv_e = een_rescaled_e_deriv_e; @@ -3870,7 +4579,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_e_deriv_e(qmckl_context context) ctx->electron.walker.num, ctx->electron.num, ctx->jastrow.cord_num, - ctx->electron.rescale_factor_kappa_ee, + ctx->jastrow.rescale_factor_ee, ctx->electron.walker.point.coord.data, ctx->electron.ee_distance, ctx->jastrow.een_rescaled_e, @@ -3894,21 +4603,21 @@ qmckl_exit_code qmckl_provide_een_rescaled_e_deriv_e(qmckl_context context) :END: #+NAME: qmckl_factor_een_rescaled_e_deriv_e_args - | 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_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_deriv_e~ | ~double[walk_num][0:cord_num][elec_num][4][elec_num]~ | out | Electron-electron rescaled distances | + | 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_deriv_e~ | ~double[walk_num][0:cord_num][elec_num][4][elec_num]~ | out | Electron-electron rescaled distances | #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_factor_een_rescaled_e_deriv_e_f( & - context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, & + context, walk_num, elec_num, cord_num, rescale_factor_ee, & coord_ee, ee_distance, een_rescaled_e, een_rescaled_e_deriv_e) & result(info) use qmckl @@ -3917,7 +4626,7 @@ integer function qmckl_compute_factor_een_rescaled_e_deriv_e_f( & 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) :: 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) @@ -3966,7 +4675,7 @@ integer function qmckl_compute_factor_een_rescaled_e_deriv_e_f( & ! prepare the actual een table do l = 1, cord_num - kappa_l = - dble(l) * rescale_factor_kappa_ee + kappa_l = - dble(l) * rescale_factor_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) @@ -4003,7 +4712,7 @@ end function qmckl_compute_factor_een_rescaled_e_deriv_e_f const int64_t walk_num, const int64_t elec_num, const int64_t cord_num, - const double rescale_factor_kappa_ee, + const double rescale_factor_ee, const double* coord_ee, const double* ee_distance, const double* een_rescaled_e, @@ -4020,7 +4729,7 @@ end function qmckl_compute_factor_een_rescaled_e_deriv_e_f walk_num, & elec_num, & cord_num, & - rescale_factor_kappa_ee, & + rescale_factor_ee, & coord_ee, & ee_distance, & een_rescaled_e, & @@ -4034,7 +4743,7 @@ end function qmckl_compute_factor_een_rescaled_e_deriv_e_f integer (c_int64_t) , intent(in) , value :: walk_num integer (c_int64_t) , intent(in) , value :: elec_num integer (c_int64_t) , intent(in) , value :: cord_num - real (c_double ) , intent(in) , value :: rescale_factor_kappa_ee + real (c_double ) , intent(in) , value :: rescale_factor_ee real (c_double ) , intent(in) :: coord_ee(elec_num,3,walk_num) real (c_double ) , intent(in) :: ee_distance(elec_num,elec_num,walk_num) real (c_double ) , intent(in) :: een_rescaled_e(elec_num,elec_num,0:cord_num,walk_num) @@ -4046,7 +4755,7 @@ end function qmckl_compute_factor_een_rescaled_e_deriv_e_f walk_num, & elec_num, & cord_num, & - rescale_factor_kappa_ee, & + rescale_factor_ee, & coord_ee, & ee_distance, & een_rescaled_e, & @@ -4137,7 +4846,6 @@ for l in range(0,cord_num+1): : een_rescaled_e_deriv_e[2, 1, 6, 2] = 0.5880599146214673 #+begin_src c :tangle (eval c_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, @@ -4152,7 +4860,596 @@ assert(fabs(een_rescaled_e_deriv_e[0][2][1][0][4] + 0.00041342909359870457) < 1. assert(fabs(een_rescaled_e_deriv_e[0][2][1][0][5] + 0.5880599146214673 ) < 1.e-12); #+end_src -** Electron-nucleus rescaled distances for each order +** Electron-nucleus rescaled distances + + ~en_distance_rescaled~ stores the matrix of the rescaled distances between + electrons and nuclei. + + \[ + C_{ij} = \left( 1 - \exp{-\kappa C_{ij}}\right)/\kappa + \] + + where \(C_{ij}\) is the matrix of electron-nucleus distances. + +*** Get + + #+begin_src c :comments org :tangle (eval h_func) :noweb yes +qmckl_exit_code qmckl_get_electron_en_distance_rescaled(qmckl_context context, double* distance_rescaled); + #+end_src + + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_get_electron_en_distance_rescaled(qmckl_context context, double* distance_rescaled) +{ + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + + qmckl_exit_code rc; + + rc = qmckl_provide_en_distance_rescaled(context); + if (rc != QMCKL_SUCCESS) return rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; + assert (ctx != NULL); + + size_t sze = ctx->electron.num * ctx->nucleus.num * ctx->electron.walker.num; + memcpy(distance_rescaled, ctx->jastrow.en_distance_rescaled, sze * sizeof(double)); + + return QMCKL_SUCCESS; +} + #+end_src + +*** Provide :noexport: + + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none +qmckl_exit_code qmckl_provide_en_distance_rescaled(qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_provide_en_distance_rescaled(qmckl_context context) +{ + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; + assert (ctx != NULL); + + if (!(ctx->nucleus.provided)) { + return QMCKL_NOT_PROVIDED; + } + + /* Compute if necessary */ + if (ctx->electron.walker.point.date > ctx->jastrow.en_distance_rescaled_date) { + + if (ctx->electron.walker.num > ctx->electron.walker_old.num) { + if (ctx->jastrow.en_distance_rescaled != NULL) { + qmckl_exit_code rc = qmckl_free(context, ctx->jastrow.en_distance_rescaled); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_provide_en_distance_rescaled", + "Unable to free ctx->jastrow.en_distance_rescaled"); + } + ctx->jastrow.en_distance_rescaled = NULL; + } + } + + /* Allocate array */ + if (ctx->jastrow.en_distance_rescaled == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->electron.num * ctx->nucleus.num * + ctx->electron.walker.num * sizeof(double); + double* en_distance_rescaled = (double*) qmckl_malloc(context, mem_info); + + if (en_distance_rescaled == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_provide_en_distance_rescaled", + NULL); + } + ctx->jastrow.en_distance_rescaled = en_distance_rescaled; + } + + qmckl_exit_code rc = + qmckl_compute_en_distance_rescaled(context, + ctx->electron.num, + ctx->nucleus.num, + ctx->jastrow.rescale_factor_en, + ctx->electron.walker.num, + ctx->electron.walker.point.coord.data, + ctx->nucleus.coord.data, + ctx->jastrow.en_distance_rescaled); + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->jastrow.en_distance_rescaled_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + #+end_src + +*** Compute + :PROPERTIES: + :Name: qmckl_compute_en_distance_rescaled + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_en_distance_rescaled_args + | 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 | + | ~rescale_factor_en~ | ~double[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 | + | | | | | + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer function qmckl_compute_en_distance_rescaled_f(context, elec_num, nucl_num, rescale_factor_en, walk_num, elec_coord, & + nucl_coord, en_distance_rescaled) & + result(info) + use qmckl + implicit none + integer(qmckl_context), intent(in) :: context + integer*8 , intent(in) :: elec_num + integer*8 , intent(in) :: nucl_num + 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(elec_num,nucl_num,walk_num) + + integer*8 :: i, k + double precision :: coord(3) + + info = QMCKL_SUCCESS + + if (context == QMCKL_NULL_CONTEXT) then + info = QMCKL_INVALID_CONTEXT + return + endif + + if (elec_num <= 0) then + info = QMCKL_INVALID_ARG_2 + return + endif + + if (nucl_num <= 0) then + info = QMCKL_INVALID_ARG_3 + return + endif + + if (walk_num <= 0) then + info = QMCKL_INVALID_ARG_5 + return + endif + + do i=1, nucl_num + coord(1:3) = nucl_coord(i,1:3) + do k=1,walk_num + info = qmckl_distance_rescaled(context, 'T', 'T', elec_num, 1_8, & + elec_coord(1,k,1), elec_num*walk_num, coord, 1_8, & + en_distance_rescaled(1,i,k), elec_num, rescale_factor_en(i)) + if (info /= QMCKL_SUCCESS) then + return + endif + end do + end do + +end function qmckl_compute_en_distance_rescaled_f + #+end_src + + #+begin_src c :tangle (eval h_private_func) :comments org :exports none +qmckl_exit_code qmckl_compute_en_distance_rescaled ( + const qmckl_context context, + const int64_t elec_num, + const int64_t nucl_num, + const double* rescale_factor_en, + const int64_t walk_num, + const double* elec_coord, + const double* nucl_coord, + double* const en_distance_rescaled ); + #+end_src + + #+CALL: generate_c_interface(table=qmckl_en_distance_rescaled_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + + #+RESULTS: + #+begin_src f90 :tangle (eval f) :comments org :exports none + integer(c_int32_t) function qmckl_compute_en_distance_rescaled & + (context, & + elec_num, & + nucl_num, & + rescale_factor_en, & + walk_num, & + elec_coord, & + nucl_coord, & + en_distance_rescaled) & + bind(C) result(info) + + use, intrinsic :: iso_c_binding + implicit none + + integer (c_int64_t) , intent(in) , value :: context + integer (c_int64_t) , intent(in) , value :: elec_num + integer (c_int64_t) , intent(in) , value :: nucl_num + real (c_double ) , intent(in) :: rescale_factor_en(nucl_num) + integer (c_int64_t) , intent(in) , value :: walk_num + real (c_double ) , intent(in) :: elec_coord(elec_num,walk_num,3) + real (c_double ) , intent(in) :: nucl_coord(elec_num,3) + real (c_double ) , intent(out) :: en_distance_rescaled(elec_num,nucl_num,walk_num) + + integer(c_int32_t), external :: qmckl_compute_en_distance_rescaled_f + info = qmckl_compute_en_distance_rescaled_f & + (context, & + elec_num, & + nucl_num, & + rescale_factor_en, & + walk_num, & + elec_coord, & + nucl_coord, & + en_distance_rescaled) + + end function qmckl_compute_en_distance_rescaled + #+end_src + +*** Test + + #+begin_src python :results output :exports none +import numpy as np + +kappa = 1.0 + +elec_1_w1 = np.array( [-0.250655104764153, 0.503070975550133 , -0.166554344502303]) +elec_2_w1 = np.array( [-0.587812193472177, -0.128751981129274 , 0.187773606533075]) +elec_5_w1 = np.array( [-0.127732483187947, -0.138975497694196 , -8.669850480215846E-002]) +elec_6_w1 = np.array( [-0.232271834949124, -1.059321673434182E-002 , -0.504862241464867]) +nucl_1 = np.array( [ 0., 0., 0. ]) +nucl_2 = np.array( [ 0., 0., 2.059801 ]) + +print ( "[0][0] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_1_w1-nucl_1)) )/kappa ) +print ( "[1][0] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_1_w1-nucl_2)) )/kappa ) +print ( "[0][1] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_2_w1-nucl_1)) )/kappa ) +print ( "[0][5] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_5_w1-nucl_1)) )/kappa ) +print ( "[1][5] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_5_w1-nucl_2)) )/kappa ) +print ( "[0][6] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_6_w1-nucl_1)) )/kappa ) + + #+end_src + + #+RESULTS: + : [0][0] : 0.4435709484118112 + : [1][0] : 0.8993601506374442 + : [0][1] : 0.46760219699910477 + : [0][5] : 0.1875631834682101 + : [1][5] : 0.8840716589810682 + : [0][6] : 0.42640469987268914 + + #+begin_src c :tangle (eval c_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.4435709484118112) < 1.e-12); + +// (1,2,1) +assert(fabs(en_distance_rescaled[0][1][0] - 0.8993601506374442) < 1.e-12); + +// (2,1,1) +assert(fabs(en_distance_rescaled[0][0][1] - 0.46760219699910477) < 1.e-12); + +// (1,1,2) +assert(fabs(en_distance_rescaled[0][0][5] - 0.1875631834682101) < 1.e-12); + +// (1,2,2) +assert(fabs(en_distance_rescaled[0][1][5] - 0.8840716589810682) < 1.e-12); + +// (2,1,2) +assert(fabs(en_distance_rescaled[0][0][6] - 0.42640469987268914) < 1.e-12); + + #+end_src + +** Electron-electron-nucleus rescaled distance gradients and laplacian with respect to electron coords + + The rescaled distances which is given as $R = (1 - \exp{-\kappa r})/\kappa$ + needs to be perturbed with respect to the nuclear coordinates. + This data is stored in the ~en_distance_rescaled_deriv_e~ tensor. The + The first three elements of this three index tensor ~[4][nucl_num][elec_num]~ gives the + derivatives in the x, y, and z directions $dx, dy, dz$ and the last index + gives the Laplacian $\partial x^2 + \partial y^2 + \partial z^2$. + +*** Get + + #+begin_src c :comments org :tangle (eval h_func) :noweb yes +qmckl_exit_code qmckl_get_electron_en_distance_rescaled_deriv_e(qmckl_context context, double* distance_rescaled_deriv_e); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_get_electron_en_distance_rescaled_deriv_e(qmckl_context context, double* distance_rescaled_deriv_e) +{ + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + + qmckl_exit_code rc; + + rc = qmckl_provide_en_distance_rescaled_deriv_e(context); + if (rc != QMCKL_SUCCESS) return rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; + assert (ctx != NULL); + + size_t sze = 4 * ctx->electron.num * ctx->nucleus.num * ctx->electron.walker.num; + memcpy(distance_rescaled_deriv_e, ctx->jastrow.en_distance_rescaled_deriv_e, sze * sizeof(double)); + + return QMCKL_SUCCESS; +} + #+end_src + +*** Provide :noexport: + + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none +qmckl_exit_code qmckl_provide_en_distance_rescaled_deriv_e(qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_provide_en_distance_rescaled_deriv_e(qmckl_context context) +{ + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; + assert (ctx != NULL); + + if (!(ctx->nucleus.provided)) { + return QMCKL_NOT_PROVIDED; + } + + /* Compute if necessary */ + if (ctx->electron.walker.point.date > ctx->jastrow.en_distance_rescaled_deriv_e_date) { + + if (ctx->electron.walker.num > ctx->electron.walker_old.num) { + if (ctx->jastrow.en_distance_rescaled_deriv_e != NULL) { + qmckl_exit_code rc = qmckl_free(context, ctx->jastrow.en_distance_rescaled_deriv_e); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_provide_en_distance_rescaled_deriv_e", + "Unable to free ctx->jastrow.en_distance_rescaled_deriv_e"); + } + ctx->jastrow.en_distance_rescaled_deriv_e = NULL; + } + } + + /* Allocate array */ + if (ctx->jastrow.en_distance_rescaled_deriv_e == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = 4 * ctx->electron.num * ctx->nucleus.num * + ctx->electron.walker.num * sizeof(double); + double* en_distance_rescaled_deriv_e = (double*) qmckl_malloc(context, mem_info); + + if (en_distance_rescaled_deriv_e == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_provide_en_distance_rescaled_deriv_e", + NULL); + } + ctx->jastrow.en_distance_rescaled_deriv_e = en_distance_rescaled_deriv_e; + } + + qmckl_exit_code rc = + qmckl_compute_en_distance_rescaled_deriv_e(context, + ctx->electron.num, + ctx->nucleus.num, + ctx->jastrow.rescale_factor_en, + ctx->electron.walker.num, + ctx->electron.walker.point.coord.data, + ctx->nucleus.coord.data, + ctx->jastrow.en_distance_rescaled_deriv_e); + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->jastrow.en_distance_rescaled_deriv_e_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + #+end_src + +*** Compute + :PROPERTIES: + :Name: qmckl_compute_en_distance_rescaled_deriv_e + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_en_distance_rescaled_deriv_e_args + | 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 | + | ~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_deriv_e~ | ~double[walk_num][nucl_num][elec_num][4]~ | out | Electron-nucleus distance derivatives | + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer function qmckl_compute_en_distance_rescaled_deriv_e_f(context, elec_num, nucl_num, & + rescale_factor_en, walk_num, elec_coord, & + nucl_coord, en_distance_rescaled_deriv_e) & + result(info) + use qmckl + implicit none + integer(qmckl_context), intent(in) :: context + integer*8 , intent(in) :: elec_num + integer*8 , intent(in) :: 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_deriv_e(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_deriv_e(context, 'T', 'T', elec_num, 1_8, & + elec_coord(1,k,1), elec_num*walk_num, coord, 1_8, & + en_distance_rescaled_deriv_e(1,1,i,k), elec_num, rescale_factor_en(i)) + if (info /= QMCKL_SUCCESS) then + return + endif + end do + end do + +end function qmckl_compute_en_distance_rescaled_deriv_e_f + #+end_src + + #+begin_src c :tangle (eval h_private_func) :comments org :exports none +qmckl_exit_code qmckl_compute_en_distance_rescaled_deriv_e ( + const qmckl_context context, + const int64_t elec_num, + const int64_t nucl_num, + const double* rescale_factor_en, + const int64_t walk_num, + const double* elec_coord, + const double* nucl_coord, + double* const en_distance_rescaled_deriv_e ); + #+end_src + + #+CALL: generate_c_interface(table=qmckl_en_distance_rescaled_deriv_e_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + + #+RESULTS: + #+begin_src f90 :tangle (eval f) :comments org :exports none + integer(c_int32_t) function qmckl_compute_en_distance_rescaled_deriv_e & + (context, & + elec_num, & + nucl_num, & + rescale_factor_en, & + walk_num, & + elec_coord, & + nucl_coord, & + en_distance_rescaled_deriv_e) & + bind(C) result(info) + + use, intrinsic :: iso_c_binding + implicit none + + integer (c_int64_t) , intent(in) , value :: context + integer (c_int64_t) , intent(in) , value :: elec_num + integer (c_int64_t) , intent(in) , value :: nucl_num + real (c_double ) , intent(in) :: rescale_factor_en(nucl_num) + integer (c_int64_t) , intent(in) , value :: walk_num + real (c_double ) , intent(in) :: elec_coord(elec_num,walk_num,3) + real (c_double ) , intent(in) :: nucl_coord(elec_num,3) + real (c_double ) , intent(out) :: en_distance_rescaled_deriv_e(4,elec_num,nucl_num,walk_num) + + integer(c_int32_t), external :: qmckl_compute_en_distance_rescaled_deriv_e_f + info = qmckl_compute_en_distance_rescaled_deriv_e_f & + (context, & + elec_num, & + nucl_num, & + rescale_factor_en, & + walk_num, & + elec_coord, & + nucl_coord, & + en_distance_rescaled_deriv_e) + + end function qmckl_compute_en_distance_rescaled_deriv_e + #+end_src + +*** Test + + #+begin_src python :results output :exports none +import numpy as np + +# TODO + #+end_src + + + #+begin_src c :tangle (eval c_test) + +assert(qmckl_electron_provided(context)); + +assert(qmckl_nucleus_provided(context)); + +double en_distance_rescaled_deriv_e[walk_num][4][nucl_num][elec_num]; + +rc = qmckl_check(context, + qmckl_get_electron_en_distance_rescaled_deriv_e(context, &(en_distance_rescaled_deriv_e[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); + + #+end_src + +** Electron-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~: @@ -4228,8 +5525,15 @@ qmckl_exit_code qmckl_provide_een_rescaled_n(qmckl_context context) if (ctx->date > ctx->jastrow.een_rescaled_n_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { - free(ctx->jastrow.een_rescaled_n); - ctx->jastrow.een_rescaled_n = NULL; + if (ctx->jastrow.een_rescaled_n != NULL) { + rc = qmckl_free(context, ctx->jastrow.een_rescaled_n); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_provide_een_rescaled_n", + "Unable to free ctx->jastrow.een_rescaled_n"); + } + ctx->jastrow.een_rescaled_n = NULL; + } } /* Allocate array */ @@ -4243,7 +5547,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_n(qmckl_context context) if (een_rescaled_n == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, - "qmckl_een_rescaled_n", + "qmckl_provide_een_rescaled_n", NULL); } ctx->jastrow.een_rescaled_n = een_rescaled_n; @@ -4254,7 +5558,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_n(qmckl_context context) ctx->electron.num, ctx->nucleus.num, ctx->jastrow.cord_num, - ctx->electron.rescale_factor_kappa_en, + ctx->jastrow.rescale_factor_en, ctx->electron.en_distance, ctx->jastrow.een_rescaled_n); if (rc != QMCKL_SUCCESS) { @@ -4276,20 +5580,20 @@ qmckl_exit_code qmckl_provide_een_rescaled_n(qmckl_context context) :END: #+NAME: qmckl_factor_een_rescaled_n_args - | 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 | + | 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_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 | #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_een_rescaled_n_f( & - context, walk_num, elec_num, nucl_num, cord_num, rescale_factor_kappa_en, & + context, walk_num, elec_num, nucl_num, cord_num, rescale_factor_en, & en_distance, een_rescaled_n) & result(info) use qmckl @@ -4299,7 +5603,7 @@ integer function qmckl_compute_een_rescaled_n_f( & 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) :: rescale_factor_en(nucl_num) 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 @@ -4341,7 +5645,7 @@ integer function qmckl_compute_een_rescaled_n_f( & 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)) + een_rescaled_n(i, a, 1, nw) = dexp(-rescale_factor_en(a) * en_distance(i, a, nw)) end do end do @@ -4364,7 +5668,7 @@ qmckl_exit_code qmckl_compute_een_rescaled_n ( const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, - const double rescale_factor_kappa_en, + const double* rescale_factor_en, const double* en_distance, double* const een_rescaled_n ) { @@ -4400,8 +5704,8 @@ qmckl_exit_code qmckl_compute_een_rescaled_n ( // prepare the actual een table //een_rescaled_n(:, :, 0, nw) = 1.0d0 een_rescaled_n[i + a * elec_num + 0 + nw * elec_num*nucl_num*(cord_num+1)] = 1.0; - //een_rescaled_n(i, a, 1, nw) = dexp(-rescale_factor_kappa_en * en_distance(i, a, nw)) - een_rescaled_n[i + a*elec_num + elec_num*nucl_num + nw*elec_num*nucl_num*(cord_num+1)] = exp(-rescale_factor_kappa_en * \ + //een_rescaled_n(i, a, 1, nw) = dexp(-rescale_factor_en(a) * en_distance(i, a, nw)) + een_rescaled_n[i + a*elec_num + elec_num*nucl_num + nw*elec_num*nucl_num*(cord_num+1)] = exp(-rescale_factor_en[a] * \ en_distance[i + a*elec_num + nw*elec_num*nucl_num]); } } @@ -4430,7 +5734,7 @@ qmckl_exit_code qmckl_compute_een_rescaled_n ( const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, - const double rescale_factor_kappa_en, + const double* rescale_factor_en, const double* en_distance, double* const een_rescaled_n ); #+end_src @@ -4570,8 +5874,15 @@ qmckl_exit_code qmckl_provide_een_rescaled_n_deriv_e(qmckl_context context) if (ctx->date > ctx->jastrow.een_rescaled_n_deriv_e_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { - free(ctx->jastrow.een_rescaled_n_deriv_e); - ctx->jastrow.een_rescaled_n_deriv_e = NULL; + if (ctx->jastrow.een_rescaled_n_deriv_e != NULL) { + rc = qmckl_free(context, ctx->jastrow.een_rescaled_n_deriv_e); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_provide_een_rescaled_n_deriv_e", + "Unable to free ctx->jastrow.een_rescaled_n_deriv_e"); + } + ctx->jastrow.een_rescaled_n_deriv_e = NULL; + } } /* Allocate array */ @@ -4585,7 +5896,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_n_deriv_e(qmckl_context context) if (een_rescaled_n_deriv_e == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, - "qmckl_een_rescaled_n_deriv_e", + "qmckl_provide_een_rescaled_n_deriv_e", NULL); } ctx->jastrow.een_rescaled_n_deriv_e = een_rescaled_n_deriv_e; @@ -4596,7 +5907,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_n_deriv_e(qmckl_context context) ctx->electron.num, ctx->nucleus.num, ctx->jastrow.cord_num, - ctx->electron.rescale_factor_kappa_en, + ctx->jastrow.rescale_factor_en, ctx->electron.walker.point.coord.data, ctx->nucleus.coord.data, ctx->electron.en_distance, @@ -4621,24 +5932,24 @@ qmckl_exit_code qmckl_provide_een_rescaled_n_deriv_e(qmckl_context context) :END: #+NAME: qmckl_compute_factor_een_rescaled_n_deriv_e_args - | 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_ee~ | ~double[walk_num][3][elec_num]~ | in | Electron coordinates | - | ~coord_en~ | ~double[3][nucl_num]~ | in | Nuclear coordinates | - | ~en_distance~ | ~double[walk_num][elec_num][nucl_num]~ | in | Electron-nucleus distances | - | ~een_rescaled_n~ | ~double[walk_num][0:cord_num][nucl_num][elec_num]~ | in | Electron-nucleus distances | - | ~een_rescaled_n_deriv_e~ | ~double[walk_num][0:cord_num][nucl_num][4][elec_num]~ | out | Electron-nucleus rescaled distances | + | 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_en~ | ~double[nucl_num]~ | in | Factor to rescale ee distances | + | ~coord_ee~ | ~double[walk_num][3][elec_num]~ | in | Electron coordinates | + | ~coord_en~ | ~double[3][nucl_num]~ | in | Nuclear coordinates | + | ~en_distance~ | ~double[walk_num][elec_num][nucl_num]~ | in | Electron-nucleus distances | + | ~een_rescaled_n~ | ~double[walk_num][0:cord_num][nucl_num][elec_num]~ | in | Electron-nucleus distances | + | ~een_rescaled_n_deriv_e~ | ~double[walk_num][0:cord_num][nucl_num][4][elec_num]~ | out | Electron-nucleus rescaled distances | #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_factor_een_rescaled_n_deriv_e_f( & context, walk_num, elec_num, nucl_num, & - cord_num, rescale_factor_kappa_en, & + cord_num, rescale_factor_en, & coord_ee, coord_en, en_distance, een_rescaled_n, een_rescaled_n_deriv_e) & result(info) use qmckl @@ -4648,7 +5959,7 @@ integer function qmckl_compute_factor_een_rescaled_n_deriv_e_f( & 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) :: rescale_factor_en(nucl_num) double precision , intent(in) :: coord_ee(elec_num,3,walk_num) double precision , intent(in) :: coord_en(nucl_num,3) double precision , intent(in) :: en_distance(elec_num,nucl_num,walk_num) @@ -4703,8 +6014,8 @@ integer function qmckl_compute_factor_een_rescaled_n_deriv_e_f( & end do do l = 0, cord_num - kappa_l = - dble(l) * rescale_factor_kappa_en do a = 1, nucl_num + kappa_l = - dble(l) * rescale_factor_en(a) 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) @@ -4741,7 +6052,7 @@ end function qmckl_compute_factor_een_rescaled_n_deriv_e_f const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, - const double rescale_factor_kappa_en, + const double* rescale_factor_en, const double* coord_ee, const double* coord_en, const double* en_distance, @@ -4759,7 +6070,7 @@ end function qmckl_compute_factor_een_rescaled_n_deriv_e_f elec_num, & nucl_num, & cord_num, & - rescale_factor_kappa_en, & + rescale_factor_en, & coord_ee, & coord_en, & en_distance, & @@ -4775,7 +6086,7 @@ end function qmckl_compute_factor_een_rescaled_n_deriv_e_f integer (c_int64_t) , intent(in) , value :: elec_num integer (c_int64_t) , intent(in) , value :: nucl_num integer (c_int64_t) , intent(in) , value :: cord_num - real (c_double ) , intent(in) , value :: rescale_factor_kappa_en + real (c_double ) , intent(in) :: rescale_factor_en(nucl_num) real (c_double ) , intent(in) :: coord_ee(elec_num,3,walk_num) real (c_double ) , intent(in) :: coord_en(nucl_num,3) real (c_double ) , intent(in) :: en_distance(nucl_num,elec_num,walk_num) @@ -4789,7 +6100,7 @@ end function qmckl_compute_factor_een_rescaled_n_deriv_e_f elec_num, & nucl_num, & cord_num, & - rescale_factor_kappa_en, & + rescale_factor_en, & coord_ee, & coord_en, & en_distance, & @@ -5172,8 +6483,15 @@ qmckl_exit_code qmckl_provide_tmp_c(qmckl_context context) if (ctx->date > ctx->jastrow.tmp_c_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { - free(ctx->jastrow.tmp_c); - ctx->jastrow.tmp_c = NULL; + if (ctx->jastrow.tmp_c != NULL) { + rc = qmckl_free(context, ctx->jastrow.tmp_c); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_provide_tmp_c", + "Unable to free ctx->jastrow.tmp_c"); + } + ctx->jastrow.tmp_c = NULL; + } } /* Allocate array */ @@ -5267,8 +6585,15 @@ qmckl_exit_code qmckl_provide_dtmp_c(qmckl_context context) if (ctx->date > ctx->jastrow.dtmp_c_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { - free(ctx->jastrow.dtmp_c); - ctx->jastrow.dtmp_c = NULL; + if (ctx->jastrow.dtmp_c != NULL) { + rc = qmckl_free(context, ctx->jastrow.dtmp_c); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_provide_dtmp_c", + "Unable to free ctx->jastrow.dtmp_c"); + } + ctx->jastrow.dtmp_c = NULL; + } } /* Allocate array */ @@ -7037,8 +8362,15 @@ qmckl_exit_code qmckl_provide_factor_een(qmckl_context context) if (ctx->date > ctx->jastrow.factor_een_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { - free(ctx->jastrow.factor_een); - ctx->jastrow.factor_een = NULL; + if (ctx->jastrow.factor_een != NULL) { + rc = qmckl_free(context, ctx->jastrow.factor_een); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_provide_factor_een", + "Unable to free ctx->jastrow.factor_een"); + } + ctx->jastrow.factor_een = NULL; + } } /* Allocate array */ @@ -7552,8 +8884,15 @@ qmckl_exit_code qmckl_provide_factor_een_deriv_e(qmckl_context context) if (ctx->date > ctx->jastrow.factor_een_deriv_e_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { - free(ctx->jastrow.factor_een_deriv_e); - ctx->jastrow.factor_een_deriv_e = NULL; + if (ctx->jastrow.factor_een_deriv_e != NULL) { + rc = qmckl_free(context, ctx->jastrow.factor_een_deriv_e); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_provide_factor_een_deriv_e", + "Unable to free ctx->jastrow.factor_een_deriv_e"); + } + ctx->jastrow.factor_een_deriv_e = NULL; + } } /* Allocate array */ diff --git a/org/qmckl_nucleus.org b/org/qmckl_nucleus.org index 159cba1..27e3a8b 100644 --- a/org/qmckl_nucleus.org +++ b/org/qmckl_nucleus.org @@ -74,15 +74,12 @@ int main() { | ~charge~ | qmckl_vector | Nuclear charges | | ~coord~ | qmckl_matrix | Nuclear coordinates, in transposed format | | ~coord_date~ | int64_t | Nuclear coordinates, date if modified | - | ~rescale_factor_kappa~ | double | The distance scaling factor | Computed data: |-----------------------------+--------------+------------------------------------------------------------| | ~nn_distance~ | qmckl_matrix | Nucleus-nucleus distances | | ~nn_distance_date~ | int64_t | Date when Nucleus-nucleus distances were computed | - | ~nn_distance_rescaled~ | qmckl_matrix | Nucleus-nucleus rescaled distances | - | ~nn_distance_rescaled_date~ | int64_t | Date when Nucleus-nucleus rescaled distances were computed | | ~repulsion~ | double | Nuclear repulsion energy | | ~repulsion_date~ | int64_t | Date when the nuclear repulsion energy was computed | @@ -93,14 +90,11 @@ typedef struct qmckl_nucleus_struct { int64_t num; int64_t repulsion_date; int64_t nn_distance_date; - int64_t nn_distance_rescaled_date; int64_t coord_date; qmckl_vector charge; qmckl_matrix coord; qmckl_matrix nn_distance; - qmckl_matrix nn_distance_rescaled; double repulsion; - double rescale_factor_kappa; int32_t uninitialized; bool provided; } qmckl_nucleus_struct; @@ -131,7 +125,6 @@ qmckl_exit_code qmckl_init_nucleus(qmckl_context context) { ctx->nucleus.uninitialized = (1 << 3) - 1; /* Default values */ - ctx->nucleus.rescale_factor_kappa = 1.0; return QMCKL_SUCCESS; } @@ -265,53 +258,6 @@ end interface #+begin_src c :comments org :tangle (eval h_func) :exports none qmckl_exit_code -qmckl_get_nucleus_rescale_factor(const qmckl_context context, - double* const rescale_factor_kappa); - #+end_src - - #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -qmckl_exit_code -qmckl_get_nucleus_rescale_factor (const qmckl_context context, - double* const rescale_factor_kappa) -{ - - if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { - return QMCKL_INVALID_CONTEXT; - } - - if (rescale_factor_kappa == NULL) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_2, - "qmckl_get_nucleus_rescale_factor", - "rescale_factor_kappa is a null pointer"); - } - - qmckl_context_struct* const ctx = (qmckl_context_struct*) context; - assert (ctx != NULL); - - assert (ctx->nucleus.rescale_factor_kappa > 0.0); - - (*rescale_factor_kappa) = ctx->nucleus.rescale_factor_kappa; - - return QMCKL_SUCCESS; -} - #+end_src - - #+begin_src f90 :tangle (eval fh_func) :comments org :exports none -interface - integer(c_int32_t) function qmckl_get_nucleus_rescale_factor(context, kappa) & - bind(C) - use, intrinsic :: iso_c_binding - import - implicit none - integer (c_int64_t) , intent(in) , value :: context - real (c_double) , intent(out) :: kappa - end function qmckl_get_nucleus_rescale_factor -end interface - #+end_src - - #+begin_src c :comments org :tangle (eval h_func) :exports none -qmckl_exit_code qmckl_get_nucleus_coord(const qmckl_context context, const char transp, double* const coord, @@ -631,55 +577,12 @@ interface end interface #+end_src - #+begin_src c :comments org :tangle (eval h_func) -qmckl_exit_code -qmckl_set_nucleus_rescale_factor(qmckl_context context, - const double kappa); - #+end_src - - Sets the rescale parameter for the nuclear distances. - - #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -qmckl_exit_code -qmckl_set_nucleus_rescale_factor(qmckl_context context, - const double rescale_factor_kappa) -{ - int32_t mask = 0; // Can be updated - - <> - - if (rescale_factor_kappa <= 0.0) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_2, - "qmckl_set_nucleus_rescale_factor", - "rescale_factor_kappa cannot be <= 0."); - } - - ctx->nucleus.rescale_factor_kappa = rescale_factor_kappa; - - return QMCKL_SUCCESS; -} - #+end_src - - #+begin_src f90 :tangle (eval fh_func) :comments org :exports none -interface - integer(c_int32_t) function qmckl_set_nucleus_rescale_factor(context, kappa) & - bind(C) - use, intrinsic :: iso_c_binding - import - implicit none - integer (c_int64_t) , intent(in) , value :: context - real (c_double) , intent(in) , value :: kappa - end function qmckl_set_nucleus_rescale_factor -end interface - #+end_src ** Test #+begin_src c :tangle (eval c_test) const double* nucl_charge = chbrclf_charge; const double* nucl_coord = &(chbrclf_nucl_coord[0][0]); -const double nucl_rescale_factor_kappa = 2.0; /* --- */ @@ -693,38 +596,25 @@ assert(rc == QMCKL_NOT_PROVIDED); rc = qmckl_set_nucleus_num (context, chbrclf_nucl_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(!qmckl_nucleus_provided(context)); rc = qmckl_get_nucleus_num (context, &n); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(n == chbrclf_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*chbrclf_nucl_num]; rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2, 3*chbrclf_nucl_num); assert(rc == QMCKL_NOT_PROVIDED); rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), 3*chbrclf_nucl_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(!qmckl_nucleus_provided(context)); rc = qmckl_get_nucleus_coord (context, 'N', nucl_coord2, 3*chbrclf_nucl_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); for (size_t k=0 ; k<3 ; ++k) { for (int64_t i=0 ; inucleus.num * ctx->nucleus.num; - if (sze > size_max) { - return qmckl_failwith(context, - QMCKL_INVALID_ARG_3, - "qmckl_get_nucleus_nn_distance_rescaled", - "Array too small"); - } - rc = qmckl_double_of_matrix(context, - ctx->nucleus.nn_distance_rescaled, - distance_rescaled, - size_max); - return rc; -} - #+end_src - - #+begin_src f90 :tangle (eval fh_func) :comments org :exports none -interface - integer(c_int32_t) function qmckl_get_nucleus_nn_distance_rescaled(context, distance_rescaled, size_max) & - bind(C) - use, intrinsic :: iso_c_binding - import - implicit none - integer (c_int64_t) , intent(in) , value :: context - real (c_double ) , intent(out) :: distance_rescaled(*) - integer (c_int64_t) , intent(in) , value :: size_max - end function -end interface - #+end_src - -*** Provide :noexport: - - #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none -qmckl_exit_code qmckl_provide_nn_distance_rescaled(qmckl_context context); - #+end_src - - #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -qmckl_exit_code qmckl_provide_nn_distance_rescaled(qmckl_context context) -{ - /* Check input parameters */ - if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { - return (char) 0; - } - - qmckl_context_struct* const ctx = (qmckl_context_struct*) context; - assert (ctx != NULL); - - if (!ctx->nucleus.provided) return QMCKL_NOT_PROVIDED; - - /* Allocate array */ - if (ctx->nucleus.nn_distance_rescaled.data == NULL) { - ctx->nucleus.nn_distance_rescaled = - qmckl_matrix_alloc(context, ctx->nucleus.num, ctx->nucleus.num); - - if (ctx->nucleus.nn_distance_rescaled.data == NULL) { - return qmckl_failwith( context, - QMCKL_ALLOCATION_FAILED, - "qmckl_nn_distance_rescaled", - NULL); - } - } - - qmckl_exit_code rc = - qmckl_compute_nn_distance_rescaled(context, - ctx->nucleus.num, - ctx->nucleus.rescale_factor_kappa, - ctx->nucleus.coord.data, - ctx->nucleus.nn_distance_rescaled.data); - if (rc != QMCKL_SUCCESS) { - return rc; - } - - ctx->nucleus.nn_distance_rescaled_date = ctx->date; - - return QMCKL_SUCCESS; -} - #+end_src - -*** Compute - - #+NAME: qmckl_nn_distance_rescaled_args - | qmckl_context | context | in | Global state | - | int64_t | nucl_num | in | Number of nuclei | - | double | coord[3][nucl_num] | in | Nuclear coordinates (au) | - | double | nn_distance_rescaled[nucl_num][nucl_num] | out | Nucleus-nucleus rescaled distances (au) | - - #+begin_src f90 :comments org :tangle (eval f) :noweb yes -integer function qmckl_compute_nn_distance_rescaled_f(context, nucl_num, rescale_factor_kappa, coord, nn_distance_rescaled) & - result(info) - use qmckl - implicit none - integer(qmckl_context), intent(in) :: context - integer*8 , intent(in) :: nucl_num - double precision , intent(in) :: rescale_factor_kappa - double precision , intent(in) :: coord(nucl_num,3) - double precision , intent(out) :: nn_distance_rescaled(nucl_num,nucl_num) - - integer*8 :: k - - 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 - - info = qmckl_distance_rescaled(context, 'T', 'T', nucl_num, nucl_num, & - coord, nucl_num, & - coord, nucl_num, & - nn_distance_rescaled, nucl_num, rescale_factor_kappa) - -end function qmckl_compute_nn_distance_rescaled_f - #+end_src - - #+begin_src c :tangle (eval h_private_func) :comments org :exports none -qmckl_exit_code qmckl_compute_nn_distance_rescaled ( - const qmckl_context context, - const int64_t nucl_num, - const double rescale_factor_kappa, - const double* coord, - double* const nn_distance_rescaled ); - #+end_src - - - #+CALL: generate_c_interface(table=qmckl_nn_distance_rescaled_args,rettyp="qmckl_exit_code",fname="qmckl_compute_nn_distance") - - #+RESULTS: - #+begin_src f90 :tangle (eval f) :comments org :exports none - integer(c_int32_t) function qmckl_compute_nn_distance_rescaled & - (context, nucl_num, rescale_factor_kappa, coord, nn_distance_rescaled) & - bind(C) result(info) - - use, intrinsic :: iso_c_binding - implicit none - - integer (c_int64_t) , intent(in) , value :: context - integer (c_int64_t) , intent(in) , value :: nucl_num - real (c_double ) , intent(in) , value :: rescale_factor_kappa - real (c_double ) , intent(in) :: coord(nucl_num,3) - real (c_double ) , intent(out) :: nn_distance_rescaled(nucl_num,nucl_num) - - integer(c_int32_t), external :: qmckl_compute_nn_distance_rescaled_f - info = qmckl_compute_nn_distance_rescaled_f & - (context, nucl_num, rescale_factor_kappa, coord, nn_distance_rescaled) - - end function qmckl_compute_nn_distance_rescaled - #+end_src - -*** Test - - #+begin_src c :tangle (eval c_test) -/* Reference input data */ -/* TODO */ - -//assert(qmckl_nucleus_provided(context)); -// -//double distance[nucl_num*nucl_num]; -//rc = qmckl_get_nucleus_nn_distance(context, distance, nucl_num*nucl_num); -//assert(distance[0] == 0.); -//assert(distance[1] == distance[nucl_num]); -//assert(fabs(distance[1]-2.070304721365169) < 1.e-12); - - #+end_src - ** Nuclear repulsion energy \[ diff --git a/org/qmckl_trexio.org b/org/qmckl_trexio.org index a94bd85..4490957 100644 --- a/org/qmckl_trexio.org +++ b/org/qmckl_trexio.org @@ -1097,37 +1097,25 @@ qmckl_trexio_read(const qmckl_context context, const char* file_name, const int6 rc = qmckl_trexio_read_electron_X(context, file); if (rc != QMCKL_SUCCESS) { trexio_close(file); - return qmckl_failwith( context, - QMCKL_FAILURE, - "qmckl_trexio_read", - "Error reading electron"); + return rc; } rc = qmckl_trexio_read_nucleus_X(context, file); if (rc != QMCKL_SUCCESS) { trexio_close(file); - return qmckl_failwith( context, - QMCKL_FAILURE, - "qmckl_trexio_read", - "Error reading nucleus"); + return rc; } rc = qmckl_trexio_read_ao_X(context, file); if (rc != QMCKL_SUCCESS) { trexio_close(file); - return qmckl_failwith( context, - QMCKL_FAILURE, - "qmckl_trexio_read", - "Error reading AOs"); + return rc; } rc = qmckl_trexio_read_mo_X(context, file); if (rc != QMCKL_SUCCESS) { trexio_close(file); - return qmckl_failwith( context, - QMCKL_FAILURE, - "qmckl_trexio_read", - "Error reading MOs"); + return rc; } trexio_close(file); @@ -1149,27 +1137,19 @@ qmckl_trexio_read(const qmckl_context context, const char* file_name, const int6 #ifdef HAVE_TREXIO qmckl_exit_code rc; -char fname[256]; -char message[256]; +char filename[256]; #ifndef QMCKL_TEST_DIR #error "QMCKL_TEST_DIR is not defined" #endif -strncpy(fname, QMCKL_TEST_DIR,255); -strncat(fname, "/chbrclf", 255); -printf("Test file: %s\n", fname); +strncpy(filename, QMCKL_TEST_DIR,255); +strncat(filename, "/chbrclf", 255); +printf("Test file: %s\n", filename); -rc = qmckl_trexio_read(context, fname, 255); +rc = qmckl_trexio_read(context, filename, 255); -if (rc != QMCKL_SUCCESS) { - printf("%s\n", qmckl_string_of_error(rc)); - qmckl_get_error(context, &rc, fname, message); - printf("%s\n", fname); - printf("%s\n", message); - } - -assert ( rc == QMCKL_SUCCESS ); +qmckl_check(context, rc); #+end_src @@ -1179,11 +1159,11 @@ assert ( rc == QMCKL_SUCCESS ); printf("Electrons\n"); int64_t up_num, dn_num; rc = qmckl_get_electron_up_num(context, &up_num); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert (up_num == chbrclf_elec_up_num); rc = qmckl_get_electron_down_num(context, &dn_num); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert (dn_num == chbrclf_elec_dn_num); #+end_src @@ -1195,13 +1175,13 @@ printf("Nuclei\n"); int64_t nucl_num; rc = qmckl_get_nucleus_num(context, &nucl_num); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert (nucl_num == chbrclf_nucl_num); printf("Nuclear charges\n"); double * charge = (double*) malloc (nucl_num * sizeof(double)); rc = qmckl_get_nucleus_charge(context, charge, nucl_num); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); for (int i=0 ; i