diff --git a/org/qmckl_electron.org b/org/qmckl_electron.org index 153f48b..03bb3e0 100644 --- a/org/qmckl_electron.org +++ b/org/qmckl_electron.org @@ -1009,6 +1009,241 @@ assert(fabs(ee_distance[elec_num*elec_num+1]-6.5517646321055665) < 1.e-12); #+end_src +** Electron-electron rescaled distances + +*** Get + + #+begin_src c :comments org :tangle (eval h_func) :noweb yes +qmckl_exit_code qmckl_get_electron_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_electron_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* const) context; + assert (ctx != NULL); + + size_t sze = ctx->electron.num * ctx->electron.num * ctx->electron.walk_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* const) context; + assert (ctx != NULL); + + + /* Compute if necessary */ + if (ctx->electron.coord_new_date > ctx->electron.ee_distance_rescaled_date) { + + /* 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.walk_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.kappa_en, + ctx->electron.walk_num, + ctx->electron.coord_new, + 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 + | qmckl_context | context | in | Global state | + | int64_t | elec_num | in | Number of electrons | + | double | kappa_ee | in | Factor to rescale ee distances | + | int64_t | walk_num | in | Number of walkers | + | double | coord[walk_num][3][elec_num] | in | Electron coordinates | + | double | ee_distance[walk_num][elec_num][elec_num] | out | Electron-electron distances | + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer function qmckl_compute_ee_distance_rescaled_f(context, elec_num, 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) :: kappa_ee + integer*8 , intent(in) :: walk_num + double precision , intent(in) :: coord(elec_num,3,walk_num) + 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,1,k), elec_num, & + coord(1,1,k), elec_num, & + ee_distance_rescaled(1,1,k), elec_num, 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 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, 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 :: 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, 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 + +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] : ", np.linalg.norm(elec_1_w1-elec_1_w1) ) +print ( "[0][1][0] : ", np.linalg.norm(elec_1_w1-elec_2_w1) ) +print ( "[1][0][0] : ", np.linalg.norm(elec_2_w1-elec_1_w1) ) +print ( "[0][0][1] : ", np.linalg.norm(elec_1_w2-elec_1_w2) ) +print ( "[0][1][1] : ", np.linalg.norm(elec_1_w2-elec_2_w2) ) +print ( "[1][0][1] : ", np.linalg.norm(elec_2_w2-elec_1_w2) ) + #+end_src + + #+RESULTS: + : [0][0][0] : 0.0 + : [0][1][0] : 7.152322512964209 + : [1][0][0] : 7.152322512964209 + : [0][0][1] : 0.0 + : [0][1][1] : 6.5517646321055665 + : [1][0][1] : 6.5517646321055665 + + #+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); + +// 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-nucleus distances *** Get