diff --git a/org/qmckl_jastrow.org b/org/qmckl_jastrow.org index 340fd36..a927c53 100644 --- a/org/qmckl_jastrow.org +++ b/org/qmckl_jastrow.org @@ -2835,7 +2835,7 @@ assert(fabs(factor_en_deriv_e[0][3][0]+0.9667363412285741 ) < 1.e-12); #+end_src - + ** Electron-electron rescaled distances for each order ~een_rescaled_e~ stores the table of the rescaled distances between all @@ -2846,7 +2846,7 @@ assert(fabs(factor_en_deriv_e[0][3][0]+0.9667363412285741 ) < 1.e-12); \] where \(C_{ij}\) is the matrix of electron-electron distances. - + *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes @@ -2993,10 +2993,12 @@ integer function qmckl_compute_een_rescaled_e_f(context, walk_num, elec_num, cor return endif + ! Prepare table of exponentiated distances raised to appropriate power + een_rescaled_e = 0.0d0 + do nw = 1, walk_num + een_rescaled_e_ij = 0.0d0 een_rescaled_e_ij(:, 1) = 1.0d0 - ! Prepare table of exponentiated distances raised to appropriate power - do nw = 1, walk_num k = 0 do j = 1, elec_num do i = 1, j - 1 @@ -3007,13 +3009,12 @@ integer function qmckl_compute_een_rescaled_e_f(context, walk_num, elec_num, cor do l = 2, cord_num do k = 1, elec_num * (elec_num - 1)/2 - een_rescaled_e_ij(k, l + 1) = een_rescaled_e_ij(k, l + 1 - 1) * een_rescaled_e_ij(k, 1) + een_rescaled_e_ij(k, l + 1) = een_rescaled_e_ij(k, l + 1 - 1) * een_rescaled_e_ij(k, 2) end do end do ! prepare the actual een table - een_rescaled_e = 0.0d0 - een_rescaled_e(0, :, :, :) = 1.0d0 + een_rescaled_e(0, :, :, nw) = 1.0d0 do l = 1, cord_num k = 0 do j = 1, elec_num @@ -3090,7 +3091,7 @@ een_rescaled_e_ij[:,0] = 1.0 k = 0 for j in range(elec_num): - for i in range(j - 1): + for i in range(j): een_rescaled_e_ij[k, 1] = np.exp(-kappa * elec_dist[i, j]) k = k + 1 @@ -3104,7 +3105,7 @@ een_rescaled_e[:,:,0] = 1.0 for l in range(1,cord_num+1): k = 0 for j in range(elec_num): - for i in range(j - 1): + for i in range(j): x = een_rescaled_e_ij[k, l] een_rescaled_e[i, j, l] = x een_rescaled_e[j, i, l] = x @@ -3134,18 +3135,314 @@ double een_rescaled_e[walk_num][elec_num][elec_num][(cord_num + 1)]; rc = qmckl_get_jastrow_een_rescaled_e(context, &(een_rescaled_e[0][0][0][0])); // value of (0,2,1) -//printf("%10.15f = \n", een_rescaled_e[0][0][2][1]); -//printf("%10.15f = \n", een_rescaled_e[0][0][3][1]); -//printf("%10.15f = \n", een_rescaled_e[0][0][4][1]); -//printf("%10.15f = \n", een_rescaled_e[0][1][3][2]); -//printf("%10.15f = \n", een_rescaled_e[0][1][4][2]); -//printf("%10.15f = \n", een_rescaled_e[0][1][5][2]); -//assert(fabs(een_rescaled_e[0][0][2][1]-) < 1.e-12); -//assert(fabs(een_rescaled_e[0][0][3][1]-) < 1.e-12); -//assert(fabs(een_rescaled_e[0][0][4][1]-) < 1.e-12); -//assert(fabs(een_rescaled_e[0][1][3][2]-) < 1.e-12); -//assert(fabs(een_rescaled_e[0][1][4][2]-) < 1.e-12); -//assert(fabs(een_rescaled_e[0][1][5][2]-) < 1.e-12); +assert(fabs(een_rescaled_e[0][0][2][1]-0.08084493981483197) < 1.e-12); +assert(fabs(een_rescaled_e[0][0][3][1]-0.1066745707571846) < 1.e-12); +assert(fabs(een_rescaled_e[0][0][4][1]-0.01754273169464735) < 1.e-12); +assert(fabs(een_rescaled_e[0][1][3][2]-0.02214680362033448) < 1.e-12); +assert(fabs(een_rescaled_e[0][1][4][2]-0.0005700154999202759) < 1.e-12); +assert(fabs(een_rescaled_e[0][1][5][2]-0.3424402276009091) < 1.e-12); + + #+end_src + +** Electron-nucleus rescaled distances for each order + + ~een_rescaled_n~ stores the table of the rescaled distances between + electrons and nucleii raised to the power \(p\) defined by ~cord_num~: + + \[ + C_{ia,p} = \left( 1 - \exp{-\kappa C_{ia}} \right)^p + \] + + where \(C_{ia}\) is the matrix of electron-nucleus distances. + +*** Get + + #+begin_src c :comments org :tangle (eval h_func) :noweb yes +qmckl_exit_code qmckl_get_jastrow_een_rescaled_n(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_een_rescaled_n(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_een_rescaled_n(context); + if (rc != QMCKL_SUCCESS) return rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + size_t sze = ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1); + memcpy(distance_rescaled, ctx->jastrow.een_rescaled_n, 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_een_rescaled_n(qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_provide_een_rescaled_n(qmckl_context context) +{ + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + /* Check if ee distance is provided */ + qmckl_exit_code rc = qmckl_provide_en_distance(context); + if(rc != QMCKL_SUCCESS) return rc; + + /* Compute if necessary */ + if (ctx->date > ctx->jastrow.een_rescaled_n_date) { + + /* Allocate array */ + if (ctx->jastrow.een_rescaled_n == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->electron.num * ctx->nucleus.num * + ctx->electron.walk_num * (ctx->jastrow.cord_num + 1) * sizeof(double); + double* een_rescaled_n = (double*) qmckl_malloc(context, mem_info); + + if (een_rescaled_n == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_een_rescaled_n", + NULL); + } + ctx->jastrow.een_rescaled_n = een_rescaled_n; + } + + qmckl_exit_code rc = + qmckl_compute_een_rescaled_n(context, + ctx->electron.walk_num, + ctx->electron.num, + ctx->nucleus.num, + ctx->jastrow.cord_num, + ctx->electron.rescale_factor_kappa_en, + ctx->electron.en_distance, + ctx->jastrow.een_rescaled_n); + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->jastrow.een_rescaled_n_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + #+end_src + +*** Compute + :PROPERTIES: + :Name: qmckl_compute_een_rescaled_n + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_factor_een_rescaled_n_args + | qmckl_context | context | in | Global state | + | int64_t | walk_num | in | Number of walkers | + | int64_t | elec_num | in | Number of electrons | + | int64_t | nucl_num | in | Number of atoms | + | int64_t | cord_num | in | Order of polynomials | + | double | rescale_factor_kappa_en | in | Factor to rescale ee distances | + | double | en_distance[walk_num][elec_num][nucl_num] | in | Electron-nucleus distances | + | double | een_rescaled_n[walk_num][elec_num][nucl_num][0:cord_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, & + en_distance, een_rescaled_n) & + result(info) + use qmckl + implicit none + integer(qmckl_context), intent(in) :: context + integer*8 , intent(in) :: walk_num + integer*8 , intent(in) :: elec_num + integer*8 , intent(in) :: nucl_num + integer*8 , intent(in) :: cord_num + double precision , intent(in) :: rescale_factor_kappa_en + double precision , intent(in) :: en_distance(elec_num,nucl_num,walk_num) + double precision , intent(out) :: een_rescaled_n(0:cord_num,nucl_num,elec_num,walk_num) + double precision :: x + integer*8 :: i, a, k, l, nw + + info = QMCKL_SUCCESS + + if (context == QMCKL_NULL_CONTEXT) then + info = QMCKL_INVALID_CONTEXT + return + endif + + if (walk_num <= 0) then + info = QMCKL_INVALID_ARG_2 + return + endif + + if (elec_num <= 0) then + info = QMCKL_INVALID_ARG_3 + return + endif + + if (nucl_num <= 0) then + info = QMCKL_INVALID_ARG_4 + return + endif + + if (cord_num <= 0) then + info = QMCKL_INVALID_ARG_5 + return + endif + + ! Prepare table of exponentiated distances raised to appropriate power + een_rescaled_n = 0.0d0 + do nw = 1, walk_num + + ! prepare the actual een table + een_rescaled_n(0, :, :, nw) = 1.0d0 + + do a = 1, nucl_num + do i = 1, elec_num + een_rescaled_n(1, a, i, nw) = dexp(-rescale_factor_kappa_en * en_distance(i, a, nw)) + end do + end do + + do l = 2, cord_num + do a = 1, nucl_num + do i = 1, elec_num + een_rescaled_n(l, a, i, nw) = een_rescaled_n(l - 1, a, i, nw) * een_rescaled_n(1, a, i, nw) + end do + end do + end do + end do + +end function qmckl_compute_een_rescaled_n_f + #+end_src + + #+CALL: generate_c_header(table=qmckl_factor_een_rescaled_n_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + + #+RESULTS: + #+begin_src c :tangle (eval h_func) :comments org + 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 ); + #+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(0:cord_num,nucl_num,elec_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 +import numpy as np + +<> + +elec_coord = np.array(elec_coord)[0] +nucl_coord = np.array(nucl_coord) +elnuc_dist = np.zeros(shape=(elec_num, nucl_num),dtype=float) +for i in range(elec_num): + for a in range(nucl_num): + elnuc_dist[i, a] = np.linalg.norm(elec_coord[i] - nucl_coord[:,a]) + +kappa = 1.0 + +een_rescaled_n = np.zeros(shape=(nucl_num, elec_num, cord_num + 1), dtype=float) +een_rescaled_n[:,:,0] = 1.0 + +for a in range(nucl_num): + for i in range(elec_num): + een_rescaled_n[a, i, 1] = np.exp(-kappa * elnuc_dist[i, a]) + +for l in range(2,cord_num+1): + for a in range(nucl_num): + for i in range(elec_num): + een_rescaled_n[a, i, l] = een_rescaled_n[a, i, l - 1] * een_rescaled_n[a, i, 1] + +print(" een_rescaled_n[0, 2, 1] = ",een_rescaled_n[0, 2, 1]) +print(" een_rescaled_n[0, 3, 1] = ",een_rescaled_n[0, 3, 1]) +print(" een_rescaled_n[0, 4, 1] = ",een_rescaled_n[0, 4, 1]) +print(" een_rescaled_n[1, 3, 2] = ",een_rescaled_n[1, 3, 2]) +print(" een_rescaled_n[1, 4, 2] = ",een_rescaled_n[1, 4, 2]) +print(" een_rescaled_n[1, 5, 2] = ",een_rescaled_n[1, 5, 2]) + #+end_src + + #+RESULTS: + : een_rescaled_n[0, 2, 1] = 0.10612983920006765 + : een_rescaled_n[0, 3, 1] = 0.135652809635553 + : een_rescaled_n[0, 4, 1] = 0.023391817607642338 + : een_rescaled_n[1, 3, 2] = 0.880957224822116 + : een_rescaled_n[1, 4, 2] = 0.027185942659395074 + : een_rescaled_n[1, 5, 2] = 0.01343938025140174 + + #+begin_src c :tangle (eval c_test) +assert(qmckl_electron_provided(context)); + +double een_rescaled_n[walk_num][elec_num][nucl_num][(cord_num + 1)]; +rc = qmckl_get_jastrow_een_rescaled_n(context, &(een_rescaled_n[0][0][0][0])); + +// value of (0,2,1) +assert(fabs(een_rescaled_n[0][2][0][1]-0.10612983920006765) < 1.e-12); +assert(fabs(een_rescaled_n[0][3][0][1]-0.135652809635553) < 1.e-12); +assert(fabs(een_rescaled_n[0][4][0][1]-0.023391817607642338) < 1.e-12); +assert(fabs(een_rescaled_n[0][3][1][2]-0.880957224822116) < 1.e-12); +assert(fabs(een_rescaled_n[0][4][1][2]-0.027185942659395074) < 1.e-12); +assert(fabs(een_rescaled_n[0][5][1][2]-0.01343938025140174) < 1.e-12); #+end_src