From 2427d1b56ebc836c603a8466ac6098fd2d8820f9 Mon Sep 17 00:00:00 2001 From: Gianfranco Abrusci Date: Tue, 22 Feb 2022 14:52:21 +0100 Subject: [PATCH] working qmckl_compute_een_rescaled_n --- org/qmckl_jastrow.org | 106 +++++++++++++++++++++++++----------------- 1 file changed, 64 insertions(+), 42 deletions(-) diff --git a/org/qmckl_jastrow.org b/org/qmckl_jastrow.org index 7b1f0b8..9867010 100644 --- a/org/qmckl_jastrow.org +++ b/org/qmckl_jastrow.org @@ -3855,6 +3855,70 @@ integer function qmckl_compute_een_rescaled_n_f(context, walk_num, elec_num, nuc end function qmckl_compute_een_rescaled_n_f #+end_src + #+begin_src c :comments org :tangle (eval c) :noweb yes +qmckl_exit_code qmckl_compute_een_rescaled_n ( + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t nucl_num, + const int64_t cord_num, + const double rescale_factor_kappa_en, + const double* en_distance, + double* const een_rescaled_n ) { + + + if (context == QMCKL_NULL_CONTEXT) { + return QMCKL_INVALID_CONTEXT; + } + + if (walk_num <= 0) { + return QMCKL_INVALID_ARG_2; + } + + if (elec_num <= 0) { + return QMCKL_INVALID_ARG_3; + } + + if (nucl_num <= 0) { + return QMCKL_INVALID_ARG_4; + } + + if (cord_num <= 0) { + return QMCKL_INVALID_ARG_5; + } + + // Prepare table of exponentiated distances raised to appropriate power + for (int i = 0; i < (walk_num*(cord_num+1)*nucl_num*elec_num); ++i) { + een_rescaled_n[i] = 17.0; + } + + for (int nw = 0; nw < walk_num; ++nw) { + for (int a = 0; a < nucl_num; ++a) { + for (int i = 0; i < elec_num; ++i) { + // 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 * \ + en_distance[i + a*elec_num + nw*elec_num*nucl_num]); + } + } + + for (int l = 2; l < (cord_num+1); ++l){ + for (int a = 0; a < nucl_num; ++a) { + for (int i = 0; i < elec_num; ++i) { + een_rescaled_n[i + a*elec_num + l*elec_num*nucl_num + nw*elec_num*nucl_num*(cord_num+1)] = een_rescaled_n[i + a*elec_num + (l-1)*elec_num*nucl_num + nw*elec_num*nucl_num*(cord_num+1)] *\ + een_rescaled_n[i + a*elec_num + elec_num*nucl_num + nw*elec_num*nucl_num*(cord_num+1)]; + } + } + } + + } + + return QMCKL_SUCCESS; +} + #+end_src + #+CALL: generate_c_header(table=qmckl_factor_een_rescaled_n_args,rettyp=get_value("CRetType"),fname=get_value("Name")) #+RESULTS: @@ -3870,47 +3934,6 @@ end function qmckl_compute_een_rescaled_n_f double* const een_rescaled_n ); #+end_src - #+CALL: generate_c_interface(table=qmckl_factor_een_rescaled_n_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_een_rescaled_n & - (context, & - walk_num, & - elec_num, & - nucl_num, & - cord_num, & - rescale_factor_kappa_en, & - en_distance, & - een_rescaled_n) & - 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 :: walk_num - integer (c_int64_t) , intent(in) , value :: elec_num - integer (c_int64_t) , intent(in) , value :: nucl_num - integer (c_int64_t) , intent(in) , value :: cord_num - real (c_double ) , intent(in) , value :: rescale_factor_kappa_en - real (c_double ) , intent(in) :: en_distance(nucl_num,elec_num,walk_num) - real (c_double ) , intent(out) :: een_rescaled_n(nucl_num,elec_num,0:cord_num,walk_num) - - integer(c_int32_t), external :: qmckl_compute_een_rescaled_n_f - info = qmckl_compute_een_rescaled_n_f & - (context, & - walk_num, & - elec_num, & - nucl_num, & - cord_num, & - rescale_factor_kappa_en, & - en_distance, & - een_rescaled_n) - - end function qmckl_compute_een_rescaled_n - #+end_src - *** Test #+begin_src python :results output :exports none :noweb yes @@ -3969,7 +3992,6 @@ assert(fabs(een_rescaled_n[0][1][0][4]-0.023391817607642338) < 1.e-12); assert(fabs(een_rescaled_n[0][2][1][3]-0.880957224822116) < 1.e-12); assert(fabs(een_rescaled_n[0][2][1][4]-0.027185942659395074) < 1.e-12); assert(fabs(een_rescaled_n[0][2][1][5]-0.01343938025140174) < 1.e-12); - #+end_src ** Electron-nucleus rescaled distances for each order and derivatives