From 3f33db6887a8bc834c0f6e3186aa86b2e873af2e Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 13 Mar 2023 15:32:35 +0100 Subject: [PATCH] Transposed en_distance --- org/qmckl_ao.org | 2 +- org/qmckl_electron.org | 43 ++++++------- org/qmckl_jastrow.org | 134 ++++++++++++++++++++--------------------- org/qmckl_point.org | 6 +- 4 files changed, 90 insertions(+), 95 deletions(-) diff --git a/org/qmckl_ao.org b/org/qmckl_ao.org index 6bc84a6..1e9bb0a 100644 --- a/org/qmckl_ao.org +++ b/org/qmckl_ao.org @@ -6305,7 +6305,7 @@ assert( fabs(ao_value[26][224] - ( 7.175045873560788e-10)) < 1.e-14 ); :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: -*** Unoptimized version +*** Reference version #+NAME: qmckl_ao_vgl_args_doc | Variable | Type | In/Out | Description | |-----------------------+-----------------------------------+--------+----------------------------------------------| diff --git a/org/qmckl_electron.org b/org/qmckl_electron.org index af06ba3..24b7630 100644 --- a/org/qmckl_electron.org +++ b/org/qmckl_electron.org @@ -95,7 +95,7 @@ int main() { |-------------------------------------+----------------------------------------+----------------------------------------------------------------------| | ~ee_distance~ | ~double[walker.num][num][num]~ | Electron-electron distances | | ~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~ | ~double[walker.num][num][nucl_num]~ | Electron-nucleus distances | | ~en_distance_date~ | ~uint64_t~ | Last modification date of the electron-electron distances | | ~ee_potential~ | ~double[walker.num]~ | Electron-electron potential energy | | ~ee_potential_date~ | ~uint64_t~ | Last modification date of the electron-electron potential | @@ -1235,7 +1235,7 @@ qmckl_exit_code qmckl_provide_en_distance(qmckl_context context) | ~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~ | ~double[walk_num][nucl_num][elec_num]~ | out | Electron-nucleus distances | + | ~en_distance~ | ~double[walk_num][elec_num][nucl_num]~ | out | Electron-nucleus distances | #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_en_distance_f(context, elec_num, nucl_num, walk_num, elec_coord, nucl_coord, en_distance) & @@ -1248,7 +1248,7 @@ integer function qmckl_compute_en_distance_f(context, elec_num, nucl_num, walk_n 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(elec_num,nucl_num,walk_num) + double precision , intent(out) :: en_distance(nucl_num,elec_num,walk_num) integer*8 :: k @@ -1274,15 +1274,10 @@ integer function qmckl_compute_en_distance_f(context, elec_num, nucl_num, walk_n return endif - do k=1,walk_num - info = qmckl_distance(context, 'T', 'T', elec_num, nucl_num, & - elec_coord(1,k,1), elec_num * walk_num, & + info = qmckl_distance(context, 'T', 'T', nucl_num, elec_num * walk_num, & nucl_coord, nucl_num, & - en_distance(1,1,k), elec_num) - if (info /= QMCKL_SUCCESS) then - exit - endif - end do + elec_coord, elec_num * walk_num, & + en_distance, nucl_num) end function qmckl_compute_en_distance_f #+end_src @@ -1315,7 +1310,7 @@ qmckl_exit_code qmckl_compute_en_distance ( 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(elec_num,nucl_num,walk_num) + real (c_double ) , intent(out) :: en_distance(nucl_num,elec_num,walk_num) integer(c_int32_t), external :: qmckl_compute_en_distance_f info = qmckl_compute_en_distance_f & @@ -1368,7 +1363,7 @@ qmckl_check(context, rc); assert(qmckl_nucleus_provided(context)); -double en_distance[walk_num][nucl_num][elec_num]; +double en_distance[walk_num][elec_num][nucl_num]; rc = qmckl_get_electron_en_distance(context, &(en_distance[0][0][0])); qmckl_check(context, rc); @@ -1378,19 +1373,19 @@ qmckl_check(context, rc); assert(fabs(en_distance[0][0][0] - 7.546738741619978) < 1.e-12); // (1,2,1) -assert(fabs(en_distance[0][1][0] - 8.77102435246984) < 1.e-12); +assert(fabs(en_distance[0][0][1] - 8.77102435246984) < 1.e-12); // (2,1,1) -assert(fabs(en_distance[0][0][1] - 3.698922010513608) < 1.e-12); +assert(fabs(en_distance[0][1][0] - 3.698922010513608) < 1.e-12); // (1,1,2) assert(fabs(en_distance[1][0][0] - 5.824059436060509) < 1.e-12); // (1,2,2) -assert(fabs(en_distance[1][1][0] - 7.080482110317645) < 1.e-12); +assert(fabs(en_distance[1][0][1] - 7.080482110317645) < 1.e-12); // (2,1,2) -assert(fabs(en_distance[1][0][1] - 3.1804527583077356) < 1.e-12); +assert(fabs(en_distance[1][1][0] - 3.1804527583077356) < 1.e-12); #+end_src @@ -1512,7 +1507,7 @@ qmckl_exit_code qmckl_provide_en_potential(qmckl_context context) | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | | ~walk_num~ | ~int64_t~ | in | Number of walkers | | ~charge~ | ~double[nucl_num]~ | in | charge of nucleus | - | ~en_distance~ | ~double[walk_num][nucl_num][elec_num]~ | in | Electron-electron rescaled distances | + | ~en_distance~ | ~double[walk_num][elec_num][nucl_num]~ | in | Electron-electron distances | | ~en_potential~ | ~double[walk_num]~ | out | Electron-electron potential | #+begin_src f90 :comments org :tangle (eval f) :noweb yes @@ -1526,7 +1521,7 @@ integer function qmckl_compute_en_potential_f(context, elec_num, nucl_num, walk_ integer*8 , intent(in) :: nucl_num integer*8 , intent(in) :: walk_num double precision , intent(in) :: charge(nucl_num) - double precision , intent(in) :: en_distance(elec_num,nucl_num,walk_num) + double precision , intent(in) :: en_distance(nucl_num,elec_num,walk_num) double precision , intent(out) :: en_potential(walk_num) integer*8 :: nw, i, j @@ -1550,10 +1545,10 @@ integer function qmckl_compute_en_potential_f(context, elec_num, nucl_num, walk_ en_potential = 0.0d0 do nw=1,walk_num - do j=1,nucl_num - do i=1,elec_num - if (dabs(en_distance(i,j,nw)) > 1e-5) then - en_potential(nw) = en_potential(nw) - charge(j)/(en_distance(i,j,nw)) + do i=1,elec_num + do j=1,nucl_num + if (dabs(en_distance(j,i,nw)) > 1.d-6) then + en_potential(nw) = en_potential(nw) - charge(j)/(en_distance(j,i,nw)) endif end do end do @@ -1592,7 +1587,7 @@ end function qmckl_compute_en_potential_f integer (c_int64_t) , intent(in) , value :: nucl_num integer (c_int64_t) , intent(in) , value :: walk_num real (c_double ) , intent(in) :: charge(nucl_num) - real (c_double ) , intent(in) :: en_distance(elec_num,nucl_num,walk_num) + real (c_double ) , intent(in) :: en_distance(nucl_num,elec_num,walk_num) real (c_double ) , intent(out) :: en_potential(walk_num) integer(c_int32_t), external :: qmckl_compute_en_potential_f diff --git a/org/qmckl_jastrow.org b/org/qmckl_jastrow.org index dc631b8..8739f44 100644 --- a/org/qmckl_jastrow.org +++ b/org/qmckl_jastrow.org @@ -11,7 +11,7 @@ \[ J(\mathbf{r},\mathbf{R}) = J_{\text{eN}}(\mathbf{r},\mathbf{R}) + J_{\text{ee}}(\mathbf{r}) + J_{\text{eeN}}(\mathbf{r},\mathbf{R}) \] - + In the following, we use the notations $r_{ij} = |\mathbf{r}_i - \mathbf{r}_j|$ and $R_{i\alpha} = |\mathbf{r}_i - \mathbf{R}_\alpha|$. @@ -19,7 +19,7 @@ \[ J_{\text{eN}}(\mathbf{r},\mathbf{R}) = - \sum_{\alpha=1}^{N_\text{nucl}} \sum_{i=1}^{N_\text{elec}} + \sum_{\alpha=1}^{N_\text{nucl}} \sum_{i=1}^{N_\text{elec}} \frac{a_{1\,\alpha}\, f_\alpha(R_{i\,\alpha})}{1+a_{2\,\alpha}\, f_\alpha(R_{i\alpha})} + \sum_{p=2}^{N_\text{ord}^a} a_{p+1\,\alpha}\, [f_\alpha(R_{i\alpha})]^p - J_{eN}^{\infty \alpha} \] @@ -640,7 +640,7 @@ qmckl_set_jastrow_a_vector(qmckl_context context, "qmckl_set_jastrow_a_vector", "aord_num not initialized"); } - + int64_t type_nucl_num = ctx->jastrow.type_nucl_num; if (type_nucl_num <= 0) { @@ -825,7 +825,7 @@ qmckl_set_jastrow_rescale_factor_ee(qmckl_context context, const double rescale_factor_ee) { int32_t mask = 1 << 8; - + <> if (rescale_factor_ee <= 0.0) { @@ -979,7 +979,7 @@ interface integer(c_int64_t), intent(in), value :: size_max double precision, intent(in) :: kappa_en(size_max) end function qmckl_set_jastrow_rescale_factor_en - + integer(qmckl_exit_code) function qmckl_set_jastrow_aord_num (context, & aord_num) bind(C) use, intrinsic :: iso_c_binding @@ -988,7 +988,7 @@ interface integer (qmckl_context) , intent(in) , value :: context integer(c_int64_t), intent(in), value :: aord_num end function qmckl_set_jastrow_aord_num - + integer(qmckl_exit_code) function qmckl_set_jastrow_bord_num (context, & bord_num) bind(C) use, intrinsic :: iso_c_binding @@ -997,7 +997,7 @@ interface integer (qmckl_context) , intent(in) , value :: context integer(c_int64_t), intent(in), value :: bord_num end function qmckl_set_jastrow_bord_num - + integer(qmckl_exit_code) function qmckl_set_jastrow_cord_num (context, & cord_num) bind(C) use, intrinsic :: iso_c_binding @@ -1015,7 +1015,7 @@ interface integer (qmckl_context) , intent(in) , value :: context integer(c_int64_t), intent(in), value :: type_nucl_num end function qmckl_set_jastrow_type_nucl_num - + integer(qmckl_exit_code) function qmckl_set_jastrow_type_nucl_vector (context, & type_nucl_vector, size_max) bind(C) use, intrinsic :: iso_c_binding @@ -1025,7 +1025,7 @@ interface integer(c_int64_t), intent(in), value :: size_max integer(c_int64_t), intent(in) :: type_nucl_vector(size_max) end function qmckl_set_jastrow_type_nucl_vector - + integer(qmckl_exit_code) function qmckl_set_jastrow_a_vector(context, & a_vector, size_max) bind(C) use, intrinsic :: iso_c_binding @@ -1035,7 +1035,7 @@ interface integer(c_int64_t), intent(in), value :: size_max double precision, intent(in) :: a_vector(size_max) end function qmckl_set_jastrow_a_vector - + integer(qmckl_exit_code) function qmckl_set_jastrow_b_vector(context, & b_vector, size_max) bind(C) use, intrinsic :: iso_c_binding @@ -1045,7 +1045,7 @@ interface integer(c_int64_t), intent(in), value :: size_max double precision, intent(in) :: b_vector(size_max) end function qmckl_set_jastrow_b_vector - + integer(qmckl_exit_code) function qmckl_set_jastrow_c_vector(context, & c_vector, size_max) bind(C) use, intrinsic :: iso_c_binding @@ -1055,7 +1055,7 @@ interface integer(c_int64_t), intent(in), value :: size_max double precision, intent(in) :: c_vector(size_max) end function qmckl_set_jastrow_c_vector - + end interface #+end_src @@ -1074,7 +1074,7 @@ qmckl_exit_code qmckl_get_jastrow_rescale_factor_ee (const qmckl_context contex 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: @@ -1455,7 +1455,7 @@ interface integer(c_int64_t), intent(in), value :: size_max double precision, intent(out) :: kappa_en(size_max) end function qmckl_get_jastrow_rescale_factor_en - + integer(qmckl_exit_code) function qmckl_get_jastrow_aord_num (context, & aord_num) bind(C) use, intrinsic :: iso_c_binding @@ -1464,7 +1464,7 @@ interface integer (qmckl_context) , intent(in), value :: context integer(c_int64_t), intent(out) :: aord_num end function qmckl_get_jastrow_aord_num - + integer(qmckl_exit_code) function qmckl_get_jastrow_bord_num (context, & bord_num) bind(C) use, intrinsic :: iso_c_binding @@ -1473,7 +1473,7 @@ interface integer (qmckl_context) , intent(in), value :: context integer(c_int64_t), intent(out) :: bord_num end function qmckl_get_jastrow_bord_num - + integer(qmckl_exit_code) function qmckl_get_jastrow_cord_num (context, & cord_num) bind(C) use, intrinsic :: iso_c_binding @@ -1491,7 +1491,7 @@ interface integer (qmckl_context) , intent(in), value :: context integer(c_int64_t), intent(out) :: type_nucl_num end function qmckl_get_jastrow_type_nucl_num - + integer(qmckl_exit_code) function qmckl_get_jastrow_type_nucl_vector (context, & type_nucl_vector, size_max) bind(C) use, intrinsic :: iso_c_binding @@ -1501,7 +1501,7 @@ interface integer(c_int64_t), intent(in), value :: size_max integer(c_int64_t), intent(out) :: type_nucl_vector(size_max) end function qmckl_get_jastrow_type_nucl_vector - + integer(qmckl_exit_code) function qmckl_get_jastrow_a_vector(context, & a_vector, size_max) bind(C) use, intrinsic :: iso_c_binding @@ -1511,7 +1511,7 @@ interface integer(c_int64_t), intent(in), value :: size_max double precision, intent(out) :: a_vector(size_max) end function qmckl_get_jastrow_a_vector - + integer(qmckl_exit_code) function qmckl_get_jastrow_b_vector(context, & b_vector, size_max) bind(C) use, intrinsic :: iso_c_binding @@ -1521,7 +1521,7 @@ interface integer(c_int64_t), intent(in), value :: size_max double precision, intent(out) :: b_vector(size_max) end function qmckl_get_jastrow_b_vector - + integer(qmckl_exit_code) function qmckl_get_jastrow_c_vector(context, & c_vector, size_max) bind(C) use, intrinsic :: iso_c_binding @@ -1531,7 +1531,7 @@ interface integer(c_int64_t), intent(in), value :: size_max double precision, intent(out) :: c_vector(size_max) end function qmckl_get_jastrow_c_vector - + end interface #+end_src @@ -1874,7 +1874,7 @@ qmckl_exit_code qmckl_compute_jastrow_asymp_jasb_doc (const qmckl_context contex const double rescale_factor_ee, double* const asymp_jasb ); #+end_src - + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_jastrow_asymp_jasb_hpc (const qmckl_context context, @@ -2067,7 +2067,7 @@ assert(fabs(asymp_jasb[1]-0.31567342786262853) < 1.e-12); \[ f_\text{ee} = \sum_{i,jjastrow.factor_ee = NULL; } } - + /* Allocate array */ if (ctx->jastrow.factor_ee == NULL) { @@ -2264,24 +2264,24 @@ integer function qmckl_compute_factor_ee_doc_f(context, walk_num, elec_num, up_n info = QMCKL_INVALID_CONTEXT return endif - + if (walk_num <= 0) then info = QMCKL_INVALID_ARG_2 return endif - + if (elec_num <= 0) then info = QMCKL_INVALID_ARG_3 return endif - + if (bord_num < 0) then info = QMCKL_INVALID_ARG_4 return endif - + factor_ee = 0.0d0 - + do nw =1, walk_num do j = 1, elec_num do i = 1, j - 1 @@ -2289,23 +2289,23 @@ integer function qmckl_compute_factor_ee_doc_f(context, walk_num, elec_num, up_n power_ser = 0.0d0 spin_fact = 1.0d0 ipar = 1 - + do p = 2, bord_num x = x * ee_distance_rescaled(i,j,nw) power_ser = power_ser + b_vector(p + 1) * x end do - + if(j <= up_num .or. i > up_num) then spin_fact = 0.5d0 ipar = 2 endif - + factor_ee(nw) = factor_ee(nw) + spin_fact * b_vector(1) * & ee_distance_rescaled(i,j,nw) / & (1.0d0 + b_vector(2) * & ee_distance_rescaled(i,j,nw)) & - + power_ser - asymp_jasb(ipar) - + + power_ser - asymp_jasb(ipar) + end do end do end do @@ -2417,8 +2417,8 @@ qmckl_exit_code qmckl_compute_factor_ee_hpc ( ipar = 1; } - factor_ee[nw] += spin_fact * b_vector[0] * - x1 / (1.0 + b_vector[1] * x1) + factor_ee[nw] += spin_fact * b_vector[0] * + x1 / (1.0 + b_vector[1] * x1) - asymp_jasb[ipar] + power_ser; } @@ -2624,7 +2624,7 @@ qmckl_exit_code qmckl_provide_jastrow_factor_ee_deriv_e(qmckl_context context) ctx->jastrow.factor_ee_deriv_e = NULL; } } - + /* Allocate array */ if (ctx->jastrow.factor_ee_deriv_e == NULL) { @@ -2804,11 +2804,11 @@ qmckl_exit_code qmckl_compute_factor_ee_deriv_e_hpc( 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; @@ -2818,7 +2818,7 @@ qmckl_exit_code qmckl_compute_factor_ee_deriv_e_hpc( return QMCKL_INVALID_ARG_4; } - + for (int nw = 0; nw < walk_num; ++nw) { for (int ii = 0; ii < 4; ++ii) { for (int j = 0; j < elec_num; ++j) { @@ -2826,7 +2826,7 @@ qmckl_exit_code qmckl_compute_factor_ee_deriv_e_hpc( } } } - + third = 1.0 / 3.0; for (int nw = 0; nw < walk_num; ++nw) { @@ -2836,14 +2836,14 @@ qmckl_exit_code qmckl_compute_factor_ee_deriv_e_hpc( if (fabs(x) < 1.0e-18) continue; for (int ii = 0; ii < 3; ++ii){ pow_ser_g[ii] = 0.0; - } + } spin_fact = 1.0; den = 1.0 + b_vector[1] * x; invden = 1.0 / den; invden2 = invden * invden; invden3 = invden2 * invden; xinv = 1.0 / (x + 1.0e-18); - + dx[0] = ee_distance_rescaled_deriv_e[0 \ + j * 4 + i * 4 * elec_num \ + nw * 4 * elec_num * elec_num]; @@ -2891,7 +2891,7 @@ qmckl_exit_code qmckl_compute_factor_ee_deriv_e_hpc( } } } - + return QMCKL_SUCCESS; } #+end_src @@ -2956,7 +2956,7 @@ integer(c_int32_t) function qmckl_compute_factor_ee_deriv_e_doc & end function qmckl_compute_factor_ee_deriv_e_doc #+end_src - + #+begin_src c :tangle (eval h_private_func) :comments org qmckl_exit_code qmckl_compute_factor_ee_deriv_e_hpc ( const qmckl_context context, @@ -2998,9 +2998,9 @@ integer(c_int32_t) function qmckl_compute_factor_ee_deriv_e_doc & double* const factor_ee_deriv_e ) { #ifdef HAVE_HPC - return qmckl_compute_factor_ee_deriv_e_hpc(context, walk_num, elec_num, up_num, bord_num, b_vector, ee_distance_rescaled, ee_distance_rescaled_deriv_e, factor_ee_deriv_e ); + return qmckl_compute_factor_ee_deriv_e_hpc(context, walk_num, elec_num, up_num, bord_num, b_vector, ee_distance_rescaled, ee_distance_rescaled_deriv_e, factor_ee_deriv_e ); #else - return qmckl_compute_factor_ee_deriv_e_doc(context, walk_num, elec_num, up_num, bord_num, b_vector, ee_distance_rescaled, ee_distance_rescaled_deriv_e, factor_ee_deriv_e ); + return qmckl_compute_factor_ee_deriv_e_doc(context, walk_num, elec_num, up_num, bord_num, b_vector, ee_distance_rescaled, ee_distance_rescaled_deriv_e, factor_ee_deriv_e ); #endif } #+end_src @@ -3295,7 +3295,7 @@ integer function qmckl_compute_jastrow_asymp_jasa_f(context, aord_num, type_nucl integer*8 :: i, j, p double precision :: kappa_inv, x, asym_one - + info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then @@ -3311,15 +3311,15 @@ integer function qmckl_compute_jastrow_asymp_jasa_f(context, aord_num, type_nucl do i=1,type_nucl_num kappa_inv = 1.0d0 / rescale_factor_en(i) - + asymp_jasa(i) = a_vector(1,i) * kappa_inv / (1.0d0 + a_vector(2,i) * kappa_inv) - + x = kappa_inv do p = 2, aord_num x = x * kappa_inv asymp_jasa(i) = asymp_jasa(i) + a_vector(p+1, i) * x end do - + end do end function qmckl_compute_jastrow_asymp_jasa_f @@ -3371,7 +3371,7 @@ qmckl_exit_code qmckl_compute_jastrow_asymp_jasa ( for (int i = 0 ; i <= type_nucl_num; ++i) { const double kappa_inv = 1.0 / rescale_factor_en[i]; asymp_jasa[i] = a_vector[aord_num*i] * kappa_inv / (1.0 + a_vector[1 + aord_num*i] * kappa_inv); - + double x = kappa_inv; for (int p = 1; p < aord_num; ++p){ x *= kappa_inv; @@ -3394,7 +3394,7 @@ qmckl_exit_code qmckl_compute_jastrow_asymp_jasa ( const int64_t type_nucl_num, const double* a_vector, const double* rescale_factor_en, - double* const asymp_jasa ); + double* const asymp_jasa ); #+end_src *** Test @@ -3577,7 +3577,7 @@ qmckl_exit_code qmckl_provide_jastrow_factor_en(qmckl_context context) if (rc != QMCKL_SUCCESS) { return rc; } - + ctx->jastrow.factor_en_date = ctx->date; } @@ -3829,7 +3829,7 @@ qmckl_exit_code qmckl_compute_factor_en ( const double* a_vector, const double* en_distance_rescaled, const double* asymp_jasa, - double* const factor_en ); + double* const factor_en ); #+end_src *** Test @@ -6501,7 +6501,7 @@ integer function qmckl_compute_een_rescaled_n_f( & integer*8 , intent(in) :: type_nucl_vector(nucl_num) integer*8 , intent(in) :: cord_num double precision , intent(in) :: rescale_factor_en(type_nucl_num) - double precision , intent(in) :: en_distance(elec_num,nucl_num,walk_num) + double precision , intent(in) :: en_distance(nucl_num,elec_num,walk_num) double precision , intent(out) :: een_rescaled_n(elec_num,nucl_num,0:cord_num,walk_num) double precision :: x integer*8 :: i, a, k, l, nw @@ -6536,16 +6536,16 @@ integer function qmckl_compute_een_rescaled_n_f( & ! Prepare table of exponentiated distances raised to appropriate power een_rescaled_n = 0.0d0 do nw = 1, walk_num - + ! prepare the actual een table een_rescaled_n(:, :, 0, nw) = 1.0d0 - + do a = 1, nucl_num do i = 1, elec_num - een_rescaled_n(i, a, 1, nw) = dexp(-rescale_factor_en(type_nucl_vector(a)) * en_distance(i, a, nw)) + een_rescaled_n(i, a, 1, nw) = dexp(-rescale_factor_en(type_nucl_vector(a)) * en_distance(a, i, nw)) end do end do - + do l = 2, cord_num do a = 1, nucl_num do i = 1, elec_num @@ -6604,7 +6604,7 @@ qmckl_exit_code qmckl_compute_een_rescaled_n ( for (int i = 0; i < elec_num; ++i) { een_rescaled_n[i + a*elec_num + nw * elec_num*nucl_num*(cord_num+1)] = 1.0; een_rescaled_n[i + a*elec_num + elec_num*nucl_num + nw*elec_num*nucl_num*(cord_num+1)] = - exp(-rescale_factor_en[type_nucl_vector[a]] * en_distance[i + a*elec_num + nw*elec_num*nucl_num]); + exp(-rescale_factor_en[type_nucl_vector[a]] * en_distance[a + i*nucl_num + nw*elec_num*nucl_num]); } } @@ -6612,7 +6612,7 @@ qmckl_exit_code qmckl_compute_een_rescaled_n ( 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 + (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)]; } } @@ -6671,7 +6671,7 @@ qmckl_exit_code qmckl_compute_een_rescaled_n ( end function qmckl_compute_een_rescaled_n #+end_src - + # #+CALL: generate_c_header(table=qmckl_factor_een_rescaled_n_args,rettyp=get_value("CRetType"),fname=get_value("Name")) #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none @@ -6917,7 +6917,7 @@ integer function qmckl_compute_factor_een_rescaled_n_deriv_e_f( & double precision , intent(in) :: rescale_factor_en(type_nucl_num) double precision , intent(in) :: coord_ee(elec_num,3,walk_num) double precision , intent(in) :: coord_en(nucl_num,3) - double precision , intent(in) :: en_distance(elec_num,nucl_num,walk_num) + double precision , intent(in) :: en_distance(nucl_num,elec_num,walk_num) double precision , intent(in) :: een_rescaled_n(elec_num,nucl_num,0:cord_num,walk_num) double precision , intent(out) :: een_rescaled_n_deriv_e(elec_num,4,nucl_num,0:cord_num,walk_num) double precision,dimension(:,:,:),allocatable :: elnuc_dist_deriv_e @@ -6960,7 +6960,7 @@ integer function qmckl_compute_factor_een_rescaled_n_deriv_e_f( & ! prepare the actual een table do a = 1, nucl_num do i = 1, elec_num - ria_inv = 1.0d0 / en_distance(i, a, nw) + ria_inv = 1.0d0 / en_distance(a, i, nw) do ii = 1, 3 elnuc_dist_deriv_e(ii, i, a) = (coord_ee(i, ii, nw) - coord_en(a, ii)) * ria_inv end do @@ -7441,7 +7441,7 @@ qmckl_exit_code qmckl_provide_dtmp_c(qmckl_context context) if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } - + qmckl_exit_code rc; qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); @@ -7487,7 +7487,7 @@ qmckl_exit_code qmckl_provide_dtmp_c(qmckl_context context) ctx->jastrow.een_rescaled_e_deriv_e, ctx->jastrow.een_rescaled_n, ctx->jastrow.dtmp_c); - + if (rc != QMCKL_SUCCESS) { return rc; } diff --git a/org/qmckl_point.org b/org/qmckl_point.org index fa30f14..a586756 100644 --- a/org/qmckl_point.org +++ b/org/qmckl_point.org @@ -19,7 +19,7 @@ coordinates of all the walkers. #include #include "qmckl_blas_private_type.h" #+end_src - + #+begin_src c :tangle (eval h_private_func) #ifndef QMCKL_POINT_HPF #define QMCKL_POINT_HPF @@ -318,8 +318,8 @@ qmckl_set_point (qmckl_context context, if (num != ctx->point.num) { if (ctx->point.coord.data != NULL) { - rc = qmckl_matrix_free(context, &(ctx->point.coord)); - assert (rc == QMCKL_SUCCESS); + rc = qmckl_matrix_free(context, &(ctx->point.coord)); + assert (rc == QMCKL_SUCCESS); } ctx->point.coord = qmckl_matrix_alloc(context, num, 3);