From 4133a6ba0e9aad68ebf7fb3080c51d84bf27d83c Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 12 Sep 2023 16:41:26 +0200 Subject: [PATCH] Fixed bug in 3-body Jastrow. Tested in CHAMP->correct --- org/qmckl_jastrow_champ.org | 238 +++++++++++++++++++++++++++--------- 1 file changed, 177 insertions(+), 61 deletions(-) diff --git a/org/qmckl_jastrow_champ.org b/org/qmckl_jastrow_champ.org index 1aa36cea..4c833f0b 100644 --- a/org/qmckl_jastrow_champ.org +++ b/org/qmckl_jastrow_champ.org @@ -42,7 +42,7 @@ \sum_{p=2}^{N_{\text{ord}}} \sum_{k=0}^{p-1} \sum_{l=0}^{p-k-2\delta_{k,0}} - c_{lkp\alpha} \left[ g_\text{ee}({r}_{ij}) \right]^k + c_{lkp\alpha} \left[ g_\text{e}({r}_{ij}) \right]^k \left[ \left[ g_\alpha({R}_{i\alpha}) \right]^l + \left[ g_\alpha({R}_{j\alpha}) \right]^l \right] \left[ g_\alpha({R}_{i\,\alpha}) \, g_\alpha({R}_{j\alpha}) \right]^{(p-k-l)/2} \] @@ -60,7 +60,7 @@ $J_{\text{ee}}$ and $J_{\text{eN}}$ have an asymptotic value of zero. The eN and eeN parameters are the same of all identical nuclei. - The types of nuclei use zero-based indexing. + Warning: The types of nuclei use zero-based indexing. * Headers :noexport: #+begin_src elisp :noexport :results none @@ -581,12 +581,12 @@ qmckl_set_jastrow_champ_type_nucl_vector(qmckl_context context, for (int i=0 ; i= type_nucl_num) { return qmckl_failwith( context, QMCKL_INVALID_ARG_2, - "qmckl_set_type_nucl_vector", + "qmckl_set_jastrow_champ_type_nucl_vector", "Inconsistent values of type_nucl_vector (>=nucl_num). Values should use 0-based indexing as in C." ); } } @@ -595,7 +595,7 @@ qmckl_set_jastrow_champ_type_nucl_vector(qmckl_context context, qmckl_exit_code rc = qmckl_free(context, ctx->jastrow_champ.type_nucl_vector); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, - "qmckl_set_type_nucl_vector", + "qmckl_set_jastrow_champ_type_nucl_vector", "Unable to free ctx->jastrow_champ.type_nucl_vector"); } } @@ -676,7 +676,7 @@ qmckl_set_jastrow_champ_a_vector(qmckl_context context, if(new_array == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, - "qmckl_set_jastrow_champ_coefficient", + "qmckl_set_jastrow_champ_a_vector", NULL); } @@ -736,7 +736,7 @@ qmckl_set_jastrow_champ_b_vector(qmckl_context context, if(new_array == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, - "qmckl_set_jastrow_champ_coefficient", + "qmckl_set_jastrow_champ_b_vector", NULL); } @@ -807,7 +807,7 @@ qmckl_set_jastrow_champ_c_vector(qmckl_context context, if(new_array == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, - "qmckl_set_jastrow_champ_coefficient", + "qmckl_set_jastrow_champ_c_vector", NULL); } @@ -2008,7 +2008,7 @@ rc = qmckl_check(context, ); assert(rc == QMCKL_SUCCESS); rc = qmckl_check(context, - qmckl_set_jastrow_champ_c_vector(context, c_vector,dim_c_vector*type_nucl_num) + qmckl_set_jastrow_champ_c_vector(context, c_vector, dim_c_vector*type_nucl_num) ); assert(rc == QMCKL_SUCCESS); @@ -7589,15 +7589,15 @@ qmckl_exit_code qmckl_provide_dtmp_c(qmckl_context context) :END: #+NAME: qmckl_factor_c_vector_full_args - | Variable | Type | In/Out | Description | - |--------------------+----------------------------------------+--------+------------------------------| - | ~context~ | ~qmckl_context~ | in | Global state | - | ~nucl_num~ | ~int64_t~ | in | Number of atoms | - | ~dim_c_vector~ | ~int64_t~ | in | dimension of cord full table | - | ~type_nucl_num~ | ~int64_t~ | in | dimension of cord full table | - | ~type_nucl_vector~ | ~int64_t[nucl_num]~ | in | dimension of cord full table | + | Variable | Type | In/Out | Description | + |--------------------+---------------------------------------+--------+------------------------------| + | ~context~ | ~qmckl_context~ | in | Global state | + | ~nucl_num~ | ~int64_t~ | in | Number of atoms | + | ~dim_c_vector~ | ~int64_t~ | in | dimension of cord full table | + | ~type_nucl_num~ | ~int64_t~ | in | dimension of cord full table | + | ~type_nucl_vector~ | ~int64_t[nucl_num]~ | in | dimension of cord full table | | ~c_vector~ | ~double[dim_c_vector][type_nucl_num]~ | in | dimension of cord full table | - | ~c_vector_full~ | ~double[dim_c_vector][nucl_num]~ | out | Full list of coefficients | + | ~c_vector_full~ | ~double[nucl_num][dim_c_vector]~ | out | Full list of coefficients | #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_c_vector_full_doc_f( & @@ -7611,8 +7611,8 @@ integer function qmckl_compute_c_vector_full_doc_f( & integer*8 , intent(in) :: dim_c_vector integer*8 , intent(in) :: type_nucl_num integer*8 , intent(in) :: type_nucl_vector(nucl_num) - double precision , intent(in) :: c_vector(type_nucl_num, dim_c_vector) - double precision , intent(out) :: c_vector_full(nucl_num,dim_c_vector) + double precision , intent(in) :: c_vector(dim_c_vector, type_nucl_num) + double precision , intent(out) :: c_vector_full(nucl_num, dim_c_vector) double precision :: x integer*8 :: i, a, k, l, nw @@ -7625,7 +7625,7 @@ integer function qmckl_compute_c_vector_full_doc_f( & if (info /= QMCKL_SUCCESS) return do a = 1, nucl_num - c_vector_full(a,1:dim_c_vector) = c_vector(type_nucl_vector(a)+1,1:dim_c_vector) + c_vector_full(a,1:dim_c_vector) = c_vector(1:dim_c_vector, type_nucl_vector(a)+1) end do end function qmckl_compute_c_vector_full_doc_f @@ -7677,7 +7677,7 @@ qmckl_exit_code qmckl_compute_c_vector_full_hpc ( for (int i=0; i < dim_c_vector; ++i) { for (int a=0; a < nucl_num; ++a){ - c_vector_full[a + i*nucl_num] = c_vector[(type_nucl_vector[a])+i*type_nucl_num]; + c_vector_full[a + i*nucl_num] = c_vector[i + type_nucl_vector[a]*dim_c_vector]; } } @@ -7746,7 +7746,7 @@ qmckl_exit_code qmckl_compute_c_vector_full ( :FRetType: qmckl_exit_code :END: - #+NAME: qmckl_factor_lkpm_combined_index_args + #+NAME: lkpm_combined_index_args | Variable | Type | In/Out | Description | |-----------------------+----------------------------+--------+-------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | @@ -7755,7 +7755,7 @@ qmckl_exit_code qmckl_compute_c_vector_full ( | ~lkpm_combined_index~ | ~int64_t[4][dim_c_vector]~ | out | Full list of combined indices | #+begin_src f90 :comments org :tangle (eval f) :noweb yes -integer function qmckl_compute_lkpm_combined_index_f( & +integer function qmckl_compute_lkpm_combined_index_doc_f( & context, cord_num, dim_c_vector, lkpm_combined_index) & result(info) use qmckl @@ -7777,13 +7777,13 @@ integer function qmckl_compute_lkpm_combined_index_f( & kk = 0 do p = 2, cord_num do k = p - 1, 0, -1 - if (k .ne. 0) then + if (k /= 0) then lmax = p - k else lmax = p - k - 2 end if do l = lmax, 0, -1 - if (iand(p - k - l, 1_8) .eq. 1) cycle + if (iand(p - k - l, 1_8) .eq. 1_8) cycle m = (p - k - l)/2 kk = kk + 1 lkpm_combined_index(kk, 1) = l @@ -7794,11 +7794,11 @@ integer function qmckl_compute_lkpm_combined_index_f( & end do end do -end function qmckl_compute_lkpm_combined_index_f +end function qmckl_compute_lkpm_combined_index_doc_f #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes -qmckl_exit_code qmckl_compute_lkpm_combined_index ( +qmckl_exit_code qmckl_compute_lkpm_combined_index_hpc ( const qmckl_context context, const int64_t cord_num, const int64_t dim_c_vector, @@ -7834,8 +7834,85 @@ qmckl_exit_code qmckl_compute_lkpm_combined_index ( } #+end_src - # #+CALL: generate_c_header(table=qmckl_factor_lkpm_combined_index_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + #+CALL: generate_c_interface(table=lkpm_combined_index_args,rettyp=get_value("CRetType"),fname="qmckl_compute_lkpm_combined_index_doc") + + #+RESULTS: + #+begin_src f90 :tangle (eval f) :comments org :exports none + integer(c_int32_t) function qmckl_compute_lkpm_combined_index_doc & + (context, cord_num, dim_c_vector, lkpm_combined_index) & + 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 :: cord_num + integer (c_int64_t) , intent(in) , value :: dim_c_vector + integer (c_int64_t) , intent(out) :: lkpm_combined_index(dim_c_vector,4) + + integer(c_int32_t), external :: qmckl_compute_lkpm_combined_index_doc_f + info = qmckl_compute_lkpm_combined_index_doc_f & + (context, cord_num, dim_c_vector, lkpm_combined_index) + + end function qmckl_compute_lkpm_combined_index_doc + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes +qmckl_exit_code qmckl_compute_lkpm_combined_index ( + const qmckl_context context, + const int64_t cord_num, + const int64_t dim_c_vector, + int64_t* const lkpm_combined_index ) { + + #ifdef HAVE_HPC + return qmckl_compute_lkpm_combined_index_hpc(context, cord_num, dim_c_vector, lkpm_combined_index); + #else + return qmckl_compute_lkpm_combined_index_doc(context, cord_num, dim_c_vector, lkpm_combined_index); + #endif +} + #+end_src + + #+CALL: generate_c_header(table=lkpm_combined_index_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none + qmckl_exit_code qmckl_compute_lkpm_combined_index ( + const qmckl_context context, + const int64_t cord_num, + const int64_t dim_c_vector, + int64_t* const lkpm_combined_index ); + #+end_src + #+CALL: generate_c_header(table=lkpm_combined_index_args,rettyp=get_value("CRetType"),fname="qmckl_compute_lkpm_combined_index_doc") + + #+RESULTS: + #+begin_src c :tangle (eval h_func) :comments org + qmckl_exit_code qmckl_compute_lkpm_combined_index_doc ( + const qmckl_context context, + const int64_t cord_num, + const int64_t dim_c_vector, + int64_t* const lkpm_combined_index ); + #+end_src + + #+CALL: generate_c_header(table=lkpm_combined_index_args,rettyp=get_value("CRetType"),fname="qmckl_compute_lkpm_combined_index_hpc") + + #+RESULTS: + #+begin_src c :tangle (eval h_func) :comments org + qmckl_exit_code qmckl_compute_lkpm_combined_index_hpc ( + const qmckl_context context, + const int64_t cord_num, + const int64_t dim_c_vector, + int64_t* const lkpm_combined_index ); + #+end_src + + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none + qmckl_exit_code qmckl_compute_lkpm_combined_index ( + const qmckl_context context, + const int64_t cord_num, + const int64_t dim_c_vector, + int64_t* const lkpm_combined_index ); + #+end_src + #+CALL: generate_c_header(table=lkpm_combined_index_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_lkpm_combined_index ( const qmckl_context context, @@ -8526,6 +8603,20 @@ qmckl_exit_code qmckl_provide_jastrow_champ_factor_een(qmckl_context context) ctx->jastrow_champ.factor_een = factor_een; } + /* + rc = qmckl_compute_jastrow_champ_factor_een_naive(context, + ctx->electron.walker.num, + ctx->electron.num, + ctx->nucleus.num, + ctx->jastrow_champ.cord_num, + ctx->jastrow_champ.dim_c_vector, + ctx->jastrow_champ.c_vector_full, + ctx->jastrow_champ.lkpm_combined_index, + ctx->jastrow_champ.een_rescaled_e, + ctx->jastrow_champ.een_rescaled_n, + ctx->jastrow_champ.factor_een); + */ + rc = qmckl_compute_jastrow_champ_factor_een(context, ctx->electron.walker.num, ctx->electron.num, @@ -8566,8 +8657,8 @@ qmckl_exit_code qmckl_provide_jastrow_champ_factor_een(qmckl_context context) | ~dim_c_vector~ | ~int64_t~ | in | dimension of full coefficient vector | | ~c_vector_full~ | ~double[dim_c_vector][nucl_num]~ | in | full coefficient vector | | ~lkpm_combined_index~ | ~int64_t[4][dim_c_vector]~ | in | combined indices | - | ~een_rescaled_e~ | ~double[walk_num][elec_num][elec_num][0:cord_num]~ | in | Electron-nucleus rescaled | - | ~een_rescaled_n~ | ~double[walk_num][elec_num][nucl_num][0:cord_num]~ | in | Electron-nucleus rescaled factor | + | ~een_rescaled_e~ | ~double[walk_num][0:cord_num][elec_num][elec_num]~ | in | Electron-nucleus rescaled | + | ~een_rescaled_n~ | ~double[walk_num][0:cord_num][nucl_num][elec_num]~ | in | Electron-nucleus rescaled factor | | ~factor_een~ | ~double[walk_num]~ | out | Electron-nucleus jastrow | #+begin_src f90 :comments org :tangle (eval f) :noweb yes @@ -8582,8 +8673,8 @@ integer function qmckl_compute_jastrow_champ_factor_een_naive_f( & integer*8 , intent(in) :: walk_num, elec_num, cord_num, nucl_num, dim_c_vector integer*8 , intent(in) :: lkpm_combined_index(dim_c_vector,4) double precision , intent(in) :: c_vector_full(nucl_num, dim_c_vector) - double precision , intent(in) :: een_rescaled_e(0:cord_num, elec_num, elec_num, walk_num) - double precision , intent(in) :: een_rescaled_n(0:cord_num, nucl_num, elec_num, walk_num) + double precision , intent(in) :: een_rescaled_e(elec_num, elec_num, 0:cord_num, walk_num) + double precision , intent(in) :: een_rescaled_n(elec_num, nucl_num, 0:cord_num, walk_num) double precision , intent(out) :: factor_een(walk_num) integer*8 :: i, a, j, l, k, p, m, n, nw @@ -8597,34 +8688,58 @@ integer function qmckl_compute_jastrow_champ_factor_een_naive_f( & if (nucl_num <= 0) info = QMCKL_INVALID_ARG_4 if (cord_num < 0) info = QMCKL_INVALID_ARG_5 if (info /= QMCKL_SUCCESS) return - - factor_een = 0.0d0 + + +! do nw =1, walk_num +! factor_een(nw) = 0.0d0 +! do n = 1, dim_c_vector +! l = lkpm_combined_index(n, 1) +! k = lkpm_combined_index(n, 2) +! p = lkpm_combined_index(n, 3) +! m = lkpm_combined_index(n, 4) + +! do a = 1, nucl_num +! accu2 = 0.0d0 +! cn = c_vector_full(a, n) +! do j = 1, elec_num +! accu = 0.0d0 +! do i = 1, elec_num + +! accu = accu + een_rescaled_e(i,j,k,nw) * & +! een_rescaled_n(i,a,m,nw) +! end do +! accu2 = accu2 + accu * een_rescaled_n(j,a,m + l,nw) +! end do +! factor_een(nw) = factor_een(nw) + accu2 * cn +! end do +! end do +! end do do nw =1, walk_num - do n = 1, dim_c_vector - l = lkpm_combined_index(n, 1) - k = lkpm_combined_index(n, 2) - p = lkpm_combined_index(n, 3) - m = lkpm_combined_index(n, 4) + factor_een(nw) = 0.d0 + do n = 1, dim_c_vector + l = lkpm_combined_index(n, 1) + k = lkpm_combined_index(n, 2) + p = lkpm_combined_index(n, 3) + m = lkpm_combined_index(n, 4) - do a = 1, nucl_num - accu2 = 0.0d0 - cn = c_vector_full(a, n) - do j = 1, elec_num - accu = 0.0d0 - do i = 1, elec_num - accu = accu + een_rescaled_e(k,i,j,nw) * & - een_rescaled_n(m,a,i,nw) - !if(nw .eq. 1) then - ! print *,l,k,p,m,j,i,een_rescaled_e(k,i,j,nw), een_rescaled_n(m,a,i,nw), accu - !endif + do a = 1, nucl_num + accu2 = 0.0d0 + cn = c_vector_full(a, n) + print *, a, l, k, p, cn + do j = 1, elec_num + accu = 0.0d0 + do i = 1, j-1 + + accu = accu + een_rescaled_e(i,j,k,nw) * & + (een_rescaled_n(i,a,l,nw) + een_rescaled_n(j,a,l,nw)) * & + (een_rescaled_n(i,a,m,nw) * een_rescaled_n(j,a,m,nw)) + end do + accu2 = accu2 + accu + end do + factor_een(nw) = factor_een(nw) + accu2 * cn end do - accu2 = accu2 + accu * een_rescaled_n(m + l,a,j,nw) - !print *, l,m,nw,accu, accu2, een_rescaled_n(m + l, a, j, nw), cn, factor_een(nw) - end do - factor_een(nw) = factor_een(nw) + accu2 * cn - end do - end do + end do end do end function qmckl_compute_jastrow_champ_factor_een_naive_f @@ -8786,11 +8901,12 @@ qmckl_compute_jastrow_champ_factor_een_doc (const qmckl_context context, const int64_t dim_c_vector, const double* c_vector_full, const int64_t* lkpm_combined_index, - const double* een_rescaled_e, + const double* tmp_c, const double* een_rescaled_n, double* const factor_een ); qmckl_exit_code + qmckl_compute_jastrow_champ_factor_een (const qmckl_context context, const int64_t walk_num, const int64_t elec_num, @@ -8799,7 +8915,7 @@ qmckl_compute_jastrow_champ_factor_een (const qmckl_context context, const int64_t dim_c_vector, const double* c_vector_full, const int64_t* lkpm_combined_index, - const double* een_rescaled_e, + const double* tmp_c, const double* een_rescaled_n, double* const factor_een ); #+end_src @@ -8814,17 +8930,17 @@ qmckl_compute_jastrow_champ_factor_een (const qmckl_context context, const int64_t dim_c_vector, const double* c_vector_full, const int64_t* lkpm_combined_index, - const double* een_rescaled_e, + const double* tmp_c, const double* een_rescaled_n, double* const factor_een ) { #ifdef HAVE_HPC return qmckl_compute_jastrow_champ_factor_een_doc (context, walk_num, elec_num, nucl_num, cord_num, dim_c_vector, - c_vector_full, lkpm_combined_index, een_rescaled_e, een_rescaled_n, + c_vector_full, lkpm_combined_index, tmp_c, een_rescaled_n, factor_een ); #else return qmckl_compute_jastrow_champ_factor_een_doc (context, walk_num, elec_num, nucl_num, cord_num, dim_c_vector, - c_vector_full, lkpm_combined_index, een_rescaled_e, een_rescaled_n, + c_vector_full, lkpm_combined_index, tmp_c, een_rescaled_n, factor_een ); #endif }