#+TITLE: CHAMP Jastrow Factor Single #+SETUPFILE: ../tools/theme.setup #+INCLUDE: ../tools/lib.org * Introduction * Headers :noexport: #+begin_src elisp :noexport :results none (org-babel-lob-ingest "../tools/lib.org") #+end_src #+begin_src c :tangle (eval h_private_func) #ifndef QMCKL_JASTROW_CHAMP_SINGLE_HPF #define QMCKL_JASTROW_CHAMP_SINGLE_HPF #+end_src #+begin_src c :tangle (eval h_private_type) #ifndef QMCKL_JASTROW_CHAMP_SINGLE_HPT #define QMCKL_JASTROW_CHAMP_SINGLE_HPT #include #+end_src #+begin_src c :tangle (eval c_test) :noweb yes #include "qmckl.h" #include #include #ifdef HAVE_CONFIG_H #include "config.h" #endif #include #include #include "n2.h" #include "qmckl_jastrow_champ_private_func.h" #include "qmckl_jastrow_champ_single_private_func.h" int main() { qmckl_context context; context = qmckl_context_create(); #+end_src #+begin_src c :tangle (eval c) #ifdef HAVE_CONFIG_H #include "config.h" #endif #ifdef HAVE_STDINT_H #include #elif HAVE_INTTYPES_H #include #endif #include #include #include #include #include #include #include "qmckl.h" #include "qmckl_context_private_type.h" #include "qmckl_memory_private_type.h" #include "qmckl_memory_private_func.h" #include "qmckl_jastrow_champ_private_type.h" #include "qmckl_jastrow_champ_private_func.h" #include "qmckl_jastrow_champ_single_private_type.h" #include "qmckl_jastrow_champ_single_private_func.h" #+end_src * Context ** Data structure #+begin_src c :comments org :tangle (eval h_private_type) typedef struct qmckl_jastrow_champ_single_struct{ int64_t num; uint64_t date; qmckl_matrix coord; double * restrict een_rescaled_single_e; uint64_t een_rescaled_single_e_date; double * restrict een_rescaled_single_n; uint64_t een_rescaled_single_n_date; uint64_t single_ee_distance_date; double* single_ee_distance; uint64_t single_en_distance_date; double* single_en_distance; double* delta_een; uint64_t delta_een_date; double* delta_p; uint64_t delta_p_date; double* ee_rescaled_single; uint64_t ee_rescaled_single_date; double* en_rescaled_single; uint64_t en_rescaled_single_date; double* delta_en; uint64_t delta_en_date; double* delta_ee; uint64_t delta_ee_date; double * restrict een_rescaled_single_e_gl; uint64_t een_rescaled_single_e_gl_date; double * restrict een_rescaled_single_n_gl; uint64_t een_rescaled_single_n_gl_date; double* delta_p_gl; uint64_t delta_p_gl_date; double* delta_een_gl; uint64_t delta_een_gl_date; double* ee_rescaled_single_gl; uint64_t ee_rescaled_single_gl_date; double* en_rescaled_single_gl; uint64_t en_rescaled_single_gl_date; double* delta_en_gl; uint64_t delta_en_gl_date; double* delta_ee_gl; uint64_t delta_ee_gl_date; } qmckl_jastrow_champ_single_struct; #+end_src ** Test #+begin_src c :tangle (eval c_test) /* Reference input data */ int64_t walk_num = n2_walk_num; int64_t elec_num = n2_elec_num; int64_t elec_up_num = n2_elec_up_num; int64_t elec_dn_num = n2_elec_dn_num; int64_t nucl_num = n2_nucl_num; double rescale_factor_ee = 0.6; double rescale_factor_en[2] = { 0.6, 0.6 }; double* elec_coord = &(n2_elec_coord[0][0][0]); const double* nucl_charge = n2_charge; double* nucl_coord = &(n2_nucl_coord[0][0]); /* Provide Electron data */ qmckl_exit_code rc; assert(!qmckl_electron_provided(context)); rc = qmckl_check(context, qmckl_set_electron_num (context, elec_up_num, elec_dn_num) ); assert(rc == QMCKL_SUCCESS); assert(qmckl_electron_provided(context)); rc = qmckl_check(context, qmckl_set_electron_coord (context, 'N', walk_num, elec_coord, walk_num*3*elec_num) ); assert(rc == QMCKL_SUCCESS); double elec_coord2[walk_num*3*elec_num]; rc = qmckl_check(context, qmckl_get_electron_coord (context, 'N', elec_coord2, walk_num*3*elec_num) ); assert(rc == QMCKL_SUCCESS); for (int64_t i=0 ; i<3*elec_num ; ++i) { assert( elec_coord[i] == elec_coord2[i] ); } /* Provide Nucleus data */ assert(!qmckl_nucleus_provided(context)); rc = qmckl_check(context, qmckl_set_nucleus_num (context, nucl_num) ); assert(rc == QMCKL_SUCCESS); assert(!qmckl_nucleus_provided(context)); double nucl_coord2[3*nucl_num]; rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2, 3*nucl_num); assert(rc == QMCKL_NOT_PROVIDED); rc = qmckl_check(context, qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), 3*nucl_num) ); assert(rc == QMCKL_SUCCESS); assert(!qmckl_nucleus_provided(context)); rc = qmckl_check(context, qmckl_get_nucleus_coord (context, 'N', nucl_coord2, nucl_num*3) ); assert(rc == QMCKL_SUCCESS); for (int64_t k=0 ; k<3 ; ++k) { for (int64_t i=0 ; isingle_point.coord.data != NULL) { rc = qmckl_matrix_free(context, &(ctx->single_point.coord)); assert (rc == QMCKL_SUCCESS); } ctx->single_point.coord = qmckl_matrix_alloc(context, 1, 3); if (ctx->single_point.coord.data == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_set_single_point", NULL); } ctx->single_point.num = num; if (transp == 'T') { double *a = ctx->single_point.coord.data; #ifdef HAVE_OPENMP #pragma omp for #endif for (int64_t i=0 ; i<3 ; ++i) { a[i] = coord[i]; } } else { qmckl_mat(ctx->single_point.coord, 0, 0) = coord[0]; qmckl_mat(ctx->single_point.coord, 0, 1) = coord[1]; qmckl_mat(ctx->single_point.coord, 0, 2) = coord[2]; } /* Increment the date of the context */ ctx->single_point.date += 1UL; return QMCKL_SUCCESS; } qmckl_exit_code qmckl_set_single_point_f (qmckl_context context, const char transp, const int64_t num, const double* coord, const int64_t size_max) { return qmckl_set_single_point(context, transp, num-1, coord, size_max); } #+end_src #+begin_src f90 :comments org :tangle (eval fh_func) :noweb yes interface integer(qmckl_exit_code) function qmckl_set_single_point(context, & transp, num, coord, size_max) bind(C, name="qmckl_set_single_point_f") use, intrinsic :: iso_c_binding import implicit none integer (c_int64_t) , intent(in) , value :: context character(c_char) , intent(in) , value :: transp integer (c_int64_t) , intent(in) , value :: num real (c_double ) , intent(in) :: coord(*) integer (c_int64_t) , intent(in) , value :: size_max end function end interface #+end_src * Electron-electron distances single point ** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_single_electron_ee_distance(qmckl_context context, double* const distance); #+end_src #+begin_src f90 :tangle (eval fh_func) :comments org :exports none interface integer(qmckl_exit_code) function qmckl_get_single_electron_ee_distance(context, distance) & bind(C) use, intrinsic :: iso_c_binding import implicit none integer (c_int64_t) , intent(in) , value :: context real (c_double ) , intent(out) :: distance(*) end function end interface #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_single_electron_ee_distance(qmckl_context context, double* const distance) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_exit_code rc; rc = qmckl_provide_single_ee_distance(context); if (rc != QMCKL_SUCCESS) return rc; qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); size_t sze = ctx->electron.num * ctx->electron.walker.num; memcpy(distance, ctx->single_point.single_ee_distance, 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_single_ee_distance(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_provide_single_ee_distance(qmckl_context context) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); /* Compute if necessary */ if (ctx->single_point.date > ctx->single_point.single_ee_distance_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { free(ctx->single_point.single_ee_distance); ctx->single_point.single_ee_distance = NULL; } /* Allocate array */ if (ctx->single_point.single_ee_distance == NULL) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = ctx->electron.num * ctx->electron.walker.num * sizeof(double); double* single_ee_distance = (double*) qmckl_malloc(context, mem_info); if (single_ee_distance == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_single_ee_distance", NULL); } ctx->single_point.single_ee_distance = single_ee_distance; } qmckl_exit_code rc = qmckl_compute_single_ee_distance(context, ctx->electron.num, ctx->electron.walker.num, ctx->electron.walker.point.coord.data, ctx->single_point.coord.data, ctx->single_point.single_ee_distance); if (rc != QMCKL_SUCCESS) { return rc; } ctx->single_point.single_ee_distance_date = ctx->date; } // printf("single_ee_distance_date %u\n", ctx->single_point.single_ee_distance_date); return QMCKL_SUCCESS; } #+end_src ** Compute :PROPERTIES: :Name: qmckl_compute_single_ee_distance :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_single_ee_distance_args | Variable | Type | In/Out | Description | |----------------------+---------------------------------+--------+-----------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~elec_num~ | ~int64_t~ | in | Number of electrons | | ~walk_num~ | ~int64_t~ | in | Number of walkers | | ~coord~ | ~double[3][walk_num][elec_num]~ | in | Electron coordinates | | ~single_coord~ | ~double[1][3]~ | in | Single electron coordinates | | ~single_ee_distance~ | ~double[walk_num][elec_num]~ | out | Electron-electron distances | #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_single_ee_distance_f(context, elec_num, walk_num, coord, single_coord, single_ee_distance) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: elec_num integer*8 , intent(in) :: walk_num double precision , intent(in) :: coord(elec_num,walk_num,3) double precision , intent(in) :: single_coord(1,3) double precision , intent(out) :: single_ee_distance(elec_num,walk_num) integer*8 :: k, i, j double precision :: x, y, z 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(context, 'T', 'T', 1_8, elec_num, & single_coord(1,1), 1_8, & coord(1,k,1), elec_num * walk_num * 1_8, & single_ee_distance(1,k), 1_8) if (info /= QMCKL_SUCCESS) then exit endif end do end function qmckl_compute_single_ee_distance_f #+end_src #+begin_src c :tangle (eval h_private_func) :comments org :exports none qmckl_exit_code qmckl_compute_single_ee_distance ( const qmckl_context context, const int64_t elec_num, const int64_t walk_num, const double* coord, const double* single_coord, double* const single_ee_distance ); #+end_src #+CALL: generate_c_interface(table=qmckl_single_ee_distance_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_single_ee_distance & (context, elec_num, walk_num, coord, single_coord, single_ee_distance) & 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 integer (c_int64_t) , intent(in) , value :: walk_num real (c_double ) , intent(in) :: coord(elec_num,3,walk_num) real (c_double ) , intent(in) :: single_coord(3) real (c_double ) , intent(out) :: single_ee_distance(elec_num,walk_num) integer(c_int32_t), external :: qmckl_compute_single_ee_distance_f info = qmckl_compute_single_ee_distance_f & (context, elec_num, walk_num, coord, single_coord, single_ee_distance) end function qmckl_compute_single_ee_distance #+end_src ** test #+begin_src c :tangle (eval c_test) /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); rc = qmckl_set_point(context, 'N', elec_num, elec_coord, walk_num*elec_num*3); assert(rc == QMCKL_SUCCESS); double ee_distance[walk_num][elec_num][elec_num]; rc = qmckl_get_electron_ee_distance(context, &ee_distance[0][0][0]); rc = qmckl_get_electron_coord(context, 'N', &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); rc = qmckl_set_single_point(context, 'N', 2, new_coords, 3); assert (rc == QMCKL_SUCCESS); double single_ee_distance[walk_num][elec_num]; rc = qmckl_get_single_electron_ee_distance(context, &single_ee_distance[0][0]); assert (rc == QMCKL_SUCCESS); coords[0][2][0] = new_coords[0]; coords[0][2][1] = new_coords[1]; coords[0][2][2] = new_coords[2]; rc = qmckl_set_point(context, 'N', elec_num, &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_electron_ee_distance(context, &ee_distance[0][0][0]); assert (rc == QMCKL_SUCCESS); for (int nw = 0; nw < walk_num; nw++){ for (int i = 0; i < elec_num; i++) { if (i == 2) continue; assert(fabs((ee_distance[nw][2][i]-single_ee_distance[nw][i])) < 1.e-12); } } #+end_src * ee distance rescaled single point ** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_een_rescaled_single_e(qmckl_context context, double* const distance_rescaled, const int64_t size_max); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_een_rescaled_single_e(qmckl_context context, double* const distance_rescaled, const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_exit_code rc; rc = qmckl_provide_een_rescaled_single_e(context); if (rc != QMCKL_SUCCESS) return rc; qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int64_t sze = ctx->electron.num * ctx->electron.walker.num * (ctx->jastrow_champ.cord_num + 1); if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "todo", "Array too small. Expected ctx->electron.num * ctx->electron.num * ctx->electron.walker.num * (ctx->jastrow_champ.cord_num + 1)"); } memcpy(distance_rescaled, ctx->single_point.een_rescaled_single_e, 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_single_e(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_provide_een_rescaled_single_e(qmckl_context context) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); /* Check if ee distance is provided */ qmckl_exit_code rc = qmckl_provide_single_ee_distance(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if rescaled ee distance is provided */ rc = qmckl_provide_een_rescaled_e(context); if(rc != QMCKL_SUCCESS) return rc; /* Compute if necessary */ if (ctx->single_point.date > ctx->single_point.een_rescaled_single_e_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { if (ctx->single_point.een_rescaled_single_e!= NULL) { rc = qmckl_free(context, ctx->single_point.een_rescaled_single_e); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_een_rescaled_single_e", "Unable to free ctx->single_point.een_rescaled_single_e"); } ctx->single_point.een_rescaled_single_e = NULL; } } /* Allocate array */ if (ctx->single_point.een_rescaled_single_e == NULL) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = ctx->electron.num * ctx->electron.walker.num * (ctx->jastrow_champ.cord_num + 1) * sizeof(double); double* een_rescaled_single_e = (double*) qmckl_malloc(context, mem_info); if (een_rescaled_single_e == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_een_rescaled_single_e", NULL); } ctx->single_point.een_rescaled_single_e = een_rescaled_single_e; } rc = qmckl_compute_een_rescaled_single_e(context, ctx->single_point.num, ctx->electron.walker.num, ctx->electron.num, ctx->jastrow_champ.cord_num, ctx->jastrow_champ.rescale_factor_ee, ctx->single_point.single_ee_distance, ctx->jastrow_champ.een_rescaled_e, ctx->single_point.een_rescaled_single_e); if (rc != QMCKL_SUCCESS) { return rc; } ctx->single_point.een_rescaled_single_e_date = ctx->single_point.date; } return QMCKL_SUCCESS; } #+end_src ** Compute :PROPERTIES: :Name: qmckl_compute_een_rescaled_single_e :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_factor_een_rescaled_single_e_args | Variable | Type | In/Out | Description | |-------------------------+----------------------------------------------------+--------+-------------------------------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~num~ | ~int64_t~ | in | Number of walkers | | ~walk_num~ | ~int64_t~ | in | Number of walkers | | ~elec_num~ | ~int64_t~ | in | Number of electrons | | ~cord_num~ | ~int64_t~ | in | Order of polynomials | | ~rescale_factor_ee~ | ~double~ | in | Factor to rescale ee distances | | ~single_ee_distance~ | ~double[walk_num][elec_num]~ | in | Single electron-electron distances for each walker | | ~een_rescaled_e~ | ~double[walk_num][0:cord_num][elec_num][elec_num]~ | in | Rescaled electron-electron distances for each walker | | ~een_rescaled_single_e~ | ~double[walk_num][0:cord_num][elec_num]~ | out | Single electron-electron rescaled distances for each walker | #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_een_rescaled_single_e_doc_f( & context, num_in, walk_num, elec_num, cord_num, rescale_factor_ee, & single_ee_distance, een_rescaled_e, een_rescaled_single_e) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: num_in integer*8 , intent(in) :: walk_num integer*8 , intent(in) :: elec_num integer*8 , intent(in) :: cord_num double precision , intent(in) :: rescale_factor_ee double precision , intent(in) :: single_ee_distance(elec_num,walk_num) double precision , intent(in) :: een_rescaled_e(elec_num,elec_num,0:cord_num,walk_num) double precision , intent(out) :: een_rescaled_single_e(elec_num,0:cord_num,walk_num) double precision,allocatable :: een_rescaled_single_e_ij(:,:) double precision :: x integer*8 :: i, j, k, l, nw, num num = num_in + 1 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 (cord_num < 0) then info = QMCKL_INVALID_ARG_4 return endif allocate(een_rescaled_single_e_ij(elec_num, cord_num + 1)) ! Prepare table of exponentiated distances raised to appropriate power do nw = 1, walk_num een_rescaled_single_e_ij(:, 1) = 1.0d0 do j = 1, elec_num een_rescaled_single_e_ij(j, 2) = dexp(-rescale_factor_ee * single_ee_distance(j, nw)) end do do l = 2, cord_num do k = 1, elec_num een_rescaled_single_e_ij(k, l + 1) = een_rescaled_single_e_ij(k, l) * een_rescaled_single_e_ij(k, 2) end do end do ! prepare the actual een table een_rescaled_single_e( :, 0, nw) = 1.0d0 do l = 1, cord_num do j = 1, elec_num x = een_rescaled_single_e_ij(j, l + 1) een_rescaled_single_e(j, l, nw) = x end do end do !een_rescaled_single_e(:,:,:) = een_rescaled_single_e(:,:,:) - een_rescaled_e(num,:,:,:) een_rescaled_single_e(num, :, :) = 0.0d0 end do end function qmckl_compute_een_rescaled_single_e_doc_f #+end_src # #+CALL: generate_c_header(table=qmckl_factor_een_rescaled_e_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_een_rescaled_single_e ( const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t cord_num, const double rescale_factor_ee, const double* single_ee_distance, const double* een_rescaled_e, double* const een_rescaled_single_e ); #+end_src #+CALL: generate_c_interface(table=qmckl_factor_een_rescaled_single_e_args,rettyp=get_value("CRetType"),fname="qmckl_compute_een_rescaled_single_e_doc") #+RESULTS: #+begin_src f90 :tangle (eval f) :comments org :exports none integer(c_int32_t) function qmckl_compute_een_rescaled_single_e_doc & (context, num, walk_num, elec_num, cord_num, rescale_factor_ee, & single_ee_distance, een_rescaled_e, een_rescaled_single_e) & 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 :: num 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 :: cord_num real (c_double ) , intent(in) , value :: rescale_factor_ee real (c_double ) , intent(in) :: single_ee_distance(elec_num,walk_num) real (c_double ) , intent(in) :: een_rescaled_e(elec_num,elec_num,0:cord_num,walk_num) real (c_double ) , intent(out) :: een_rescaled_single_e(elec_num,0:cord_num,walk_num) integer(c_int32_t), external :: qmckl_compute_een_rescaled_single_e_doc_f info = qmckl_compute_een_rescaled_single_e_doc_f & (context, num, walk_num, elec_num, cord_num, rescale_factor_ee, single_ee_distance, een_rescaled_e, een_rescaled_single_e) end function qmckl_compute_een_rescaled_single_e_doc #+end_src #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_een_rescaled_single_e ( const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t cord_num, const double rescale_factor_ee, const double* single_ee_distance, const double* een_rescaled_e, double* const een_rescaled_single_e ); #+end_src #+begin_src c :tangle (eval h_private_func) :comments org qmckl_exit_code qmckl_compute_een_rescaled_single_e_doc ( const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t cord_num, const double rescale_factor_ee, const double* single_ee_distance, const double* een_rescaled_e, double* const een_rescaled_single_e ); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes qmckl_exit_code qmckl_compute_een_rescaled_single_e ( const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t cord_num, const double rescale_factor_ee, const double* single_ee_distance, const double* een_rescaled_e, double* const een_rescaled_single_e ) { #ifdef HAVE_HPC return qmckl_compute_een_rescaled_single_e_doc #else return qmckl_compute_een_rescaled_single_e_doc #endif (context, num, walk_num, elec_num, cord_num, rescale_factor_ee, single_ee_distance, een_rescaled_e, een_rescaled_single_e); } #+end_src ** test #+begin_src c :tangle (eval c_test) /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); rc = qmckl_set_point(context, 'N', elec_num, elec_coord, walk_num*elec_num*3); assert(rc == QMCKL_SUCCESS); double rescaled_een_ee_distance[walk_num][cord_num+1][elec_num][elec_num]; rc = qmckl_get_jastrow_champ_een_rescaled_e(context, &rescaled_een_ee_distance[0][0][0][0], walk_num*(cord_num+1)*elec_num*elec_num); rc = qmckl_get_electron_coord(context, 'N', &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); rc = qmckl_set_single_point(context, 'N', 2, new_coords, 3); assert (rc == QMCKL_SUCCESS); double single_rescaled_een_ee_distance[walk_num][cord_num+1][elec_num]; rc = qmckl_get_een_rescaled_single_e(context, &single_rescaled_een_ee_distance[0][0][0], walk_num*(cord_num+1)*elec_num); assert (rc == QMCKL_SUCCESS); coords[0][2][0] = new_coords[0]; coords[0][2][1] = new_coords[1]; coords[0][2][2] = new_coords[2]; rc = qmckl_set_point(context, 'N', elec_num, &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_jastrow_champ_een_rescaled_e(context, &rescaled_een_ee_distance[0][0][0][0], walk_num*(cord_num+1)*elec_num*elec_num); assert (rc == QMCKL_SUCCESS); for (int nw = 0; nw < walk_num; nw++){ for (int l = 0; l <= cord_num; l++){ for (int i = 0; i < elec_num; i++) { if (i == 2) continue; assert(fabs((rescaled_een_ee_distance[nw][l][2][i]-single_rescaled_een_ee_distance[nw][l][i])) < 1.e-12); } } } #+end_src * Electron-nucleus distances single point ** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_single_electron_en_distance(qmckl_context context, double* distance); #+end_src #+begin_src f90 :tangle (eval fh_func) :comments org :exports none interface integer(qmckl_exit_code) function qmckl_get_single_electron_en_distance(context, distance) & bind(C) use, intrinsic :: iso_c_binding import implicit none integer (c_int64_t) , intent(in) , value :: context real (c_double ) , intent(out) :: distance(*) end function end interface #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_single_electron_en_distance(qmckl_context context, double* distance) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_exit_code rc; rc = qmckl_provide_single_en_distance(context); if (rc != QMCKL_SUCCESS) return rc; qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); size_t sze = ctx->nucleus.num; memcpy(distance, ctx->single_point.single_en_distance, 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_single_en_distance(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_provide_single_en_distance(qmckl_context context) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); if (!(ctx->nucleus.provided)) { return qmckl_failwith( context, QMCKL_NOT_PROVIDED, "qmckl_provide_single_en_distance", NULL); } /* Compute if necessary */ if (ctx->single_point.date > ctx->single_point.single_en_distance_date) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = ctx->nucleus.num * sizeof(double); if (ctx->single_point.single_en_distance != NULL) { qmckl_memory_info_struct mem_info_test = qmckl_memory_info_struct_zero; qmckl_exit_code rc = qmckl_get_malloc_info(context, ctx->single_point.single_en_distance, &mem_info_test); /* if rc != QMCKL_SUCCESS, we are maybe in an _inplace function because the memory was not allocated with qmckl_malloc */ if ((rc == QMCKL_SUCCESS) && (mem_info_test.size != mem_info.size)) { rc = qmckl_free(context, ctx->single_point.single_en_distance); assert (rc == QMCKL_SUCCESS); ctx->single_point.single_en_distance = NULL; } } /* Allocate array */ if (ctx->single_point.single_en_distance == NULL) { double* single_en_distance = (double*) qmckl_malloc(context, mem_info); if (single_en_distance == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_single_en_distance", NULL); } ctx->single_point.single_en_distance = single_en_distance; } qmckl_exit_code rc = qmckl_compute_single_en_distance(context, ctx->nucleus.num, ctx->single_point.coord.data, ctx->nucleus.coord.data, ctx->single_point.single_en_distance); if (rc != QMCKL_SUCCESS) { return rc; } ctx->single_point.single_en_distance_date = ctx->single_point.date; } // printf("single_en_distance_date %u\n", ctx->single_point.single_en_distance_date); return QMCKL_SUCCESS; } #+end_src ** Compute :PROPERTIES: :Name: qmckl_compute_single_en_distance :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_single_en_distance_args | Variable | Type | In/Out | Description | |----------------------+-----------------------+--------+----------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | | ~elec_coord~ | ~double[3]~ | in | Electron coordinates | | ~nucl_coord~ | ~double[3][nucl_num]~ | in | Nuclear coordinates | | ~single_en_distance~ | ~double[nucl_num]~ | out | Electron-nucleus distances | #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_single_en_distance_f(context, nucl_num, elec_coord, nucl_coord, single_en_distance) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: nucl_num double precision , intent(in) :: elec_coord(3) double precision , intent(in) :: nucl_coord(nucl_num,3) double precision , intent(out) :: single_en_distance(nucl_num) integer*8 :: k info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif if (nucl_num <= 0) then info = QMCKL_INVALID_ARG_2 return endif info = qmckl_distance(context, 'T', 'T', nucl_num, 1_8, & nucl_coord, nucl_num, & elec_coord, 1_8, & single_en_distance, nucl_num) end function qmckl_compute_single_en_distance_f #+end_src #+begin_src c :tangle (eval h_private_func) :comments org :exports none qmckl_exit_code qmckl_compute_single_en_distance ( const qmckl_context context, const int64_t nucl_num, const double* elec_coord, const double* nucl_coord, double* const single_en_distance ); #+end_src #+CALL: generate_c_interface(table=qmckl_single_en_distance_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_single_en_distance & (context, nucl_num, elec_coord, nucl_coord, single_en_distance) & 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 :: nucl_num real (c_double ) , intent(in) :: elec_coord(3) real (c_double ) , intent(in) :: nucl_coord(nucl_num,3) real (c_double ) , intent(out) :: single_en_distance(nucl_num) integer(c_int32_t), external :: qmckl_compute_single_en_distance_f info = qmckl_compute_single_en_distance_f & (context, nucl_num, elec_coord, nucl_coord, single_en_distance) end function qmckl_compute_single_en_distance #+end_src ** test #+begin_src c :tangle (eval c_test) /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); rc = qmckl_set_point(context, 'N', elec_num, elec_coord, walk_num*elec_num*3); assert(rc == QMCKL_SUCCESS); double en_distance[elec_num][nucl_num]; rc = qmckl_get_electron_en_distance(context, &en_distance[0][0]); rc = qmckl_get_electron_coord(context, 'N', &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); rc = qmckl_set_single_point(context, 'N', 2, new_coords, 3); assert (rc == QMCKL_SUCCESS); double single_en_distance[nucl_num]; rc = qmckl_get_single_electron_en_distance(context, &single_en_distance[0]); assert (rc == QMCKL_SUCCESS); coords[0][2][0] = new_coords[0]; coords[0][2][1] = new_coords[1]; coords[0][2][2] = new_coords[2]; rc = qmckl_set_point(context, 'N', elec_num, &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_electron_en_distance(context, &en_distance[0][0]); assert (rc == QMCKL_SUCCESS); for (int a = 0; a < nucl_num; a++){ assert(fabs((en_distance[2][a]-single_en_distance[a])) < 1.e-12); } #+end_src * en distance rescaled single point ** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_een_rescaled_single_n(qmckl_context context, double* const distance_rescaled, const int64_t size_max); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_een_rescaled_single_n(qmckl_context context, double* const distance_rescaled, const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_exit_code rc; rc = qmckl_provide_een_rescaled_single_n(context); if (rc != QMCKL_SUCCESS) return rc; qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int64_t sze = ctx->nucleus.num * ctx->electron.walker.num * (ctx->jastrow_champ.cord_num + 1); if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "mckl_get_een_rescaled_single_n", "Array too small. Expected ctx->nucleus.num * ctx->electron.walker.num * (ctx->jastrow_champ.cord_num + 1)"); } memcpy(distance_rescaled, ctx->single_point.een_rescaled_single_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_single_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_single_n(qmckl_context context) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); /* Check if ee distance is provided */ qmckl_exit_code rc = qmckl_provide_single_en_distance(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if een rescaled distance is provided */ rc = qmckl_provide_een_rescaled_n(context); if(rc != QMCKL_SUCCESS) return rc; /* Compute if necessary */ if (ctx->single_point.date > ctx->single_point.een_rescaled_single_n_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { if (ctx->single_point.een_rescaled_single_n != NULL) { rc = qmckl_free(context, ctx->single_point.een_rescaled_single_n); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_een_rescaled_single_n", "Unable to free ctx->single_point.een_rescaled_single_n"); } ctx->single_point.een_rescaled_single_n = NULL; } } /* Allocate array */ if (ctx->single_point.een_rescaled_single_n == NULL) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = ctx->nucleus.num * ctx->electron.walker.num * (ctx->jastrow_champ.cord_num + 1) * sizeof(double); double* een_rescaled_single_n = (double*) qmckl_malloc(context, mem_info); if (een_rescaled_single_n == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_een_rescaled_single_n", NULL); } ctx->single_point.een_rescaled_single_n = een_rescaled_single_n; } rc = qmckl_compute_een_rescaled_single_n(context, ctx->single_point.num, ctx->electron.walker.num, ctx->electron.num, ctx->nucleus.num, ctx->jastrow_champ.type_nucl_num, ctx->jastrow_champ.type_nucl_vector, ctx->jastrow_champ.cord_num, ctx->jastrow_champ.rescale_factor_en, ctx->single_point.single_en_distance, ctx->jastrow_champ.een_rescaled_n, ctx->single_point.een_rescaled_single_n); if (rc != QMCKL_SUCCESS) { return rc; } ctx->single_point.een_rescaled_single_n_date = ctx->single_point.date; } return QMCKL_SUCCESS; } #+end_src ** Compute :PROPERTIES: :Name: qmckl_compute_een_rescaled_single_n :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_factor_een_rescaled_single_n_args | Variable | Type | In/Out | Description | |-------------------------+----------------------------------------------------+--------+--------------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~num~ | ~int64_t~ | in | Number of single electron | | ~walk_num~ | ~int64_t~ | in | Number of walkers | | ~elec_num~ | ~int64_t~ | in | Number of atoms | | ~nucl_num~ | ~int64_t~ | in | Number of atoms | | ~type_nucl_num~ | ~int64_t~ | in | Number of atom types | | ~type_nucl_vector~ | ~int64_t[nucl_num]~ | in | Types of atoms | | ~cord_num~ | ~int64_t~ | in | Order of polynomials | | ~rescale_factor_en~ | ~double[nucl_num]~ | in | Factor to rescale ee distances | | ~single_en_distance~ | ~double[walk_num][nucl_num]~ | in | Electron-nucleus distances | | ~een_rescaled_n~ | ~double[walk_num][0:cord_num][nucl_num][elec_num]~ | in | Electron-nucleus rescaled distances | | ~een_rescaled_single_n~ | ~double[walk_num][0:cord_num][nucl_num]~ | out | Single electron-nucleus rescaled distances | #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_een_rescaled_single_n_f( & context, num_in, walk_num, elec_num, nucl_num, & type_nucl_num, type_nucl_vector, cord_num, rescale_factor_en, & single_en_distance, een_rescaled_n, een_rescaled_single_n) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: num_in integer*8 , intent(in) :: walk_num integer*8 , intent(in) :: elec_num integer*8 , intent(in) :: nucl_num integer*8 , intent(in) :: type_nucl_num 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) :: single_en_distance(nucl_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_single_n(nucl_num,0:cord_num,walk_num) double precision :: x integer*8 :: i, a, k, l, nw, num num = num_in + 1 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 (nucl_num <= 0) then info = QMCKL_INVALID_ARG_3 return endif if (cord_num < 0) then info = QMCKL_INVALID_ARG_4 return endif do nw = 1, walk_num ! prepare the actual een table een_rescaled_single_n(:, 0, nw) = 1.0d0 do a = 1, nucl_num een_rescaled_single_n(a, 1, nw) = dexp(-rescale_factor_en(type_nucl_vector(a)+1) * single_en_distance(a, nw)) end do do l = 2, cord_num do a = 1, nucl_num een_rescaled_single_n(a, l, nw) = een_rescaled_single_n(a, l - 1, nw) * een_rescaled_single_n(a, 1, nw) end do end do !een_rescaled_single_n(:,:,:) = een_rescaled_single_n(:,:,:) - een_rescaled_n(num,:,:,:) end do end function qmckl_compute_een_rescaled_single_n_f #+end_src #+CALL: generate_c_interface(table=qmckl_factor_een_rescaled_single_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_single_n & (context, & num, & walk_num, & elec_num, & nucl_num, & type_nucl_num, & type_nucl_vector, & cord_num, & rescale_factor_en, & single_en_distance, & een_rescaled_n, & een_rescaled_single_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 :: num 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 :: type_nucl_num integer (c_int64_t) , intent(in) :: type_nucl_vector(nucl_num) integer (c_int64_t) , intent(in) , value :: cord_num real (c_double ) , intent(in) :: rescale_factor_en(nucl_num) real (c_double ) , intent(in) :: single_en_distance(nucl_num,walk_num) real (c_double ) , intent(in) :: een_rescaled_n(elec_num,nucl_num,0:cord_num,walk_num) real (c_double ) , intent(out) :: een_rescaled_single_n(nucl_num,0:cord_num,walk_num) integer(c_int32_t), external :: qmckl_compute_een_rescaled_single_n_f info = qmckl_compute_een_rescaled_single_n_f & (context, & num, & walk_num, & elec_num, & nucl_num, & type_nucl_num, & type_nucl_vector, & cord_num, & rescale_factor_en, & single_en_distance, & een_rescaled_n, & een_rescaled_single_n) end function qmckl_compute_een_rescaled_single_n #+end_src #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_een_rescaled_single_n ( const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t type_nucl_num, int64_t* const type_nucl_vector, const int64_t cord_num, const double* rescale_factor_en, const double* single_en_distance, const double* een_rescaled_n, double* const een_rescaled_single_n ); #+end_src ** test #+begin_src c :tangle (eval c_test) /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); rc = qmckl_set_point(context, 'N', elec_num, elec_coord, walk_num*elec_num*3); assert(rc == QMCKL_SUCCESS); double rescaled_een_en_distance[walk_num][cord_num+1][nucl_num][elec_num]; rc = qmckl_get_jastrow_champ_een_rescaled_n(context, &rescaled_een_en_distance[0][0][0][0], walk_num*(cord_num+1)*nucl_num*elec_num); rc = qmckl_get_electron_coord(context, 'N', &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); rc = qmckl_set_single_point(context, 'N', 2, new_coords, 3); assert (rc == QMCKL_SUCCESS); double single_rescaled_een_en_distance[walk_num][cord_num+1][nucl_num]; rc = qmckl_get_een_rescaled_single_n(context, &single_rescaled_een_en_distance[0][0][0], walk_num*(cord_num+1)*nucl_num); assert (rc == QMCKL_SUCCESS); coords[0][2][0] = new_coords[0]; coords[0][2][1] = new_coords[1]; coords[0][2][2] = new_coords[2]; rc = qmckl_set_point(context, 'N', elec_num, &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_jastrow_champ_een_rescaled_n(context, &rescaled_een_en_distance[0][0][0][0], walk_num*(cord_num+1)*nucl_num*elec_num); assert (rc == QMCKL_SUCCESS); for (int nw = 0; nw < walk_num; nw++){ for (int l = 0; l <= cord_num; l++){ for (int a = 0; a < nucl_num; a++) { assert(fabs((rescaled_een_en_distance[nw][l][a][2]-single_rescaled_een_en_distance[nw][l][a])) < 1.e-12); } } } #+end_src * Delta p ** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_jastrow_champ_delta_p(qmckl_context context, double* const delta_p, const int64_t size_max); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_jastrow_champ_delta_p(qmckl_context context, double* const delta_p, const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_exit_code rc; rc = qmckl_provide_jastrow_champ_delta_p(context); if (rc != QMCKL_SUCCESS) return rc; qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int64_t sze = ctx->electron.walker.num * ctx->jastrow_champ.cord_num * (ctx->jastrow_champ.cord_num + 1) * ctx->nucleus.num * ctx->electron.num; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_jastrow_champ_delta_p", "Array too small."); } memcpy(delta_p, ctx->single_point.delta_p, sze * sizeof(double)); return QMCKL_SUCCESS; } #+end_src ** provide #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_provide_jastrow_champ_delta_p(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_provide_jastrow_champ_delta_p(qmckl_context context) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); qmckl_exit_code rc = qmckl_provide_een_rescaled_single_e(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_een_rescaled_single_n(context); if(rc != QMCKL_SUCCESS) return rc; /* Compute if necessary */ if (ctx->single_point.date > ctx->single_point.delta_p_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { if (ctx->single_point.delta_p != NULL) { rc = qmckl_free(context, ctx->single_point.delta_p); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_jastrow_champ_delta_p", "Unable to free ctx->single_point.delta_p"); } ctx->single_point.delta_p = NULL; } } /* Allocate array */ if (ctx->single_point.delta_p == NULL) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = ctx->electron.walker.num * ctx->jastrow_champ.cord_num * (ctx->jastrow_champ.cord_num + 1) * ctx->nucleus.num * ctx->electron.num * sizeof(double); double* delta_p = (double*) qmckl_malloc(context, mem_info); if (delta_p == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_jastrow_champ_delta_p", NULL); } ctx->single_point.delta_p = delta_p; } rc = qmckl_compute_jastrow_champ_delta_p_doc(context, ctx->single_point.num, ctx->electron.walker.num, ctx->electron.num, ctx->nucleus.num, ctx->jastrow_champ.cord_num, ctx->jastrow_champ.een_rescaled_n, ctx->jastrow_champ.een_rescaled_e, ctx->single_point.een_rescaled_single_n, ctx->single_point.een_rescaled_single_e, ctx->single_point.delta_p); if (rc != QMCKL_SUCCESS) { return rc; } ctx->single_point.delta_p_date = ctx->single_point.date; } return QMCKL_SUCCESS; } #+end_src ** Compute :PROPERTIES: :Name: qmckl_compute_jastrow_champ_delta_p_doc :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_factor_delta_p_args | Variable | Type | In/Out | Description | |-------------------------+------------------------------------------------------------------+--------+---------------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~num~ | ~int64_t~ | in | Single point number | | ~walk_num~ | ~int64_t~ | in | Number of walkers | | ~elec_num~ | ~int64_t~ | in | Number of electrons | | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | | ~cord_num~ | ~int64_t~ | in | order of polynomials | | ~een_rescaled_n~ | ~double[walk_num][0:cord_num][nucl_num][elec_num]~ | in | Electron-nucleus rescaled distances | | ~een_rescaled_e~ | ~double[walk_num][0:cord_num][elec_num][elec_num]~ | in | Electron-electron rescaled distances | | ~een_rescaled_single_n~ | ~double[walk_num][0:cord_num][nucl_num]~ | in | Electron-nucleus single rescaled distances | | ~een_rescaled_single_e~ | ~double[walk_num][0:cord_num][elec_num]~ | in | Electron-electron single rescaled distances | | ~delta_p~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]~ | out | Electron-nucleus jastrow | #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_jastrow_champ_delta_p_doc_f( & context, num_in, walk_num, elec_num, nucl_num, cord_num, & een_rescaled_n, een_rescaled_e, een_rescaled_single_n, een_rescaled_single_e, delta_p) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: num_in, walk_num, elec_num, cord_num, nucl_num double precision , intent(in) :: een_rescaled_n(elec_num, nucl_num, 0:cord_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_single_n(nucl_num, 0:cord_num, walk_num) double precision , intent(in) :: een_rescaled_single_e(elec_num, 0:cord_num, walk_num) double precision , intent(out) :: delta_p(elec_num, nucl_num,0:cord_num, 0:cord_num-1, walk_num) double precision :: delta_c(nucl_num,0:cord_num, 0:cord_num-1, walk_num) double precision :: delta_c2(elec_num, nucl_num,0:cord_num, 0:cord_num-1, walk_num) double precision :: een_rescaled_delta_n(nucl_num, 0:cord_num, walk_num) double precision :: een_rescaled_delta_e(elec_num, 0:cord_num, walk_num) integer*8 :: i, a, j, l, k, p, m, n, nw, num double precision :: accu, accu2, cn integer*8 :: LDA, LDB, LDC num = num_in + 1 info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT if (walk_num <= 0) info = QMCKL_INVALID_ARG_3 if (elec_num <= 0) info = QMCKL_INVALID_ARG_4 if (nucl_num <= 0) info = QMCKL_INVALID_ARG_5 if (cord_num < 0) info = QMCKL_INVALID_ARG_6 if (info /= QMCKL_SUCCESS) return if (cord_num == 0) return een_rescaled_delta_n(:,:,:) = een_rescaled_single_n(:,:,:) - een_rescaled_n(num,:,:,:) een_rescaled_delta_e(:,:,:) = een_rescaled_single_e(:,:,:) - een_rescaled_e(:,num,:,:) ! First calculate delta p do nw=1, walk_num do i=0, cord_num-1 info = qmckl_dgemm(context, 'T', 'N', 1_8, nucl_num * (cord_num+1_8), elec_num, 1.0d0, & een_rescaled_delta_e(1,i,nw),elec_num, & een_rescaled_n(1,1,0,nw),elec_num, & 0.0d0, & delta_c(1,0,i,nw),1_8) info = qmckl_dgemm(context, 'N', 'N', elec_num, nucl_num * (cord_num+1), 1_8, 1.0d0, & een_rescaled_delta_e(1,i,nw),elec_num, & een_rescaled_n(num,1,0,nw),elec_num, & 0.0d0, & delta_c2(1,1,0,i,nw),elec_num) delta_c2(num,:,:,i,nw) = delta_c2(num,:,:,i,nw) + delta_c(:,:,i,nw) delta_p(:,:,:,i,nw) = delta_c2(:,:,:,i,nw) info = qmckl_dgemm(context, 'N', 'T', elec_num, nucl_num * (cord_num+1), 1_8, 1.0d0, & een_rescaled_e(:,num,i,nw),elec_num, & een_rescaled_delta_n(:,:,nw),nucl_num* (cord_num+1), & 0.0d0, & delta_c2(:,:,:,i,nw),elec_num) delta_p(:,:,:,i,nw) = delta_p(:,:,:,i,nw) + delta_c2(:,:,:,i,nw) info = qmckl_dgemm(context, 'N', 'T', elec_num, nucl_num * (cord_num+1), 1_8, 1.0d0, & een_rescaled_delta_e(:,i,nw),elec_num, & een_rescaled_delta_n(:,:,nw),nucl_num* (cord_num+1), & 0.0d0, & delta_c2(:,:,:,i,nw),elec_num) delta_p(:,:,:,i,nw) = delta_p(:,:,:,i,nw) + delta_c2(:,:,:,i,nw) end do end do end function qmckl_compute_jastrow_champ_delta_p_doc_f #+end_src # #+CALL: generate_c_header(table=qmckl_factor_een_args,rettyp=qmckl_exit_code),fname=get_value("Name")) #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_jastrow_champ_delta_p_doc (const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, const double* een_rescaled_n, const double* een_rescaled_e, const double* een_rescaled_single_n, const double* een_rescaled_single_e, double* const delta_p ); qmckl_exit_code qmckl_compute_jastrow_champ_delta_p (const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, const double* een_rescaled_n, const double* een_rescaled_e, const double* een_rescaled_single_n, const double* een_rescaled_single_e, double* const delta_p ); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_compute_jastrow_champ_delta_p (const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, const double* een_rescaled_n, const double* een_rescaled_e, const double* een_rescaled_single_n, const double* een_rescaled_single_e, double* const delta_p ) { #ifdef HAVE_HPC return qmckl_compute_jastrow_champ_delta_p_doc #else return qmckl_compute_jastrow_champ_delta_p_doc #endif (context, num, walk_num, elec_num, nucl_num, cord_num, een_rescaled_n, een_rescaled_e, een_rescaled_single_n, een_rescaled_single_e, delta_p ); } #+end_src #+CALL: generate_c_interface(table=qmckl_factor_delta_p_args,rettyp=get_value("CRetType"),fname="qmckl_compute_jastrow_champ_delta_p_doc")) #+RESULTS: #+begin_src f90 :tangle (eval f) :comments org :exports none integer(c_int32_t) function qmckl_compute_jastrow_champ_delta_p_doc & (context, & num, & walk_num, & elec_num, & nucl_num, & cord_num, & een_rescaled_n, & een_rescaled_e, & een_rescaled_single_n, & een_rescaled_single_e, & delta_p) & 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 :: num 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) :: een_rescaled_n(elec_num,nucl_num,0:cord_num,walk_num) real (c_double ) , intent(in) :: een_rescaled_e(elec_num,elec_num,0:cord_num,walk_num) real (c_double ) , intent(in) :: een_rescaled_single_n(nucl_num,0:cord_num,walk_num) real (c_double ) , intent(in) :: een_rescaled_single_e(elec_num,0:cord_num,walk_num) real (c_double ) , intent(out) :: delta_p(elec_num,nucl_num,0:cord_num,0:cord_num-1,walk_num) integer(c_int32_t), external :: qmckl_compute_jastrow_champ_delta_p_doc_f info = qmckl_compute_jastrow_champ_delta_p_doc_f & (context, & num, & walk_num, & elec_num, & nucl_num, & cord_num, & een_rescaled_n, & een_rescaled_e, & een_rescaled_single_n, & een_rescaled_single_e, & delta_p) end function qmckl_compute_jastrow_champ_delta_p_doc #+end_src ** test #+begin_src c :tangle (eval c_test) /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); rc = qmckl_set_point(context, 'N', elec_num, elec_coord, walk_num*elec_num*3); assert(rc == QMCKL_SUCCESS); double p_old[walk_num][cord_num][cord_num+1][nucl_num][elec_num]; rc = qmckl_get_jastrow_champ_tmp_c(context, &p_old[0][0][0][0][0]); rc = qmckl_get_electron_coord(context, 'N', &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); rc = qmckl_set_single_point(context, 'N', 2, new_coords, 3); assert (rc == QMCKL_SUCCESS); double delta_p[walk_num][cord_num][cord_num+1][nucl_num][elec_num]; rc = qmckl_get_jastrow_champ_delta_p(context, &delta_p[0][0][0][0][0], walk_num*cord_num*(cord_num+1)*nucl_num*elec_num); assert (rc == QMCKL_SUCCESS); coords[0][2][0] = new_coords[0]; coords[0][2][1] = new_coords[1]; coords[0][2][2] = new_coords[2]; rc = qmckl_set_point(context, 'N', elec_num, &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); double p_new[walk_num][cord_num][cord_num+1][nucl_num][elec_num]; rc = qmckl_get_jastrow_champ_tmp_c(context, &p_new[0][0][0][0][0]); assert (rc == QMCKL_SUCCESS); for (int nw = 0; nw < walk_num; nw++){ for (int l = 0; l < cord_num; l++){ for (int m = 0; m <= cord_num; m++){ for (int a = 0; a < nucl_num; a++) { for (int i = 0; i < elec_num; i++){ assert(fabs(((p_new[nw][l][m][a][i]-p_old[nw][l][m][a][i])-delta_p[nw][l][m][a][i])) < 1.e-12); } } } } } #+end_src * Delta e-e-n ** get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_jastrow_champ_single_een(qmckl_context context, double* const delta_een, const int64_t size_max); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_jastrow_champ_single_een(qmckl_context context, double* const delta_een, const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_exit_code rc; rc = qmckl_provide_jastrow_champ_single_een(context); if (rc != QMCKL_SUCCESS) return rc; qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int64_t sze = ctx->electron.walker.num; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_jastrow_champ_single_eenn", "Array too small. Expected ctx->electron.walker.num"); } memcpy(delta_een, ctx->single_point.delta_een, sze * sizeof(double)); return QMCKL_SUCCESS; } #+end_src #+begin_src f90 :tangle (eval fh_func) :comments org interface integer(qmckl_exit_code) function qmckl_get_jastrow_champ_single_een (context, & delta_een, size_max) bind(C) use, intrinsic :: iso_c_binding import implicit none integer (qmckl_context) , intent(in), value :: context integer(c_int64_t), intent(in), value :: size_max real(c_double), intent(out) :: delta_een(size_max) end function end interface #+end_src ** Provide :noexport: #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_provide_jastrow_champ_single_een(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_provide_jastrow_champ_single_een(qmckl_context context) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); qmckl_exit_code rc = qmckl_provide_jastrow_champ_delta_p(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_jastrow_champ_factor_een(context); if(rc != QMCKL_SUCCESS) return rc; /* Compute if necessary */ if (ctx->single_point.date > ctx->single_point.delta_een_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { if (ctx->single_point.delta_een != NULL) { rc = qmckl_free(context, ctx->single_point.delta_een); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_jastrow_champ_single_een", "Unable to free ctx->single_point.delta_een"); } ctx->single_point.delta_een = NULL; } } /* Allocate array */ if (ctx->single_point.delta_een == NULL) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = ctx->electron.walker.num * sizeof(double); double* delta_een = (double*) qmckl_malloc(context, mem_info); if (delta_een == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_jastrow_champ_single_een", NULL); } ctx->single_point.delta_een = delta_een; } rc = qmckl_compute_jastrow_champ_factor_single_een_doc(context, ctx->single_point.num, 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.tmp_c, ctx->single_point.delta_p, ctx->jastrow_champ.een_rescaled_n, ctx->jastrow_champ.een_rescaled_e, ctx->single_point.een_rescaled_single_n, ctx->single_point.een_rescaled_single_e, ctx->single_point.delta_een); if (rc != QMCKL_SUCCESS) { return rc; } ctx->single_point.delta_een_date = ctx->single_point.date; } //printf("ctx date %u\n", ctx->date); //printf("single point date %u\n", ctx->single_point.date); //printf("jastrow champ tmp_c date %u\n", ctx->jastrow_champ.tmp_c_date); //printf("delta p date %u\n", ctx->single_point.delta_p_date); return QMCKL_SUCCESS; } #+end_src ** Compute :PROPERTIES: :Name: qmckl_compute_jastrow_champ_factor_single_een_doc :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_factor_single_een_args | Variable | Type | In/Out | Description | |-------------------------+------------------------------------------------------------------+--------+---------------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~num~ | ~int64_t~ | in | Single point number | | ~walk_num~ | ~int64_t~ | in | Number of walkers | | ~elec_num~ | ~int64_t~ | in | Number of electrons | | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | | ~cord_num~ | ~int64_t~ | in | order of polynomials | | ~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 | | ~tmp_c~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]~ | in | vector of non-zero coefficients | | ~delta_p~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]~ | in | vector of non-zero coefficients | | ~een_rescaled_n~ | ~double[walk_num][0:cord_num][nucl_num][elec_num]~ | in | Electron-nucleus rescaled distances | | ~een_rescaled_e~ | ~double[walk_num][0:cord_num][elec_num][elec_num]~ | in | Electron-electron rescaled distances | | ~een_rescaled_single_n~ | ~double[walk_num][0:cord_num][nucl_num]~ | in | Electron-nucleus single rescaled distances | | ~een_rescaled_single_e~ | ~double[walk_num][0:cord_num][elec_num]~ | in | Electron-electron single rescaled distances | | ~delta_een~ | ~double[walk_num]~ | out | Electron-nucleus jastrow | #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_jastrow_champ_factor_single_een_doc_f( & context, num_in, walk_num, elec_num, nucl_num, cord_num, & dim_c_vector, c_vector_full, lkpm_combined_index, & tmp_c, delta_p, een_rescaled_n, een_rescaled_e, een_rescaled_single_n, & een_rescaled_single_e, delta_een) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: num_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) :: tmp_c(elec_num, nucl_num,0:cord_num, 0:cord_num-1, walk_num) double precision , intent(in) :: delta_p(elec_num, nucl_num,0:cord_num, 0:cord_num-1, walk_num) double precision , intent(in) :: een_rescaled_n(elec_num, nucl_num, 0:cord_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_single_n(nucl_num, 0:cord_num, walk_num) double precision , intent(in) :: een_rescaled_single_e(elec_num, 0:cord_num, walk_num) double precision , intent(out) :: delta_een(walk_num) double precision :: delta_c(nucl_num,0:cord_num, 0:cord_num-1, walk_num) double precision :: delta_c2(elec_num, nucl_num,0:cord_num, 0:cord_num-1, walk_num) double precision :: een_rescaled_delta_n(nucl_num, 0:cord_num, walk_num) integer*8 :: i, a, j, l, k, p, m, n, nw, num double precision :: accu, accu2, cn integer*8 :: LDA, LDB, LDC num = num_in + 1 info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT if (walk_num <= 0) info = QMCKL_INVALID_ARG_3 if (elec_num <= 0) info = QMCKL_INVALID_ARG_4 if (nucl_num <= 0) info = QMCKL_INVALID_ARG_5 if (cord_num < 0) info = QMCKL_INVALID_ARG_6 if (info /= QMCKL_SUCCESS) return delta_een = 0.0d0 een_rescaled_delta_n(:,:,:) = een_rescaled_single_n(:,:,:) - een_rescaled_n(num,:,:,:) if (cord_num == 0) return 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) do a = 1, nucl_num cn = c_vector_full(a, n) if(cn == 0.d0) cycle accu = 0.0d0 do j = 1, elec_num accu = accu + een_rescaled_n(j,a,m,nw) * delta_p(j,a,m+l,k,nw) end do accu = accu + een_rescaled_delta_n(a,m,nw) * (tmp_c(num,a,m+l,k,nw) + delta_p(num,a,m+l,k,nw)) delta_een(nw) = delta_een(nw) + accu * cn end do end do end do end function qmckl_compute_jastrow_champ_factor_single_een_doc_f #+end_src # #+CALL: generate_c_header(table=qmckl_factor_een_args,rettyp=qmckl_exit_code),fname=get_value("Name")) #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_jastrow_champ_factor_single_een_doc (const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, const int64_t dim_c_vector, const double* c_vector_full, const int64_t* lkpm_combined_index, const double* tmp_c, const double* delta_p, const double* een_rescaled_n, const double* een_rescaled_e, const double* een_rescaled_single_n, const double* een_rescaled_single_e, double* const delta_een ); qmckl_exit_code qmckl_compute_jastrow_champ_factor_single_een (const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, const int64_t dim_c_vector, const double* c_vector_full, const int64_t* lkpm_combined_index, const double* tmp_c, const double* delta_p, const double* een_rescaled_n, const double* een_rescaled_e, const double* een_rescaled_single_n, const double* een_rescaled_single_e, double* const delta_een ); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_compute_jastrow_champ_factor_single_een (const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, const int64_t dim_c_vector, const double* c_vector_full, const int64_t* lkpm_combined_index, const double* tmp_c, const double* delta_p, const double* een_rescaled_n, const double* een_rescaled_e, const double* een_rescaled_single_n, const double* een_rescaled_single_e, double* const delta_een ) { #ifdef HAVE_HPC return qmckl_compute_jastrow_champ_factor_single_een_doc #else return qmckl_compute_jastrow_champ_factor_single_een_doc #endif (context, num, walk_num, elec_num, nucl_num, cord_num, dim_c_vector, c_vector_full, lkpm_combined_index, tmp_c, delta_p, een_rescaled_n, een_rescaled_e, een_rescaled_single_n, een_rescaled_single_e, delta_een ); } #+end_src #+CALL: generate_c_interface(table=qmckl_factor_single_een_args,rettyp=get_value("CRetType"),fname="qmckl_compute_jastrow_champ_factor_single_een_doc")) #+RESULTS: #+begin_src f90 :tangle (eval f) :comments org :exports none integer(c_int32_t) function qmckl_compute_jastrow_champ_factor_single_een_doc & (context, & num, & walk_num, & elec_num, & nucl_num, & cord_num, & dim_c_vector, & c_vector_full, & lkpm_combined_index, & tmp_c, & delta_p, & een_rescaled_n, & een_rescaled_e, & een_rescaled_single_n, & een_rescaled_single_e, & delta_een) & 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 :: num 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 integer (c_int64_t) , intent(in) , value :: dim_c_vector real (c_double ) , intent(in) :: c_vector_full(nucl_num,dim_c_vector) integer (c_int64_t) , intent(in) :: lkpm_combined_index(dim_c_vector,4) real (c_double ) , intent(in) :: tmp_c(elec_num,nucl_num,0:cord_num,0:cord_num-1,walk_num) real (c_double ) , intent(in) :: delta_p(elec_num,nucl_num,0:cord_num,0:cord_num-1,walk_num) real (c_double ) , intent(in) :: een_rescaled_n(elec_num,nucl_num,0:cord_num,walk_num) real (c_double ) , intent(in) :: een_rescaled_e(elec_num,elec_num,0:cord_num,walk_num) real (c_double ) , intent(in) :: een_rescaled_single_n(nucl_num,0:cord_num,walk_num) real (c_double ) , intent(in) :: een_rescaled_single_e(elec_num,0:cord_num,walk_num) real (c_double ) , intent(out) :: delta_een(walk_num) integer(c_int32_t), external :: qmckl_compute_jastrow_champ_factor_single_een_doc_f info = qmckl_compute_jastrow_champ_factor_single_een_doc_f & (context, & num, & walk_num, & elec_num, & nucl_num, & cord_num, & dim_c_vector, & c_vector_full, & lkpm_combined_index, & tmp_c, & delta_p, & een_rescaled_n, & een_rescaled_e, & een_rescaled_single_n, & een_rescaled_single_e, & delta_een) end function qmckl_compute_jastrow_champ_factor_single_een_doc #+end_src ** test #+begin_src c :tangle (eval c_test) /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); //rc = qmckl_check(context, // qmckl_set_electron_coord (context, 'N', walk_num, elec_coord, walk_num*3*elec_num) // ); qmckl_set_point(context, 'N', elec_num, elec_coord, walk_num*elec_num*3); assert(rc == QMCKL_SUCCESS); double jastrow_een_old[walk_num]; rc = qmckl_get_jastrow_champ_factor_een(context, &jastrow_een_old[0], walk_num); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_electron_coord(context, 'N', &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); rc = qmckl_set_single_point(context, 'N', 2, new_coords, 3); assert (rc == QMCKL_SUCCESS); double delta_een[walk_num]; rc = qmckl_get_jastrow_champ_single_een(context, &delta_een[0], walk_num); assert (rc == QMCKL_SUCCESS); coords[0][2][0] = new_coords[0]; coords[0][2][1] = new_coords[1]; coords[0][2][2] = new_coords[2]; rc = qmckl_set_point(context, 'N', elec_num, &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); //rc = qmckl_check(context, // qmckl_set_electron_coord (context, 'N', walk_num, &coords[0][0][0], walk_num*3*elec_num) // ); //assert(rc == QMCKL_SUCCESS); double jastrow_een_new[walk_num]; rc = qmckl_get_jastrow_champ_factor_een(context, &jastrow_een_new[0], walk_num); assert (rc == QMCKL_SUCCESS); for (int i = 0; i < walk_num; i++) { assert(fabs((jastrow_een_new[i]-jastrow_een_old[i])-delta_een[i]) < 1.e-12); } #+end_src * ee distance rescaled single point ** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_ee_rescaled_single(qmckl_context context, double* const distance_rescaled, const int64_t size_max); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_ee_rescaled_single(qmckl_context context, double* const distance_rescaled, const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_exit_code rc; rc = qmckl_provide_ee_rescaled_single(context); if (rc != QMCKL_SUCCESS) return rc; qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int64_t sze = ctx->electron.num * ctx->electron.walker.num; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "todo", "Array too small. Expected ctx->electron.num * ctx->electron.walker.num "); } memcpy(distance_rescaled, ctx->single_point.ee_rescaled_single, 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_rescaled_single(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_provide_ee_rescaled_single(qmckl_context context) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); /* Check if ee distance is provided */ qmckl_exit_code rc = qmckl_provide_single_ee_distance(context); if(rc != QMCKL_SUCCESS) return rc; /* Compute if necessary */ if (ctx->single_point.date > ctx->single_point.ee_rescaled_single_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { if (ctx->single_point.ee_rescaled_single!= NULL) { rc = qmckl_free(context, ctx->single_point.ee_rescaled_single); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_ee_rescaled_single", "Unable to free ctx->single_point.ee_rescaled_single"); } ctx->single_point.ee_rescaled_single = NULL; } } /* Allocate array */ if (ctx->single_point.ee_rescaled_single == NULL) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = ctx->electron.num * ctx->electron.walker.num * sizeof(double); double* ee_rescaled_single = (double*) qmckl_malloc(context, mem_info); if (ee_rescaled_single == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_ee_rescaled_single", NULL); } ctx->single_point.ee_rescaled_single = ee_rescaled_single; } rc = qmckl_compute_ee_rescaled_single(context, ctx->electron.num, ctx->jastrow_champ.rescale_factor_ee, ctx->electron.walker.num, ctx->single_point.single_ee_distance, ctx->single_point.ee_rescaled_single); if (rc != QMCKL_SUCCESS) { return rc; } ctx->single_point.ee_rescaled_single_date = ctx->single_point.date; } return QMCKL_SUCCESS; } #+end_src ** Compute :PROPERTIES: :Name: qmckl_compute_ee_rescaled_single :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_ee_rescaled_single_args | Variable | Type | In/Out | Description | |----------------------+------------------------------+--------+--------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~elec_num~ | ~int64_t~ | in | Number of electrons | | ~rescale_factor_ee~ | ~double~ | in | Factor to rescale ee distances | | ~walk_num~ | ~int64_t~ | in | Number of walkers | | ~single_ee_distance~ | ~double[walk_num][elec_num]~ | in | Electron coordinates | | ~ee_rescaled_single~ | ~double[walk_num][elec_num]~ | out | Electron-electron rescaled distances | #+begin_src f90 :comments org :tangle (eval f) :noweb yes function qmckl_compute_ee_rescaled_single_doc(context, & elec_num, rescale_factor_ee, walk_num, & single_ee_distance, ee_rescaled_single) & bind(C) result(info) use qmckl implicit none integer(qmckl_context), intent(in), value :: context integer (c_int64_t) , intent(in) , value :: elec_num real (c_double ) , intent(in) , value :: rescale_factor_ee integer (c_int64_t) , intent(in) , value :: walk_num real (c_double ) , intent(in) :: single_ee_distance(elec_num,walk_num) real (c_double ) , intent(out) :: ee_rescaled_single(elec_num,walk_num) integer(qmckl_exit_code) :: info integer*8 :: k, i real (c_double) :: inverse_rescale_factor_ee inverse_rescale_factor_ee = 1.0d0 / rescale_factor_ee 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 do i=1,elec_num ee_rescaled_single(i,k) = (1.0d0 - dexp(-rescale_factor_ee * single_ee_distance(i,k))) * inverse_rescale_factor_ee enddo end do end function qmckl_compute_ee_rescaled_single_doc #+end_src #+begin_src c :tangle (eval h_private_func) :comments org :exports none qmckl_exit_code qmckl_compute_ee_rescaled_single_doc ( const qmckl_context context, const int64_t elec_num, const double rescale_factor_ee, const int64_t walk_num, const double* single_ee_distance, double* const ee_rescaled_single ); qmckl_exit_code qmckl_compute_ee_rescaled_single( const qmckl_context context, const int64_t elec_num, const double rescale_factor_ee, const int64_t walk_num, const double* single_ee_distance, double* const ee_rescaled_single ); #+end_src #+begin_src c :tangle (eval c) :comments org :exports none qmckl_exit_code qmckl_compute_ee_rescaled_single ( const qmckl_context context, const int64_t elec_num, const double rescale_factor_ee, const int64_t walk_num, const double* single_ee_distance, double* const ee_rescaled_single ) { #ifdef HAVE_HPC return qmckl_compute_ee_rescaled_single_doc #else return qmckl_compute_ee_rescaled_single_doc #endif (context, elec_num, rescale_factor_ee, walk_num, single_ee_distance, ee_rescaled_single); } #+end_src ** test #+begin_src c :tangle (eval c_test) /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); rc = qmckl_set_point(context, 'N', elec_num, elec_coord, walk_num*elec_num*3); assert(rc == QMCKL_SUCCESS); double ee_rescaled[walk_num][elec_num][elec_num]; rc = qmckl_get_jastrow_champ_ee_distance_rescaled(context, &ee_rescaled[0][0][0]); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_electron_coord(context, 'N', &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); rc = qmckl_set_single_point(context, 'N', 2, new_coords, 3); assert (rc == QMCKL_SUCCESS); double single_ee_rescaled[walk_num][elec_num]; rc = qmckl_get_ee_rescaled_single(context, &single_ee_rescaled[0][0], walk_num*elec_num); assert (rc == QMCKL_SUCCESS); coords[0][2][0] = new_coords[0]; coords[0][2][1] = new_coords[1]; coords[0][2][2] = new_coords[2]; rc = qmckl_set_point(context, 'N', elec_num, &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_jastrow_champ_ee_distance_rescaled(context, &ee_rescaled[0][0][0]); assert (rc == QMCKL_SUCCESS); for (int nw = 0; nw < walk_num; nw++) { for (int i = 0; i < elec_num; i++){ if (i == 2) continue; assert(fabs(ee_rescaled[nw][2][i]-single_ee_rescaled[nw][i]) < 1.e-12); } } #+end_src * en distance rescaled single point ** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_en_rescaled_single(qmckl_context context, double* const distance_rescaled, const int64_t size_max); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_en_rescaled_single(qmckl_context context, double* const distance_rescaled, const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_exit_code rc; rc = qmckl_provide_en_rescaled_single(context); if (rc != QMCKL_SUCCESS) return rc; qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int64_t sze = ctx->nucleus.num * ctx->electron.walker.num; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "todo", "Array too small. Expected ctx->nucleus.num * ctx->electron.walker.num "); } memcpy(distance_rescaled, ctx->single_point.en_rescaled_single, 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_en_rescaled_single(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_provide_en_rescaled_single(qmckl_context context) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); /* Check if ee distance is provided */ qmckl_exit_code rc = qmckl_provide_single_en_distance(context); if(rc != QMCKL_SUCCESS) return rc; /* Compute if necessary */ if (ctx->single_point.date > ctx->single_point.en_rescaled_single_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { if (ctx->single_point.en_rescaled_single!= NULL) { rc = qmckl_free(context, ctx->single_point.en_rescaled_single); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_en_rescaled_single", "Unable to free ctx->single_point.en_rescaled_single"); } ctx->single_point.en_rescaled_single = NULL; } } /* Allocate array */ if (ctx->single_point.en_rescaled_single == NULL) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = ctx->nucleus.num * ctx->electron.walker.num * sizeof(double); double* en_rescaled_single = (double*) qmckl_malloc(context, mem_info); if (en_rescaled_single == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_en_rescaled_single", NULL); } ctx->single_point.en_rescaled_single = en_rescaled_single; } rc = qmckl_compute_en_rescaled_single(context, ctx->nucleus.num, ctx->jastrow_champ.type_nucl_num, ctx->jastrow_champ.type_nucl_vector, ctx->jastrow_champ.rescale_factor_en, ctx->electron.walker.num, ctx->single_point.single_en_distance, ctx->single_point.en_rescaled_single); if (rc != QMCKL_SUCCESS) { return rc; } ctx->single_point.en_rescaled_single_date = ctx->single_point.date; } return QMCKL_SUCCESS; } #+end_src ** Compute :PROPERTIES: :Name: qmckl_compute_en_rescaled_single :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_en_rescaled_single_args | Variable | Type | In/Out | Description | |----------------------+------------------------------+--------+--------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~elec_num~ | ~int64_t~ | in | Number of electrons | | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | | ~type_nucl_num~ | ~int64_t~ | in | Number of types of nuclei | | ~type_nucl_vector~ | ~int64_t[nucl_num]~ | in | Number of types of nuclei | | ~rescale_factor_en~ | ~double[type_nucl_num]~ | in | The factor for rescaled distances | | ~walk_num~ | ~int64_t~ | in | Number of walkers | | ~single_en_distance~ | ~double[walk_num][nucl_num]~ | in | Electron coordinates | | ~en_rescaled_single~ | ~double[walk_num][nucl_num]~ | out | Electron-electron rescaled distances | #+begin_src f90 :comments org :tangle (eval f) :noweb yes function qmckl_compute_en_rescaled_single_doc(context, & nucl_num, type_nucl_num, & type_nucl_vector, rescale_factor_en, walk_num, single_en_distance, en_rescaled_single) & bind(C) result(info) use qmckl implicit none integer (qmckl_context), intent(in), value :: context integer (c_int64_t) , intent(in) , value :: nucl_num integer (c_int64_t) , intent(in) , value :: type_nucl_num integer (c_int64_t) , intent(in) :: type_nucl_vector(nucl_num) real (c_double ) , intent(in) :: rescale_factor_en(type_nucl_num) integer (c_int64_t) , intent(in) , value :: walk_num real (c_double ) , intent(in) :: single_en_distance(nucl_num,walk_num) real (c_double ) , intent(out) :: en_rescaled_single(nucl_num,walk_num) integer(qmckl_exit_code) :: info integer*8 :: i, k double precision :: coord(3) info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif if (nucl_num <= 0) then info = QMCKL_INVALID_ARG_2 return endif if (walk_num <= 0) then info = QMCKL_INVALID_ARG_4 return endif do i=1, nucl_num do k=1,walk_num en_rescaled_single(i,k) = (1.0d0 - dexp(-rescale_factor_en(type_nucl_vector(i)+1) * & single_en_distance(i,k))) / rescale_factor_en(type_nucl_vector(i)+1) end do end do end function qmckl_compute_en_rescaled_single_doc #+end_src #+begin_src c :tangle (eval h_private_func) :comments org :exports none qmckl_exit_code qmckl_compute_en_rescaled_single_doc ( const qmckl_context context, const int64_t nucl_num, const int64_t type_nucl_num, const int64_t* type_nucl_vector, const double* rescale_factor_en, const int64_t walk_num, const double* single_en_distance, double* const en_rescaled_single ); qmckl_exit_code qmckl_compute_en_rescaled_single ( const qmckl_context context, const int64_t nucl_num, const int64_t type_nucl_num, const int64_t* type_nucl_vector, const double* rescale_factor_en, const int64_t walk_num, const double* single_en_distance, double* const en_rescaled_single ); #+end_src #+begin_src c :tangle (eval c) :comments org :exports none qmckl_exit_code qmckl_compute_en_rescaled_single( const qmckl_context context, const int64_t nucl_num, const int64_t type_nucl_num, const int64_t* type_nucl_vector, const double* rescale_factor_en, const int64_t walk_num, const double* single_en_distance, double* const en_rescaled_single ) { #ifdef HAVE_HPC return qmckl_compute_en_rescaled_single_doc #else return qmckl_compute_en_rescaled_single_doc #endif (context, nucl_num, type_nucl_num, type_nucl_vector, rescale_factor_en, walk_num, single_en_distance, en_rescaled_single ); } #+end_src ** test #+begin_src c :tangle (eval c_test) /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); rc = qmckl_set_point(context, 'N', elec_num, elec_coord, walk_num*elec_num*3); assert(rc == QMCKL_SUCCESS); double en_rescaled[walk_num][nucl_num][elec_num]; rc = qmckl_get_electron_en_distance_rescaled(context, &en_rescaled[0][0][0]); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_electron_coord(context, 'N', &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); rc = qmckl_set_single_point(context, 'N', 2, new_coords, 3); assert (rc == QMCKL_SUCCESS); double single_en_rescaled[walk_num][nucl_num]; rc = qmckl_get_en_rescaled_single(context, &single_en_rescaled[0][0], walk_num*nucl_num); assert (rc == QMCKL_SUCCESS); coords[0][2][0] = new_coords[0]; coords[0][2][1] = new_coords[1]; coords[0][2][2] = new_coords[2]; rc = qmckl_set_point(context, 'N', elec_num, &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_electron_en_distance_rescaled(context, &en_rescaled[0][0][0]); assert (rc == QMCKL_SUCCESS); for (int nw = 0; nw < walk_num; nw++) { for (int a = 0; a < nucl_num; a++){ assert(fabs(en_rescaled[nw][a][2]-single_en_rescaled[nw][a]) < 1.e-12); } } #+end_src * Delta ee ** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_jastrow_champ_single_ee(qmckl_context context, double* const delta_ee, const int64_t size_max); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_jastrow_champ_single_ee(qmckl_context context, double* const delta_ee, const int64_t size_max) { qmckl_exit_code rc; if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return qmckl_failwith( context, QMCKL_INVALID_CONTEXT, "qmckl_get_jastrow_champ_single_ee", NULL); } qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); rc = qmckl_provide_jastrow_champ_single_ee(context); if (rc != QMCKL_SUCCESS) return rc; int64_t sze=ctx->electron.walker.num; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_jastrow_champ_single_ee", "Array too small. Expected walker.num"); } memcpy(delta_ee, ctx->single_point.delta_ee, sze*sizeof(double)); return QMCKL_SUCCESS; } #+end_src #+begin_src f90 :tangle (eval fh_func) :comments org interface integer(qmckl_exit_code) function qmckl_get_jastrow_champ_single_ee (context, & delta_ee, size_max) bind(C) use, intrinsic :: iso_c_binding import implicit none integer (qmckl_context) , intent(in), value :: context integer(c_int64_t), intent(in), value :: size_max real(c_double), intent(out) :: delta_ee(size_max) end function qmckl_get_jastrow_champ_single_ee end interface #+end_src ** Provide :noexport: #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_provide_jastrow_champ_single_ee(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_provide_jastrow_champ_single_ee(qmckl_context context) { qmckl_exit_code rc; if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return qmckl_failwith( context, QMCKL_INVALID_CONTEXT, "qmckl_provide_jastrow_champ_single_ee", NULL); } qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); if (!ctx->jastrow_champ.provided) { return qmckl_failwith( context, QMCKL_NOT_PROVIDED, "qmckl_provide_jastrow_champ_single_ee", NULL); } rc = qmckl_provide_ee_distance_rescaled(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_ee_rescaled_single(context); if(rc != QMCKL_SUCCESS) return rc; /* Compute if necessary */ if (ctx->single_point.date > ctx->single_point.delta_ee_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { if (ctx->single_point.delta_ee != NULL) { rc = qmckl_free(context, ctx->single_point.delta_ee); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_jastrow_champ_single_ee", "Unable to free ctx->single_point.delta_ee"); } ctx->single_point.delta_ee = NULL; } } /* Allocate array */ if (ctx->single_point.delta_ee == NULL) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = ctx->electron.walker.num * sizeof(double); double* delta_ee = (double*) qmckl_malloc(context, mem_info); if (delta_ee == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_jastrow_champ_single_ee", NULL); } ctx->single_point.delta_ee = delta_ee; } rc = qmckl_compute_jastrow_champ_single_ee(context, ctx->single_point.num, ctx->electron.walker.num, ctx->electron.num, ctx->electron.up_num, ctx->jastrow_champ.bord_num, ctx->jastrow_champ.b_vector, ctx->jastrow_champ.ee_distance_rescaled, ctx->single_point.ee_rescaled_single, ctx->jastrow_champ.spin_independent, ctx->single_point.delta_ee); if (rc != QMCKL_SUCCESS) { return rc; } ctx->single_point.delta_ee_date = ctx->date; } return QMCKL_SUCCESS; } #+end_src ** Compute :PROPERTIES: :Name: qmckl_compute_jastrow_champ_single_ee :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_single_ee_args | Variable | Type | In/Out | Description | |------------------------+----------------------------------------+--------+-----------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~num~ | ~int64_t~ | in | Number of walkers | | ~walk_num~ | ~int64_t~ | in | Number of walkers | | ~elec_num~ | ~int64_t~ | in | Number of electrons | | ~up_num~ | ~int64_t~ | in | Number of alpha electrons | | ~bord_num~ | ~int64_t~ | in | Number of coefficients | | ~b_vector~ | ~double[bord_num+1]~ | in | List of coefficients | | ~ee_distance_rescaled~ | ~double[walk_num][elec_num][elec_num]~ | in | Electron-electron distances | | ~ee_rescaled_single~ | ~double[walk_num][][elec_num]~ | in | Electron-electron distances | | ~delta_ee~ | ~double[walk_num]~ | out | $f_{ee}$ | # #+CALL: generate_c_interface(table=qmckl_factor_ee_args,rettyp=get_value("CRetType"),fname=get_value("Name")) #+begin_src f90 :comments org :tangle (eval f) :noweb yes function qmckl_compute_jastrow_champ_single_ee_doc(context, & num_in, walk_num, elec_num, up_num, bord_num, b_vector, & ee_distance_rescaled, ee_rescaled_single, spin_independent, delta_ee) & bind(C) result(info) use qmckl implicit none integer (qmckl_context), intent(in), value :: context integer (c_int64_t) , intent(in), value :: num_in 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 :: up_num integer (c_int64_t) , intent(in), value :: bord_num real (c_double ) , intent(in) :: b_vector(bord_num+1) real (c_double ) , intent(in) :: ee_distance_rescaled(elec_num,elec_num,walk_num) real (c_double ) , intent(in) :: ee_rescaled_single(elec_num,walk_num) integer (c_int32_t) , intent(in), value :: spin_independent real (c_double ) , intent(out) :: delta_ee(walk_num) integer(qmckl_exit_code) :: info integer*8 :: i, j, k, nw, num double precision :: x, xk, y, yk num = num_in + 1 info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif if (walk_num <= 0) then info = QMCKL_INVALID_ARG_3 return endif if (elec_num <= 0) then info = QMCKL_INVALID_ARG_4 return endif if (bord_num < 0) then info = QMCKL_INVALID_ARG_5 return endif do nw =1, walk_num delta_ee(nw) = 0.0d0 do i=1,elec_num !print *,i, ee_rescaled_single(i,nw) !print *, i, ee_distance_rescaled(i,num,nw) !print *, ' ' if (i.ne.num) then x = ee_distance_rescaled(i,num,nw) y = ee_rescaled_single(i,nw) if (spin_independent == 1) then delta_ee(nw) = delta_ee(nw) - (b_vector(1) * x / (1.d0 + b_vector(2) * x)) & + (b_vector(1) * y / (1.d0 + b_vector(2) * y)) else if ( (i <= up_num).or.(num > up_num) ) then delta_ee(nw) = delta_ee(nw) - (0.5d0 * b_vector(1) * x / (1.d0 + b_vector(2) * x)) & + (0.5d0 * b_vector(1) * y / (1.d0 + b_vector(2) * y)) else delta_ee(nw) = delta_ee(nw) - (b_vector(1) * x / (1.d0 + b_vector(2) * x)) & + (b_vector(1) * y / (1.d0 + b_vector(2) * y)) endif endif xk = x yk = y do k=2,bord_num xk = xk * x yk = yk * y delta_ee(nw) = delta_ee(nw) - (b_vector(k+1) * xk) + (b_vector(k+1) * yk) end do endif end do end do end function qmckl_compute_jastrow_champ_single_ee_doc #+end_src #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_jastrow_champ_single_ee_doc (const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t up_num, const int64_t bord_num, const double* b_vector, const double* ee_distance_rescaled, const double* ee_rescaled_singe, const int32_t spin_independent, double* const delta_ee ); #+end_src # #+CALL: generate_c_header(table=qmckl_factor_ee_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_jastrow_champ_single_ee (const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t up_num, const int64_t bord_num, const double* b_vector, const double* ee_distance_rescaled, const double* ee_rescaled_single, const int32_t spin_independent, double* const delta_ee ); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes qmckl_exit_code qmckl_compute_jastrow_champ_single_ee (const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t up_num, const int64_t bord_num, const double* b_vector, const double* ee_distance_rescaled, const double* ee_rescaled_single, const int32_t spin_independent, double* const delta_ee ) { #ifdef HAVE_HPC return qmckl_compute_jastrow_champ_single_ee_doc #else return qmckl_compute_jastrow_champ_single_ee_doc #endif (context, num, walk_num, elec_num, up_num, bord_num, b_vector, ee_distance_rescaled, ee_rescaled_single, spin_independent, delta_ee); } #+end_src ** test #+begin_src c :tangle (eval c_test) /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); rc = qmckl_set_point(context, 'N', elec_num, elec_coord, walk_num*elec_num*3); assert(rc == QMCKL_SUCCESS); double jastrow_ee_old[walk_num]; rc = qmckl_get_jastrow_champ_factor_ee(context, &jastrow_ee_old[0], walk_num); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_electron_coord(context, 'N', &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); rc = qmckl_set_single_point(context, 'N', 2, new_coords, 3); assert (rc == QMCKL_SUCCESS); double delta_ee[walk_num]; rc = qmckl_get_jastrow_champ_single_ee(context, &delta_ee[0], walk_num); assert (rc == QMCKL_SUCCESS); coords[0][2][0] = new_coords[0]; coords[0][2][1] = new_coords[1]; coords[0][2][2] = new_coords[2]; rc = qmckl_set_point(context, 'N', elec_num, &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); double jastrow_ee_new[walk_num]; rc = qmckl_get_jastrow_champ_factor_ee(context, &jastrow_ee_new[0], walk_num); assert (rc == QMCKL_SUCCESS); for (int nw = 0; nw < walk_num; nw++) { assert(fabs((jastrow_ee_new[nw] - jastrow_ee_old[nw]) - delta_ee[nw]) < 1.e-12); } #+end_src * Delta en ** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_jastrow_champ_single_en(qmckl_context context, double* const delta_en, const int64_t size_max); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_jastrow_champ_single_en(qmckl_context context, double* const delta_en, const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return qmckl_failwith( context, QMCKL_INVALID_CONTEXT, "qmckl_get_jastrow_champ_single_en", NULL); } qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); qmckl_exit_code rc; rc = qmckl_provide_jastrow_champ_single_en(context); if (rc != QMCKL_SUCCESS) return rc; int64_t sze=ctx->electron.walker.num; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_jastrow_champ_single_en", "Array too small. Expected walker.num"); } memcpy(delta_en, ctx->single_point.delta_en, sze*sizeof(double)); return QMCKL_SUCCESS; } #+end_src #+begin_src f90 :tangle (eval fh_func) :comments org interface integer(qmckl_exit_code) function qmckl_get_jastrow_champ_single_en (context, & delta_en, size_max) bind(C) use, intrinsic :: iso_c_binding import implicit none integer (qmckl_context) , intent(in), value :: context integer(c_int64_t), intent(in), value :: size_max real(c_double), intent(out) :: delta_en(size_max) end function qmckl_get_jastrow_champ_single_en end interface #+end_src ** Provide :noexport: #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_provide_jastrow_champ_single_en(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_provide_jastrow_champ_single_en(qmckl_context context) { qmckl_exit_code rc; if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return qmckl_failwith( context, QMCKL_INVALID_CONTEXT, "qmckl_provide_jastrow_champ_single_en", NULL); } qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); if (!ctx->jastrow_champ.provided) { return qmckl_failwith( context, QMCKL_NOT_PROVIDED, "qmckl_provide_jastrow_champ_single_en", NULL); } /* Check if en rescaled distance is provided */ rc = qmckl_provide_en_distance_rescaled(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_en_rescaled_single(context); if(rc != QMCKL_SUCCESS) return rc; /* Compute if necessary */ if (ctx->single_point.date > ctx->single_point.delta_en_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { if (ctx->single_point.delta_en != NULL) { rc = qmckl_free(context, ctx->single_point.delta_en); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_jastrow_champ_single_en", "Unable to free ctx->single_point.delta_en"); } ctx->single_point.delta_en = NULL; } } /* Allocate array */ if (ctx->single_point.delta_en == NULL) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = ctx->electron.walker.num * sizeof(double); double* delta_en = (double*) qmckl_malloc(context, mem_info); if (delta_en == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_jastrow_champ_single_en", NULL); } ctx->single_point.delta_en = delta_en; } rc = qmckl_compute_jastrow_champ_single_en(context, ctx->single_point.num, ctx->electron.walker.num, ctx->electron.num, ctx->nucleus.num, ctx->jastrow_champ.type_nucl_num, ctx->jastrow_champ.type_nucl_vector, ctx->jastrow_champ.aord_num, ctx->jastrow_champ.a_vector, ctx->jastrow_champ.en_distance_rescaled, ctx->single_point.en_rescaled_single, ctx->single_point.delta_en); if (rc != QMCKL_SUCCESS) { return rc; } ctx->single_point.delta_en_date = ctx->single_point.date; } return QMCKL_SUCCESS; } #+end_src ** Compute :PROPERTIES: :Name: qmckl_compute_jastrow_champ_single_en_doc :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_single_en_args | Variable | Type | In/Out | Description | |------------------------+----------------------------------------+--------+----------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~num~ | ~int64_t~ | in | Number of single point | | ~walk_num~ | ~int64_t~ | in | Number of walkers | | ~elec_num~ | ~int64_t~ | in | Number of electrons | | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | | ~type_nucl_num~ | ~int64_t~ | in | Number of unique nuclei | | ~type_nucl_vector~ | ~int64_t[nucl_num]~ | in | IDs of unique nuclei | | ~aord_num~ | ~int64_t~ | in | Number of coefficients | | ~a_vector~ | ~double[type_nucl_num][aord_num+1]~ | in | List of coefficients | | ~en_distance_rescaled~ | ~double[walk_num][nucl_num][elec_num]~ | in | Electron-nucleus distances | | ~en_rescaled_single ~ | ~double[walk_num][nucl_num]~ | in | Electron-nucleus distances | | ~delta_en~ | ~double[walk_num]~ | out | Electron-nucleus jastrow | #+begin_src f90 :comments org :tangle (eval f) :noweb yes function qmckl_compute_jastrow_champ_single_en_doc( & context, num_in, walk_num, elec_num, nucl_num, type_nucl_num, & type_nucl_vector, aord_num, a_vector, & en_distance_rescaled, en_rescaled_single, delta_en) & bind(C) result(info) use qmckl implicit none integer (qmckl_context), intent(in), value :: context integer (c_int64_t) , intent(in) , value :: num_in 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 :: type_nucl_num integer (c_int64_t) , intent(in) :: type_nucl_vector(nucl_num) integer (c_int64_t) , intent(in) , value :: aord_num real (c_double ) , intent(in) :: a_vector(aord_num+1,type_nucl_num) real (c_double ) , intent(in) :: en_distance_rescaled(elec_num,nucl_num,walk_num) real (c_double ) , intent(in) :: en_rescaled_single(nucl_num,walk_num) real (c_double ) , intent(out) :: delta_en(walk_num) integer(qmckl_exit_code) :: info integer*8 :: i, a, p, nw, num double precision :: x, power_ser, y num = num_in + 1 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 (type_nucl_num <= 0) then info = QMCKL_INVALID_ARG_4 return endif if (aord_num < 0) then info = QMCKL_INVALID_ARG_7 return endif do nw =1, walk_num delta_en(nw) = 0.0d0 do a = 1, nucl_num x = en_distance_rescaled(num, a, nw) y = en_rescaled_single(a, nw) delta_en(nw) = delta_en(nw) - a_vector(1, type_nucl_vector(a)+1) * x / (1.0d0 + a_vector(2, type_nucl_vector(a)+1) * x) delta_en(nw) = delta_en(nw) + a_vector(1, type_nucl_vector(a)+1) * y / (1.0d0 + a_vector(2, type_nucl_vector(a)+1) * y) do p = 2, aord_num x = x * en_distance_rescaled(num, a, nw) y = y * en_rescaled_single(a, nw) delta_en(nw) = delta_en(nw) - a_vector(p + 1, type_nucl_vector(a)+1) * x + a_vector(p + 1, type_nucl_vector(a)+1) * y end do end do end do end function qmckl_compute_jastrow_champ_single_en_doc #+end_src #+CALL: generate_c_header(table=qmckl_single_en_args,rettyp=get_value("CRetType"),fname=get_value("Name")) #+RESULTS: #+begin_src c :tangle (eval h_func) :comments org qmckl_exit_code qmckl_compute_jastrow_champ_single_en ( const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t type_nucl_num, const int64_t* type_nucl_vector, const int64_t aord_num, const double* a_vector, const double* en_distance_rescaled, const double* en_rescaled_single, double* const delta_en ); qmckl_exit_code qmckl_compute_jastrow_champ_single_en_doc ( const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t type_nucl_num, const int64_t* type_nucl_vector, const int64_t aord_num, const double* a_vector, const double* en_distance_rescaled, const double* en_rescaled_single, double* const delta_en ); #+end_src #+begin_src c :tangle (eval c) :comments org qmckl_exit_code qmckl_compute_jastrow_champ_single_en ( const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t type_nucl_num, const int64_t* type_nucl_vector, const int64_t aord_num, const double* a_vector, const double* en_distance_rescaled, const double* en_rescaled_single, double* const delta_en ) { #ifdef HAVE_HPC return qmckl_compute_jastrow_champ_single_en_doc #else return qmckl_compute_jastrow_champ_single_en_doc #endif (context, num, walk_num, elec_num, nucl_num, type_nucl_num, type_nucl_vector, aord_num, a_vector, en_distance_rescaled, en_rescaled_single, delta_en ); } #+end_src ** test #+begin_src c :tangle (eval c_test) /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); rc = qmckl_set_point(context, 'N', elec_num, elec_coord, walk_num*elec_num*3); assert(rc == QMCKL_SUCCESS); double jastrow_en_old[walk_num]; rc = qmckl_get_jastrow_champ_factor_en(context, &jastrow_en_old[0], walk_num); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_electron_coord(context, 'N', &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); rc = qmckl_set_single_point(context, 'N', 2, new_coords, 3); assert (rc == QMCKL_SUCCESS); double delta_en[walk_num]; rc = qmckl_get_jastrow_champ_single_en(context, &delta_en[0], walk_num); assert (rc == QMCKL_SUCCESS); coords[0][2][0] = new_coords[0]; coords[0][2][1] = new_coords[1]; coords[0][2][2] = new_coords[2]; rc = qmckl_set_point(context, 'N', elec_num, &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); double jastrow_en_new[walk_num]; rc = qmckl_get_jastrow_champ_factor_en(context, &jastrow_en_new[0], walk_num); assert (rc == QMCKL_SUCCESS); for (int nw = 0; nw < walk_num; nw++) { assert(fabs((jastrow_en_new[nw] - jastrow_en_old[nw]) - delta_en[nw]) < 1.e-12); } #+end_src * En rescaled derivative een ** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_een_rescaled_single_n_gl(qmckl_context context, double* const distance_rescaled, const int64_t size_max); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_een_rescaled_single_n_gl(qmckl_context context, double* const distance_rescaled, const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_exit_code rc; rc = qmckl_provide_een_rescaled_single_n_gl(context); if (rc != QMCKL_SUCCESS) return rc; qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int64_t sze = 4 * ctx->nucleus.num * ctx->electron.walker.num * (ctx->jastrow_champ.cord_num + 1); if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_jastrow_champ_single_een_gl", "Array too small. Expected 4 * ctx->nucleus.num * ctx->electron.walker.num * (ctx->jastrow_champ.cord_num + 1)"); } memcpy(distance_rescaled, ctx->single_point.een_rescaled_single_n_gl, 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_single_n_gl(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_provide_een_rescaled_single_n_gl(qmckl_context context) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); /* Check if ee distance is provided */ qmckl_exit_code rc = qmckl_provide_single_en_distance(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if ee distance is provided */ rc = qmckl_provide_een_rescaled_single_n(context); if(rc != QMCKL_SUCCESS) return rc; /* Compute if necessary */ if (ctx->single_point.date > ctx->single_point.een_rescaled_single_n_gl_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { if (ctx->single_point.een_rescaled_single_n_gl != NULL) { rc = qmckl_free(context, ctx->single_point.een_rescaled_single_n_gl); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_een_rescaled_single_n_gl", "Unable to free ctx->single_pont.een_rescaled_single_n_gl"); } ctx->single_point.een_rescaled_single_n_gl = NULL; } } /* Allocate array */ if (ctx->single_point.een_rescaled_single_n_gl == NULL) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = 4 * ctx->nucleus.num * ctx->electron.walker.num * (ctx->jastrow_champ.cord_num + 1) * sizeof(double); double* een_rescaled_single_n_gl = (double*) qmckl_malloc(context, mem_info); if (een_rescaled_single_n_gl == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_een_rescaled_single_n_gl", NULL); } ctx->single_point.een_rescaled_single_n_gl = een_rescaled_single_n_gl; } rc = qmckl_compute_een_rescaled_single_n_gl(context, ctx->electron.walker.num, ctx->nucleus.num, ctx->jastrow_champ.type_nucl_num, ctx->jastrow_champ.type_nucl_vector, ctx->jastrow_champ.cord_num, ctx->jastrow_champ.rescale_factor_en, ctx->single_point.coord.data, ctx->nucleus.coord.data, ctx->single_point.single_en_distance, ctx->single_point.een_rescaled_single_n, ctx->single_point.een_rescaled_single_n_gl); if (rc != QMCKL_SUCCESS) { return rc; } ctx->single_point.een_rescaled_single_n_gl_date = ctx->single_point.date; } return QMCKL_SUCCESS; } #+end_src ** Compute :PROPERTIES: :Name: qmckl_compute_een_rescaled_single_n_gl :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_compute_een_rescaled_single_n_gl_args | Variable | Type | In/Out | Description | |----------------------------+---------------------------------------------+--------+-------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~walk_num~ | ~int64_t~ | in | Number of walkers | | ~nucl_num~ | ~int64_t~ | in | Number of atoms | | ~type_nucl_num~ | ~int64_t~ | in | Number of atom types | | ~type_nucl_vector~ | ~int64_t[nucl_num]~ | in | Types of atoms | | ~cord_num~ | ~int64_t~ | in | Order of polynomials | | ~rescale_factor_en~ | ~double[nucl_num]~ | in | Factor to rescale ee distances | | ~coord_ee~ | ~double[walk_num][3]~ | in | Electron coordinates | | ~coord_n~ | ~double[3][nucl_num]~ | in | Nuclear coordinates | | ~single_en_distance~ | ~double[walk_num][nucl_num]~ | in | Electron-nucleus distances | | ~een_rescaled_single_n~ | ~double[walk_num][0:cord_num][nucl_num]~ | in | Electron-nucleus distances | | ~een_rescaled_single_n_gl~ | ~double[walk_num][0:cord_num][nucl_num][4]~ | out | Electron-nucleus rescaled distances | #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_een_rescaled_single_n_gl_f( & context, walk_num, nucl_num, type_nucl_num, type_nucl_vector, & cord_num, rescale_factor_en, & coord_ee, coord_n, single_en_distance, een_rescaled_single_n, een_rescaled_single_n_gl) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: walk_num integer*8 , intent(in) :: nucl_num integer*8 , intent(in) :: type_nucl_num 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) :: coord_ee(3) double precision , intent(in) :: coord_n(nucl_num,3) double precision , intent(in) :: single_en_distance(nucl_num,walk_num) double precision , intent(in) :: een_rescaled_single_n(nucl_num,0:cord_num,walk_num) double precision , intent(out) :: een_rescaled_single_n_gl(4,nucl_num,0:cord_num,walk_num) double precision,dimension(:,:),allocatable :: elnuc_dist_gl double precision :: x, ria_inv, kappa_l integer*8 :: i, a, k, l, nw, ii allocate(elnuc_dist_gl(4, nucl_num)) 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 (nucl_num <= 0) then info = QMCKL_INVALID_ARG_3 return endif if (cord_num < 0) then info = QMCKL_INVALID_ARG_4 return endif ! Prepare table of exponentiated distances raised to appropriate power een_rescaled_single_n_gl = 0.0d0 do nw = 1, walk_num ! prepare the actual een table do a = 1, nucl_num ria_inv = 1.0d0 / single_en_distance(a, nw) do ii = 1, 3 elnuc_dist_gl(ii, a) = (coord_ee(ii) - coord_n(a, ii)) * ria_inv end do elnuc_dist_gl(4, a) = 2.0d0 * ria_inv end do do l = 0, cord_num do a = 1, nucl_num kappa_l = - dble(l) * rescale_factor_en(type_nucl_vector(a)+1) een_rescaled_single_n_gl(1, a, l, nw) = kappa_l * elnuc_dist_gl(1, a) een_rescaled_single_n_gl(2, a, l, nw) = kappa_l * elnuc_dist_gl(2, a) een_rescaled_single_n_gl(3, a, l, nw) = kappa_l * elnuc_dist_gl(3, a) een_rescaled_single_n_gl(4, a, l, nw) = kappa_l * elnuc_dist_gl(4, a) een_rescaled_single_n_gl(4, a, l, nw) = een_rescaled_single_n_gl(4, a, l, nw) & + een_rescaled_single_n_gl(1, a, l, nw) * een_rescaled_single_n_gl(1, a, l, nw) & + een_rescaled_single_n_gl(2, a, l, nw) * een_rescaled_single_n_gl(2, a, l, nw) & + een_rescaled_single_n_gl(3, a, l, nw) * een_rescaled_single_n_gl(3, a, l, nw) een_rescaled_single_n_gl(1, a, l, nw) = een_rescaled_single_n_gl(1, a, l, nw) * & een_rescaled_single_n(a, l, nw) een_rescaled_single_n_gl(2, a, l, nw) = een_rescaled_single_n_gl(2, a, l, nw) * & een_rescaled_single_n(a, l, nw) een_rescaled_single_n_gl(3, a, l, nw) = een_rescaled_single_n_gl(3, a, l, nw) * & een_rescaled_single_n(a, l, nw) een_rescaled_single_n_gl(4, a, l, nw) = een_rescaled_single_n_gl(4, a, l, nw) * & een_rescaled_single_n(a, l, nw) end do end do end do end function qmckl_compute_een_rescaled_single_n_gl_f #+end_src # #+CALL: generate_c_header(table=qmckl_compute_jastrow_champ_factor_een_rescaled_n_gl_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_een_rescaled_single_n_gl ( const qmckl_context context, const int64_t walk_num, const int64_t nucl_num, const int64_t type_nucl_num, int64_t* const type_nucl_vector, const int64_t cord_num, const double* rescale_factor_en, const double* coord_ee, const double* coord_n, const double* single_en_distance, const double* een_rescaled_single_n, double* const een_rescaled_single_n_gl ); #+end_src #+CALL: generate_c_interface(table=qmckl_compute_een_rescaled_single_n_gl_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_single_n_gl & (context, & walk_num, & nucl_num, & type_nucl_num, & type_nucl_vector, & cord_num, & rescale_factor_en, & coord_ee, & coord_n, & single_en_distance, & een_rescaled_single_n, & een_rescaled_single_n_gl) & 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 :: nucl_num integer (c_int64_t) , intent(in) , value :: type_nucl_num integer (c_int64_t) , intent(in) :: type_nucl_vector(nucl_num) integer (c_int64_t) , intent(in) , value :: cord_num real (c_double ) , intent(in) :: rescale_factor_en(nucl_num) real (c_double ) , intent(in) :: coord_ee(3) real (c_double ) , intent(in) :: coord_n(nucl_num,3) real (c_double ) , intent(in) :: single_en_distance(nucl_num,walk_num) real (c_double ) , intent(in) :: een_rescaled_single_n(nucl_num,0:cord_num,walk_num) real (c_double ) , intent(out) :: een_rescaled_single_n_gl(4,nucl_num,0:cord_num,walk_num) integer(c_int32_t), external :: qmckl_compute_een_rescaled_single_n_gl_f info = qmckl_compute_een_rescaled_single_n_gl_f & (context, & walk_num, & nucl_num, & type_nucl_num, & type_nucl_vector, & cord_num, & rescale_factor_en, & coord_ee, & coord_n, & single_en_distance, & een_rescaled_single_n, & een_rescaled_single_n_gl) end function qmckl_compute_een_rescaled_single_n_gl #+end_src ** test #+begin_src c :tangle (eval c_test) /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); qmckl_set_point(context, 'N', elec_num, elec_coord, walk_num*elec_num*3); assert(rc == QMCKL_SUCCESS); double een_rescaled_en_gl[walk_num][cord_num+1][nucl_num][4][elec_num]; rc = qmckl_get_jastrow_champ_een_rescaled_n_gl(context, &een_rescaled_en_gl[0][0][0][0][0], walk_num*(cord_num+1)*nucl_num*elec_num*4); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_electron_coord(context, 'N', &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); rc = qmckl_set_single_point(context, 'N', 2, new_coords, 3); assert (rc == QMCKL_SUCCESS); double een_rescaled_single_n_gl[walk_num][cord_num+1][nucl_num][4]; rc = qmckl_get_een_rescaled_single_n_gl(context, &een_rescaled_single_n_gl[0][0][0][0], walk_num*(cord_num+1)*nucl_num*4); assert (rc == QMCKL_SUCCESS); coords[0][2][0] = new_coords[0]; coords[0][2][1] = new_coords[1]; coords[0][2][2] = new_coords[2]; rc = qmckl_set_point(context, 'N', elec_num, &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_jastrow_champ_een_rescaled_n_gl(context, &een_rescaled_en_gl[0][0][0][0][0], walk_num*(cord_num+1)*nucl_num*elec_num*4); assert (rc == QMCKL_SUCCESS); for (int l = 0; l < cord_num+1; l++) { for (int nw = 0; nw < walk_num; nw++) { for (int a = 0; a < nucl_num; a++) { for (int m = 0; m < 4; m++) { assert(fabs(een_rescaled_en_gl[nw][l][a][m][2] - een_rescaled_single_n_gl[nw][l][a][m]) < 1.e-12); } } } } #+end_src * EE rescaled distance derivative for een ** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_een_rescaled_single_e_gl(qmckl_context context, double* const distance_rescaled, const int64_t size_max); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_een_rescaled_single_e_gl(qmckl_context context, double* const distance_rescaled, const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_exit_code rc; rc = qmckl_provide_een_rescaled_single_e_gl(context); if (rc != QMCKL_SUCCESS) return rc; qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int64_t sze = 4 * ctx->electron.num * ctx->electron.walker.num * (ctx->jastrow_champ.cord_num + 1); if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_jastrow_champ_factor_een_gl", "Array too small. Expected 4 * ctx->electron.num * ctx->electron.walker.num * (ctx->jastrow_champ.cord_num + 1)"); } memcpy(distance_rescaled, ctx->single_point.een_rescaled_single_e_gl, 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_single_e_gl(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_provide_een_rescaled_single_e_gl(qmckl_context context) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); /* Check if rescaled een-ee distance is provided */ qmckl_exit_code rc = qmckl_provide_een_rescaled_single_e(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_single_en_distance(context); if(rc != QMCKL_SUCCESS) return rc; /* Compute if necessary */ if (ctx->single_point.date > ctx->single_point.een_rescaled_single_e_gl_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { if (ctx->single_point.een_rescaled_single_e_gl != NULL) { rc = qmckl_free(context, ctx->single_point.een_rescaled_single_e_gl); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_een_rescaled_e_gl", "Unable to free ctx->single_point.een_rescaled_single_e_gl"); } ctx->single_point.een_rescaled_single_e_gl = NULL; } } /* Allocate array */ if (ctx->single_point.een_rescaled_single_e_gl == NULL) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = 4 * ctx->electron.num * ctx->electron.walker.num * (ctx->jastrow_champ.cord_num + 1) * sizeof(double); double* een_rescaled_single_e_gl = (double*) qmckl_malloc(context, mem_info); if (een_rescaled_single_e_gl == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_een_rescaled_single_e_gl", NULL); } ctx->single_point.een_rescaled_single_e_gl = een_rescaled_single_e_gl; } rc = qmckl_compute_een_rescaled_single_e_gl(context, ctx->single_point.num, ctx->electron.walker.num, ctx->electron.num, ctx->jastrow_champ.cord_num, ctx->jastrow_champ.rescale_factor_ee, ctx->single_point.coord.data, ctx->electron.walker.point.coord.data, ctx->single_point.single_ee_distance, ctx->single_point.een_rescaled_single_e, ctx->single_point.een_rescaled_single_e_gl); if (rc != QMCKL_SUCCESS) { return rc; } ctx->single_point.een_rescaled_single_e_gl_date = ctx->single_point.date; } return QMCKL_SUCCESS; } #+end_src ** Compute :PROPERTIES: :Name: qmckl_compute_een_rescaled_single_e_gl :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_een_rescaled_single_e_gl_args | Variable | Type | In/Out | Description | |----------------------------+---------------------------------------------+--------+--------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~num~ | ~int64_t~ | in | Number of walkers | | ~walk_num~ | ~int64_t~ | in | Number of walkers | | ~elec_num~ | ~int64_t~ | in | Number of electrons | | ~cord_num~ | ~int64_t~ | in | Order of polynomials | | ~rescale_factor_ee~ | ~double~ | in | Factor to rescale ee distances | | ~coord~ | ~double[3]~ | in | Electron coordinates | | ~coord_ee~ | ~double[walk_num][3][elec_num]~ | in | Electron coordinates | | ~single_ee_distance~ | ~double[walk_num][elec_num]~ | in | Electron-electron distances | | ~een_rescaled_single_e~ | ~double[walk_num][0:cord_num][elec_num]~ | in | Electron-electron distances | | ~een_rescaled_single_e_gl~ | ~double[walk_num][0:cord_num][elec_num][4]~ | out | Electron-electron rescaled distances | #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_een_rescaled_single_e_gl_f( & context, num_in, walk_num, elec_num, cord_num, rescale_factor_ee, & coord, coord_ee, single_ee_distance, een_rescaled_single_e, een_rescaled_single_e_gl) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: num_in integer*8 , intent(in) :: walk_num integer*8 , intent(in) :: elec_num integer*8 , intent(in) :: cord_num double precision , intent(in) :: rescale_factor_ee double precision , intent(in) :: coord(3) double precision , intent(in) :: coord_ee(elec_num,3,walk_num) double precision , intent(in) :: single_ee_distance(elec_num,walk_num) double precision , intent(in) :: een_rescaled_single_e(elec_num,0:cord_num,walk_num) double precision , intent(out) :: een_rescaled_single_e_gl(4,elec_num,0:cord_num,walk_num) double precision,dimension(:,:),allocatable :: elec_dist_gl double precision :: x, rij_inv, kappa_l integer*8 :: i, j, k, l, nw, ii, num num = num_in + 1 allocate(elec_dist_gl(4, elec_num)) info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif if (walk_num <= 0) then info = QMCKL_INVALID_ARG_3 return endif if (elec_num <= 0) then info = QMCKL_INVALID_ARG_4 return endif if (cord_num < 0) then info = QMCKL_INVALID_ARG_5 return endif ! Not necessary: should be set to zero by qmckl_malloc ! een_rescaled_single_e_gl = 0.0d0 ! Prepare table of exponentiated distances raised to appropriate power do nw = 1, walk_num do i = 1, elec_num rij_inv = 1.0d0 / single_ee_distance(i, nw) do ii = 1, 3 elec_dist_gl(ii, i) = (coord(ii) - coord_ee(i, ii, nw)) * rij_inv end do elec_dist_gl(4, i) = 2.0d0 * rij_inv end do elec_dist_gl(:, num) = 0.0d0 do l = 1, cord_num kappa_l = - dble(l) * rescale_factor_ee do i = 1, elec_num een_rescaled_single_e_gl(1, i, l, nw) = kappa_l * elec_dist_gl(1, i) een_rescaled_single_e_gl(2, i, l, nw) = kappa_l * elec_dist_gl(2, i) een_rescaled_single_e_gl(3, i, l, nw) = kappa_l * elec_dist_gl(3, i) een_rescaled_single_e_gl(4, i, l, nw) = kappa_l * elec_dist_gl(4, i) een_rescaled_single_e_gl(4, i, l, nw) = een_rescaled_single_e_gl(4, i, l, nw) & + een_rescaled_single_e_gl(1, i, l, nw) * een_rescaled_single_e_gl(1, i, l, nw) & + een_rescaled_single_e_gl(2, i, l, nw) * een_rescaled_single_e_gl(2, i, l, nw) & + een_rescaled_single_e_gl(3, i, l, nw) * een_rescaled_single_e_gl(3, i, l, nw) een_rescaled_single_e_gl(1,i,l,nw) = een_rescaled_single_e_gl(1,i,l,nw) * een_rescaled_single_e(i,l,nw) een_rescaled_single_e_gl(2,i,l,nw) = een_rescaled_single_e_gl(2,i,l,nw) * een_rescaled_single_e(i,l,nw) een_rescaled_single_e_gl(3,i,l,nw) = een_rescaled_single_e_gl(3,i,l,nw) * een_rescaled_single_e(i,l,nw) een_rescaled_single_e_gl(4,i,l,nw) = een_rescaled_single_e_gl(4,i,l,nw) * een_rescaled_single_e(i,l,nw) end do end do end do end function qmckl_compute_een_rescaled_single_e_gl_f #+end_src # #+CALL: generate_c_header(table=qmckl_factor_een_rescaled_e_gl_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_een_rescaled_single_e_gl ( const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t cord_num, const double rescale_factor_ee, const double* coord, const double* coord_ee, const double* single_ee_distance, const double* een_rescaled_single_e, double* const een_rescaled_single_e_gl ); #+end_src #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_een_rescaled_single_e_gl_doc ( const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t cord_num, const double rescale_factor_ee, const double* coord, const double* coord_ee, const double* single_ee_distance, const double* een_rescaled_single_e, double* const een_rescaled_single_e_gl ); #+end_src #+CALL: generate_c_interface(table=qmckl_een_rescaled_single_e_gl_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_single_e_gl_doc & (context, & num, & walk_num, & elec_num, & cord_num, & rescale_factor_ee, & coord, & coord_ee, & single_ee_distance, & een_rescaled_single_e, & een_rescaled_single_e_gl) & 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 :: num 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 :: cord_num real (c_double ) , intent(in) , value :: rescale_factor_ee real (c_double ) , intent(in) :: coord(3) real (c_double ) , intent(in) :: coord_ee(elec_num, 3, walk_num) real (c_double ) , intent(in) :: single_ee_distance(elec_num,walk_num) real (c_double ) , intent(in) :: een_rescaled_single_e(elec_num,0:cord_num,walk_num) real (c_double ) , intent(out) :: een_rescaled_single_e_gl(4,elec_num,0:cord_num,walk_num) integer(c_int32_t), external :: qmckl_compute_een_rescaled_single_e_gl_f info = qmckl_compute_een_rescaled_single_e_gl_f & (context, & num, & walk_num, & elec_num, & cord_num, & rescale_factor_ee, & coord, & coord_ee, & single_ee_distance, & een_rescaled_single_e, & een_rescaled_single_e_gl) end function qmckl_compute_een_rescaled_single_e_gl_doc #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes qmckl_exit_code qmckl_compute_een_rescaled_single_e_gl ( const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t cord_num, const double rescale_factor_ee, const double* coord, const double* coord_ee, const double* single_ee_distance, const double* een_rescaled_single_e, double* const een_rescaled_single_e_gl ) { #ifdef HAVE_HPC return qmckl_compute_een_rescaled_single_e_gl_doc #else return qmckl_compute_een_rescaled_single_e_gl_doc #endif (context, num, walk_num, elec_num, cord_num, rescale_factor_ee, coord, coord_ee, single_ee_distance, een_rescaled_single_e, een_rescaled_single_e_gl ); } #+end_src ** test #+begin_src c :tangle (eval c_test) /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); //rc = qmckl_check(context, // qmckl_set_electron_coord (context, 'N', walk_num, elec_coord, walk_num*3*elec_num) // ); qmckl_set_point(context, 'N', elec_num, elec_coord, walk_num*elec_num*3); assert(rc == QMCKL_SUCCESS); double een_rescaled_ee_gl[walk_num][cord_num+1][elec_num][4][elec_num]; rc = qmckl_get_jastrow_champ_een_rescaled_e_gl(context, &een_rescaled_ee_gl[0][0][0][0][0], walk_num*(cord_num+1)*elec_num*elec_num*4); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_electron_coord(context, 'N', &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); rc = qmckl_set_single_point(context, 'N', 2, new_coords, 3); assert (rc == QMCKL_SUCCESS); double een_rescaled_single_e_gl[walk_num][cord_num+1][elec_num][4]; rc = qmckl_get_een_rescaled_single_e_gl(context, &een_rescaled_single_e_gl[0][0][0][0], walk_num*(cord_num+1)*elec_num*4); assert (rc == QMCKL_SUCCESS); coords[0][2][0] = new_coords[0]; coords[0][2][1] = new_coords[1]; coords[0][2][2] = new_coords[2]; rc = qmckl_set_point(context, 'N', elec_num, &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); //rc = qmckl_check(context, // qmckl_set_electron_coord (context, 'N', walk_num, &coords[0][0][0], walk_num*3*elec_num) // ); //assert(rc == QMCKL_SUCCESS); rc = qmckl_get_jastrow_champ_een_rescaled_e_gl(context, &een_rescaled_ee_gl[0][0][0][0][0], walk_num*(cord_num+1)*elec_num*elec_num*4); assert (rc == QMCKL_SUCCESS); for (int l = 0; l < cord_num+1; l++) { for (int nw = 0; nw < walk_num; nw++) { for (int i = 0; i < elec_num; i++) { for (int m = 0; m < 4; m++) { //printf("een_rescaled_ee_gl[nw][l][i][m][2] %i %i %i %f \n", l, m ,i, een_rescaled_ee_gl[nw][l][i][m][2]); //printf("een_rescaled_single_e_gl[nw][l][i][m] %i %i %i %f\n", l, m, i,een_rescaled_single_e_gl[nw][l][i][m]); //if (m == 3) { // assert(fabs(een_rescaled_ee_gl[nw][l][2][m][i] - een_rescaled_single_e_gl[nw][l][m][i]) < 1.e-12); //} else{ // assert(fabs(een_rescaled_ee_gl[nw][l][2][m][i] + een_rescaled_single_e_gl[nw][l][m][i]) < 1.e-12); //} assert(fabs(een_rescaled_ee_gl[nw][l][i][m][2] - een_rescaled_single_e_gl[nw][l][i][m]) < 1.e-12); } } } } #+end_src * gl delta p ** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_jastrow_champ_delta_p_gl(qmckl_context context, double* const delta_p_gl, const int64_t size_max); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_jastrow_champ_delta_p_gl(qmckl_context context, double* const delta_p_gl, const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_exit_code rc; rc = qmckl_provide_jastrow_champ_delta_p_gl(context); if (rc != QMCKL_SUCCESS) return rc; qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int64_t sze = ctx->electron.walker.num * ctx->jastrow_champ.cord_num * (ctx->jastrow_champ.cord_num + 1) * ctx->nucleus.num * ctx->electron.num * 4; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_jastrow_champ_delta_p_gl", "Array too small."); } memcpy(delta_p_gl, ctx->single_point.delta_p_gl, sze * sizeof(double)); return QMCKL_SUCCESS; } #+end_src ** provide #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_provide_jastrow_champ_delta_p_gl(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_provide_jastrow_champ_delta_p_gl(qmckl_context context) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); qmckl_exit_code rc = qmckl_provide_een_rescaled_single_e(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_een_rescaled_single_n(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_een_rescaled_single_e_gl(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_een_rescaled_single_n_gl(context); if(rc != QMCKL_SUCCESS) return rc; /* Compute if necessary */ if (ctx->single_point.date > ctx->single_point.delta_p_gl_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { if (ctx->single_point.delta_p_gl != NULL) { rc = qmckl_free(context, ctx->single_point.delta_p_gl); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_jastrow_champ_delta_p_gl", "Unable to free ctx->single_point.delta_p_gl"); } ctx->single_point.delta_p_gl = NULL; } } /* Allocate array */ if (ctx->single_point.delta_p_gl == NULL) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = 4 * ctx->electron.walker.num * ctx->jastrow_champ.cord_num * (ctx->jastrow_champ.cord_num + 1) * ctx->nucleus.num * ctx->electron.num * sizeof(double); double* delta_p_gl = (double*) qmckl_malloc(context, mem_info); if (delta_p_gl == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_jastrow_champ_delta_p_gl", NULL); } ctx->single_point.delta_p_gl = delta_p_gl; } rc = qmckl_compute_jastrow_champ_delta_p_gl(context, ctx->single_point.num, ctx->electron.walker.num, ctx->electron.num, ctx->nucleus.num, ctx->jastrow_champ.cord_num, ctx->jastrow_champ.een_rescaled_n, ctx->jastrow_champ.een_rescaled_e, ctx->single_point.een_rescaled_single_n, ctx->single_point.een_rescaled_single_e, ctx->jastrow_champ.een_rescaled_n_gl, ctx->jastrow_champ.een_rescaled_e_gl, ctx->single_point.een_rescaled_single_n_gl, ctx->single_point.een_rescaled_single_e_gl, ctx->single_point.delta_p_gl); if (rc != QMCKL_SUCCESS) { return rc; } ctx->single_point.delta_p_gl_date = ctx->single_point.date; } return QMCKL_SUCCESS; } #+end_src ** Compute :PROPERTIES: :Name: qmckl_compute_jastrow_champ_delta_p_gl_doc :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_factor_delta_p_gl_args | Variable | Type | In/Out | Description | |----------------------------+---------------------------------------------------------------------+--------+---------------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~num~ | ~int64_t~ | in | Single point number | | ~walk_num~ | ~int64_t~ | in | Number of walkers | | ~elec_num~ | ~int64_t~ | in | Number of electrons | | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | | ~cord_num~ | ~int64_t~ | in | order of polynomials | | ~een_rescaled_n~ | ~double[walk_num][0:cord_num][nucl_num][elec_num]~ | in | Electron-nucleus rescaled distances | | ~een_rescaled_e~ | ~double[walk_num][0:cord_num][elec_num][elec_num]~ | in | Electron-electron rescaled distances | | ~een_rescaled_single_n~ | ~double[walk_num][0:cord_num][nucl_num]~ | in | Electron-nucleus single rescaled distances | | ~een_rescaled_single_e~ | ~double[walk_num][0:cord_num][elec_num]~ | in | Electron-electron single rescaled distances | | ~een_rescaled_n_gl~ | ~double[walk_num][0:cord_num][nucl_num][4][elec_num]~ | in | Electron-nucleus rescaled distances | | ~een_rescaled_e_gl~ | ~double[walk_num][0:cord_num][elec_num][4][elec_num]~ | in | Electron-electron rescaled distances | | ~een_rescaled_single_n_gl~ | ~double[walk_num][0:cord_num][nucl_num][4]~ | in | Electron-nucleus single rescaled distances | | ~een_rescaled_single_e_gl~ | ~double[walk_num][0:cord_num][elec_num][4]~ | in | Electron-electron single rescaled distances | | ~delta_p_gl~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][4][elec_num]~ | out | Electron-nucleus jastrow | #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_jastrow_champ_delta_p_gl_doc_f( & context, num_in, walk_num, elec_num, nucl_num, cord_num, & een_rescaled_n, een_rescaled_e, een_rescaled_single_n, een_rescaled_single_e, & een_rescaled_n_gl, een_rescaled_e_gl, een_rescaled_single_n_gl, een_rescaled_single_e_gl, delta_p_gl) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: num_in, walk_num, elec_num, cord_num, nucl_num double precision , intent(in) :: een_rescaled_n(elec_num, nucl_num, 0:cord_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_single_n(nucl_num, 0:cord_num, walk_num) double precision , intent(in) :: een_rescaled_single_e(elec_num, 0:cord_num, walk_num) double precision , intent(in) :: een_rescaled_n_gl(elec_num, 4, nucl_num, 0:cord_num, walk_num) double precision , intent(in) :: een_rescaled_e_gl(elec_num, 4, elec_num, 0:cord_num, walk_num) double precision , intent(in) :: een_rescaled_single_n_gl(4, nucl_num, 0:cord_num, walk_num) double precision , intent(in) :: een_rescaled_single_e_gl(4,elec_num, 0:cord_num, walk_num) double precision , intent(out) :: delta_p_gl(elec_num,4,nucl_num,0:cord_num, 0:cord_num-1, walk_num) double precision :: delta_e_gl(4,elec_num, 0:cord_num, walk_num) double precision :: delta_e_gl_2(elec_num, 0:cord_num, walk_num) double precision :: een_rescaled_e_gl_2(elec_num, elec_num, 0:cord_num, walk_num) double precision :: een_rescaled_delta_n(nucl_num, 0:cord_num, walk_num) double precision :: delta_c(nucl_num,0:cord_num, 0:cord_num-1, walk_num) double precision :: delta_c2(elec_num, nucl_num,0:cord_num, 0:cord_num-1, walk_num) integer*8 :: i, a, j, l, k, p, m, n, nw, num double precision :: accu, accu2, cn integer*8 :: LDA, LDB, LDC num = num_in + 1 info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT if (walk_num <= 0) info = QMCKL_INVALID_ARG_3 if (elec_num <= 0) info = QMCKL_INVALID_ARG_4 if (nucl_num <= 0) info = QMCKL_INVALID_ARG_5 if (cord_num < 0) info = QMCKL_INVALID_ARG_6 if (info /= QMCKL_SUCCESS) return if (cord_num == 0) return delta_e_gl(:,:,:,:) = een_rescaled_single_e_gl(:,:,:,:) - een_rescaled_e_gl(num, :, :, :, :) delta_e_gl(:, num, :, :) = 0.0d0 een_rescaled_delta_n(:,:,:) = een_rescaled_single_n(:,:,:) - een_rescaled_n(num, :, :, :) delta_c = 0.0d0 delta_c2 = 0.0d0 delta_p_gl = 0.0d0 if (.false.) then ! Compute using loops do nw=1, walk_num do m=0, cord_num-1 do a = 1, nucl_num do l=0, cord_num do j = 1, elec_num do k = 1, 3 delta_p_gl(num,k,a,l,m,nw) = delta_p_gl(num,k,a,l,m,nw) + & delta_e_gl(k,j,m,nw) * een_rescaled_n(j,a,l,nw) delta_p_gl(j,k,a,l,m,nw) = delta_p_gl(j,k,a,l,m,nw) - & delta_e_gl(k,j,m,nw) * een_rescaled_n(num,a,l,nw) delta_p_gl(j,k,a,l,m,nw) = delta_p_gl(j,k,a,l,m,nw) + & een_rescaled_e_gl(j,k,num,m,nw) * een_rescaled_delta_n(a,l,nw) delta_p_gl(j,k,a,l,m,nw) = delta_p_gl(j,k,a,l,m,nw) - & delta_e_gl(k,j,m,nw) * een_rescaled_delta_n(a,l,nw) end do delta_p_gl(num,4,a,l,m,nw) = delta_p_gl(num,4,a,l,m,nw) + & delta_e_gl(4,j,m,nw) * een_rescaled_n(j,a,l,nw) delta_p_gl(j,4,a,l,m,nw) = delta_p_gl(j,4,a,l,m,nw) + & delta_e_gl(4,j,m,nw) * een_rescaled_n(num,a,l,nw) delta_p_gl(j,4,a,l,m,nw) = delta_p_gl(j,4,a,l,m,nw) + & een_rescaled_e_gl(num,4,j,m,nw) * een_rescaled_delta_n(a,l,nw) delta_p_gl(j,4,a,l,m,nw) = delta_p_gl(j,4,a,l,m,nw) + & delta_e_gl(4,j,m,nw) * een_rescaled_delta_n(a,l,nw) end do end do end do end do end do else ! Use DGEMM do nw=1, walk_num do m=0, cord_num-1 do k = 1, 3 delta_e_gl_2(:,:,:) = delta_e_gl(k, :,:,:) een_rescaled_e_gl_2(:,:,:,:) = een_rescaled_e_gl(:,k, :,:,:) info = qmckl_dgemm(context, 'T', 'N', 1_8, nucl_num * (cord_num+1), elec_num, 1.0d0, & delta_e_gl_2(1,m,nw),elec_num, & een_rescaled_n(1,1,0,nw),elec_num, & 0.0d0, & delta_c(1,0,m,nw),1_8) info = qmckl_dgemm(context, 'N', 'N', elec_num, nucl_num * (cord_num+1), 1_8, -1.0d0, & delta_e_gl_2(1,m,nw),elec_num, & een_rescaled_n(num,1,0,nw),elec_num, & 0.0d0, & delta_c2(1,1,0,m,nw),elec_num) delta_c2(num,:,:,m,nw) = delta_c2(num,:,:,m,nw) + delta_c(:,:,m,nw) delta_p_gl(:,k,:,:,m,nw) = delta_c2(:,:,:,m,nw) info = qmckl_dgemm(context, 'N', 'T', elec_num, nucl_num * (cord_num+1), 1_8, 1.0d0, & een_rescaled_e_gl_2(:,num,m,nw),elec_num, & een_rescaled_delta_n(:,:,nw),nucl_num* (cord_num+1), & 0.0d0, & delta_c2(:,:,:,m,nw),elec_num) delta_p_gl(:,k,:,:,m,nw) = delta_p_gl(:,k,:,:,m,nw) + delta_c2(:,:,:,m,nw) info = qmckl_dgemm(context, 'N', 'T', elec_num, nucl_num * (cord_num+1), 1_8, -1.0d0, & delta_e_gl_2(:,m,nw),elec_num, & een_rescaled_delta_n(:,:,nw),nucl_num* (cord_num+1), & 0.0d0, & delta_c2(:,:,:,m,nw),elec_num) delta_p_gl(:,k,:,:,m,nw) = delta_p_gl(:,k,:,:,m,nw) + delta_c2(:,:,:,m,nw) end do k = 4 delta_e_gl_2(:,:,:) = delta_e_gl(k, :,:,:) een_rescaled_e_gl_2(:,:,:,:) = een_rescaled_e_gl(:,k, :,:,:) info = qmckl_dgemm(context, 'T', 'N', 1_8, nucl_num * (cord_num+1), elec_num, 1.0d0, & delta_e_gl_2(1,m,nw),elec_num, & een_rescaled_n(1,1,0,nw),elec_num, & 0.0d0, & delta_c(1,0,m,nw),1_8) info = qmckl_dgemm(context, 'N', 'N', elec_num, nucl_num * (cord_num+1), 1_8, 1.0d0, & delta_e_gl_2(1,m,nw),elec_num, & een_rescaled_n(num,1,0,nw),elec_num, & 0.0d0, & delta_c2(1,1,0,m,nw),elec_num) delta_c2(num,:,:,m,nw) = delta_c2(num,:,:,m,nw) + delta_c(:,:,m,nw) delta_p_gl(:,k,:,:,m,nw) = delta_c2(:,:,:,m,nw) info = qmckl_dgemm(context, 'N', 'T', elec_num, nucl_num * (cord_num+1), 1_8, 1.0d0, & een_rescaled_e_gl_2(:,num,m,nw),elec_num, & een_rescaled_delta_n(:,:,nw),nucl_num* (cord_num+1), & 0.0d0, & delta_c2(:,:,:,m,nw),elec_num) delta_p_gl(:,k,:,:,m,nw) = delta_p_gl(:,k,:,:,m,nw) + delta_c2(:,:,:,m,nw) info = qmckl_dgemm(context, 'N', 'T', elec_num, nucl_num * (cord_num+1), 1_8, 1.0d0, & delta_e_gl_2(:,m,nw),elec_num, & een_rescaled_delta_n(:,:,nw),nucl_num* (cord_num+1), & 0.0d0, & delta_c2(:,:,:,m,nw),elec_num) delta_p_gl(:,k,:,:,m,nw) = delta_p_gl(:,k,:,:,m,nw) + delta_c2(:,:,:,m,nw) end do end do end if end function qmckl_compute_jastrow_champ_delta_p_gl_doc_f #+end_src # #+CALL: generate_c_header(table=qmckl_factor_een_args,rettyp=qmckl_exit_code),fname=get_value("Name")) #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_jastrow_champ_delta_p_gl_doc (const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, const double* een_rescaled_n, const double* een_rescaled_e, const double* een_rescaled_single_n, const double* een_rescaled_single_e, const double* een_rescaled_n_gl, const double* een_rescaled_e_gl, const double* een_rescaled_single_n_gl, const double* een_rescaled_single_e_gl, double* const delta_p_gl ); qmckl_exit_code qmckl_compute_jastrow_champ_delta_p_gl (const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, const double* een_rescaled_n, const double* een_rescaled_e, const double* een_rescaled_single_n, const double* een_rescaled_single_e, const double* een_rescaled_n_gl, const double* een_rescaled_e_gl, const double* een_rescaled_single_n_gl, const double* een_rescaled_single_e_gl, double* const delta_p_gl ); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_compute_jastrow_champ_delta_p_gl (const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, const double* een_rescaled_n, const double* een_rescaled_e, const double* een_rescaled_single_n, const double* een_rescaled_single_e, const double* een_rescaled_n_gl, const double* een_rescaled_e_gl, const double* een_rescaled_single_n_gl, const double* een_rescaled_single_e_gl, double* const delta_p_gl ) { #ifdef HAVE_HPC return qmckl_compute_jastrow_champ_delta_p_gl_doc #else return qmckl_compute_jastrow_champ_delta_p_gl_doc #endif (context, num, walk_num, elec_num, nucl_num, cord_num, een_rescaled_n, een_rescaled_e, een_rescaled_single_n, een_rescaled_single_e, een_rescaled_n_gl, een_rescaled_e_gl, een_rescaled_single_n_gl, een_rescaled_single_e_gl, delta_p_gl); } #+end_src #+CALL: generate_c_interface(table=qmckl_factor_delta_p_gl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_jastrow_champ_delta_p_gl_doc")) #+RESULTS: #+begin_src f90 :tangle (eval f) :comments org :exports none integer(c_int32_t) function qmckl_compute_jastrow_champ_delta_p_gl_doc & (context, & num, & walk_num, & elec_num, & nucl_num, & cord_num, & een_rescaled_n, & een_rescaled_e, & een_rescaled_single_n, & een_rescaled_single_e, & een_rescaled_n_gl, & een_rescaled_e_gl, & een_rescaled_single_n_gl, & een_rescaled_single_e_gl, & delta_p_gl) & 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 :: num 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) :: een_rescaled_n(elec_num,nucl_num,0:cord_num,walk_num) real (c_double ) , intent(in) :: een_rescaled_e(elec_num,elec_num,0:cord_num,walk_num) real (c_double ) , intent(in) :: een_rescaled_single_n(nucl_num,0:cord_num,walk_num) real (c_double ) , intent(in) :: een_rescaled_single_e(elec_num,0:cord_num,walk_num) real (c_double ) , intent(in) :: een_rescaled_n_gl(elec_num,4,nucl_num,0:cord_num,walk_num) real (c_double ) , intent(in) :: een_rescaled_e_gl(elec_num,4,elec_num,0:cord_num,walk_num) real (c_double ) , intent(in) :: een_rescaled_single_n_gl(4, nucl_num,0:cord_num,walk_num) real (c_double ) , intent(in) :: een_rescaled_single_e_gl(4,elec_num,0:cord_num,walk_num) real (c_double ) , intent(out) :: delta_p_gl(elec_num,4,nucl_num,0:cord_num,0:cord_num-1,walk_num) integer(c_int32_t), external :: qmckl_compute_jastrow_champ_delta_p_gl_doc_f info = qmckl_compute_jastrow_champ_delta_p_gl_doc_f & (context, & num, & walk_num, & elec_num, & nucl_num, & cord_num, & een_rescaled_n, & een_rescaled_e, & een_rescaled_single_n, & een_rescaled_single_e, & een_rescaled_n_gl, & een_rescaled_e_gl, & een_rescaled_single_n_gl, & een_rescaled_single_e_gl, & delta_p_gl) end function qmckl_compute_jastrow_champ_delta_p_gl_doc #+end_src ** test #+begin_src c :tangle (eval c_test) /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); rc = qmckl_set_point(context, 'N', elec_num, elec_coord, walk_num*elec_num*3); assert(rc == QMCKL_SUCCESS); double p_gl_old[walk_num][cord_num][cord_num+1][nucl_num][4][elec_num]; rc = qmckl_get_jastrow_champ_dtmp_c(context, &p_gl_old[0][0][0][0][0][0]); rc = qmckl_get_electron_coord(context, 'N', &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); rc = qmckl_set_single_point(context, 'N', 2, new_coords, 3); assert (rc == QMCKL_SUCCESS); double delta_p_gl[walk_num][cord_num][cord_num+1][nucl_num][4][elec_num]; rc = qmckl_get_jastrow_champ_delta_p_gl(context, &delta_p_gl[0][0][0][0][0][0], 4*walk_num*cord_num*(cord_num+1)*nucl_num*elec_num); assert (rc == QMCKL_SUCCESS); coords[0][2][0] = new_coords[0]; coords[0][2][1] = new_coords[1]; coords[0][2][2] = new_coords[2]; rc = qmckl_set_point(context, 'N', elec_num, &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); double p_gl_new[walk_num][cord_num][cord_num+1][nucl_num][4][elec_num]; rc = qmckl_get_jastrow_champ_dtmp_c(context, &p_gl_new[0][0][0][0][0][0]); assert (rc == QMCKL_SUCCESS); for (int nw = 0; nw < walk_num; nw++){ for (int l = 0; l < cord_num; l++){ for (int m = 0; m <= cord_num; m++){ for (int a = 0; a < nucl_num; a++) { for (int i = 0; i < elec_num; i++){ for (int k = 0; k < 4; k++){ if (fabs(((p_gl_new[nw][l][m][a][k][i]-p_gl_old[nw][l][m][a][k][i])-delta_p_gl[nw][l][m][a][k][i])) > 1.e-12) { printf("p_gl[%d][%d][%d][%d][%d][%d] = %f\n", nw, l, m, a, k, i, p_gl_new[nw][l][m][a][k][i] - p_gl_old[nw][l][m][a][k][i]); printf("delta_p_gl[%d][%d][%d][%d][%d][%d] = %f\n", nw, l, m, a, k, i, delta_p_gl[nw][l][m][a][k][i]); } assert(fabs(((p_gl_new[nw][l][m][a][k][i]-p_gl_old[nw][l][m][a][k][i])-delta_p_gl[nw][l][m][a][k][i])) < 1.e-12); } } } } } } //assert(0); #+end_src * Delta grad e-e-n ** get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_jastrow_champ_single_een_gl(qmckl_context context, double* const delta_een_gl, const int64_t size_max); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_jastrow_champ_single_een_gl(qmckl_context context, double* const delta_een_gl, const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_exit_code rc; rc = qmckl_provide_jastrow_champ_single_een_gl(context); if (rc != QMCKL_SUCCESS) return rc; qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int64_t sze = 4 * ctx->electron.num * ctx->electron.walker.num; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_jastrow_champ_single_een_gl", "Array too small. Expected 4 * ctx->electron.num * ctx->electron.walker.num"); } memcpy(delta_een_gl, ctx->single_point.delta_een_gl, sze * sizeof(double)); return QMCKL_SUCCESS; } #+end_src #+begin_src f90 :tangle (eval fh_func) :comments org interface integer(qmckl_exit_code) function qmckl_get_jastrow_champ_single_een_gl (context, & delta_een_gl, size_max) bind(C) use, intrinsic :: iso_c_binding import implicit none integer (qmckl_context) , intent(in), value :: context integer(c_int64_t), intent(in), value :: size_max real(c_double), intent(out) :: delta_een_gl(size_max) end function end interface #+end_src ** Provide :noexport: #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_provide_jastrow_champ_single_een_gl(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_provide_jastrow_champ_single_een_gl(qmckl_context context) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); qmckl_exit_code rc = qmckl_provide_jastrow_champ_delta_p(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_jastrow_champ_factor_een(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_jastrow_champ_delta_p_gl(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_jastrow_champ_factor_een_gl(context); if(rc != QMCKL_SUCCESS) return rc; /* Compute if necessary */ if (ctx->single_point.date > ctx->single_point.delta_een_gl_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { if (ctx->single_point.delta_een_gl != NULL) { rc = qmckl_free(context, ctx->single_point.delta_een_gl); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_jastrow_champ_single_een_gl", "Unable to free ctx->single_point.delta_een_gl"); } ctx->single_point.delta_een_gl = NULL; } } /* Allocate array */ if (ctx->single_point.delta_een_gl == NULL) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = 4 * ctx->electron.num * ctx->electron.walker.num * sizeof(double); double* delta_een_gl = (double*) qmckl_malloc(context, mem_info); if (delta_een_gl == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_jastrow_champ_single_een_gl", NULL); } ctx->single_point.delta_een_gl = delta_een_gl; } rc = qmckl_compute_jastrow_champ_factor_single_een_gl_doc(context, ctx->single_point.num, 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.tmp_c, ctx->jastrow_champ.dtmp_c, ctx->single_point.delta_p, ctx->single_point.delta_p_gl, ctx->jastrow_champ.een_rescaled_n, ctx->single_point.een_rescaled_single_n, ctx->jastrow_champ.een_rescaled_n_gl, ctx->single_point.een_rescaled_single_n_gl, ctx->single_point.delta_een_gl); if (rc != QMCKL_SUCCESS) { return rc; } ctx->single_point.delta_een_gl_date = ctx->single_point.date; } return QMCKL_SUCCESS; } #+end_src ** Compute :PROPERTIES: :Name: qmckl_compute_jastrow_champ_factor_single_een_gl_doc :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_factor_single_een_gl_args | Variable | Type | In/Out | Description | |----------------------------+---------------------------------------------------------------------+--------+--------------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~num~ | ~int64_t~ | in | Single point number | | ~walk_num~ | ~int64_t~ | in | Number of walkers | | ~elec_num~ | ~int64_t~ | in | Number of electrons | | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | | ~cord_num~ | ~int64_t~ | in | order of polynomials | | ~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 | | ~tmp_c~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]~ | in | vector of non-zero coefficients | | ~dtmp_c~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][4][elec_num]~ | in | vector of non-zero coefficients | | ~delta_p~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]~ | in | vector of non-zero coefficients | | ~delta_p_gl~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][4][elec_num]~ | in | vector of non-zero coefficients | | ~een_rescaled_n~ | ~double[walk_num][0:cord_num][nucl_num][elec_num]~ | in | Electron-nucleus rescaled distances | | ~een_rescaled_single_n~ | ~double[walk_num][0:cord_num][nucl_num]~ | in | Electron-nucleus single rescaled distances | | ~een_rescaled_n_gl~ | ~double[walk_num][0:cord_num][nucl_num][4][elec_num]~ | in | Electron-nucleus rescaled distances | | ~een_rescaled_single_n_gl~ | ~double[walk_num][0:cord_num][nucl_num][4]~ | in | Electron-nucleus single rescaled distances | | ~delta_een_gl~ | ~double[walk_num][elec_num][4]~ | out | Electron-nucleus jastrow | #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_jastrow_champ_factor_single_een_gl_doc_f( & context, num_in, walk_num, elec_num, nucl_num, cord_num, & dim_c_vector, c_vector_full, lkpm_combined_index, & tmp_c, dtmp_c, delta_p, delta_p_gl, een_rescaled_n, een_rescaled_single_n, & een_rescaled_n_gl, een_rescaled_single_n_gl, delta_een_gl) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: num_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) :: tmp_c(elec_num, nucl_num,0:cord_num, 0:cord_num-1, walk_num) double precision , intent(in) :: dtmp_c(elec_num, 4, nucl_num,0:cord_num, 0:cord_num-1, walk_num) double precision , intent(in) :: delta_p(elec_num, nucl_num,0:cord_num, 0:cord_num-1, walk_num) double precision , intent(in) :: delta_p_gl(elec_num, 4, nucl_num,0:cord_num, 0:cord_num-1, walk_num) double precision , intent(in) :: een_rescaled_n(elec_num, nucl_num, 0:cord_num, walk_num) double precision , intent(in) :: een_rescaled_single_n(nucl_num, 0:cord_num, walk_num) double precision , intent(in) :: een_rescaled_n_gl(elec_num, 4, nucl_num, 0:cord_num, walk_num) double precision , intent(in) :: een_rescaled_single_n_gl(4, nucl_num, 0:cord_num, walk_num) double precision , intent(out) :: delta_een_gl(4, elec_num, walk_num) integer*8 :: i, a, j, l, k, p, m, n, nw, kk, num double precision :: accu, accu2, cn integer*8 :: LDA, LDB, LDC double precision :: een_rescaled_delta_n_gl(4, nucl_num, 0:cord_num, walk_num) double precision :: een_rescaled_delta_n(nucl_num, 0:cord_num, walk_num) num = num_in + 1 info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT if (walk_num <= 0) info = QMCKL_INVALID_ARG_3 if (elec_num <= 0) info = QMCKL_INVALID_ARG_4 if (nucl_num <= 0) info = QMCKL_INVALID_ARG_5 if (cord_num < 0) info = QMCKL_INVALID_ARG_6 if (info /= QMCKL_SUCCESS) return delta_een_gl = 0.0d0 een_rescaled_delta_n(:,:,:) = een_rescaled_single_n(:,:,:) - een_rescaled_n(num, :, :, :) een_rescaled_delta_n_gl(:,:,:,:) = een_rescaled_single_n_gl(:,:,:,:) - een_rescaled_n_gl(num, :,:,:,:) if (cord_num == 0) return 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) do a = 1, nucl_num cn = c_vector_full(a, n) if(cn == 0.d0) cycle do i = 1, elec_num do kk = 1, 4 delta_een_gl(kk,i,nw) = delta_een_gl(kk,i,nw) + ( & delta_p_gl(i,kk,a,m ,k,nw) * een_rescaled_n(i,a,m+l,nw) + & delta_p_gl(i,kk,a,m+l,k,nw) * een_rescaled_n(i,a,m ,nw) + & delta_p(i,a,m ,k,nw) * een_rescaled_n_gl(i,kk,a,m+l,nw) + & delta_p(i,a,m+l,k,nw) * een_rescaled_n_gl(i,kk,a,m ,nw) ) * cn end do end do do kk = 1, 4 delta_een_gl(kk,num,nw) = delta_een_gl(kk,num,nw) + ( & (dtmp_c(num,kk,a,m ,k,nw) + delta_p_gl(num,kk,a,m ,k,nw)) * een_rescaled_delta_n(a,m+l,nw) + & (dtmp_c(num,kk,a,m+l,k,nw) + delta_p_gl(num,kk,a,m+l,k,nw)) * een_rescaled_delta_n(a,m ,nw) + & (tmp_c(num,a,m ,k,nw) + delta_p(num,a,m ,k,nw)) * een_rescaled_delta_n_gl(kk,a,m+l,nw) + & (tmp_c(num,a,m+l,k,nw) + delta_p(num,a,m+l,k,nw)) * een_rescaled_delta_n_gl(kk,a,m ,nw) ) & ,* cn end do cn = cn + cn do i = 1, elec_num delta_een_gl(4,i,nw) = delta_een_gl(4,i,nw) + ( & delta_p_gl(i,1,a,m ,k,nw) * een_rescaled_n_gl(i,1,a,m+l,nw) + & delta_p_gl(i,1,a,m+l,k,nw) * een_rescaled_n_gl(i,1,a,m ,nw) + & delta_p_gl(i,2,a,m ,k,nw) * een_rescaled_n_gl(i,2,a,m+l,nw) + & delta_p_gl(i,2,a,m+l,k,nw) * een_rescaled_n_gl(i,2,a,m ,nw) + & delta_p_gl(i,3,a,m ,k,nw) * een_rescaled_n_gl(i,3,a,m+l,nw) + & delta_p_gl(i,3,a,m+l,k,nw) * een_rescaled_n_gl(i,3,a,m ,nw) ) * cn end do delta_een_gl(4,num,nw) = delta_een_gl(4,num,nw) + ( & (delta_p_gl(num,1,a,m ,k,nw) + dtmp_c(num,1,a,m ,k,nw)) * een_rescaled_delta_n_gl(1,a,m+l,nw) + & (delta_p_gl(num,1,a,m+l,k,nw) + dtmp_c(num,1,a,m+l,k,nw)) * een_rescaled_delta_n_gl(1,a,m ,nw) + & (delta_p_gl(num,2,a,m ,k,nw) + dtmp_c(num,2,a,m ,k,nw)) * een_rescaled_delta_n_gl(2,a,m+l,nw) + & (delta_p_gl(num,2,a,m+l,k,nw) + dtmp_c(num,2,a,m+l,k,nw)) * een_rescaled_delta_n_gl(2,a,m ,nw) + & (delta_p_gl(num,3,a,m ,k,nw) + dtmp_c(num,3,a,m ,k,nw)) * een_rescaled_delta_n_gl(3,a,m+l,nw) + & (delta_p_gl(num,3,a,m+l,k,nw) + dtmp_c(num,3,a,m+l,k,nw)) * een_rescaled_delta_n_gl(3,a,m ,nw) ) * cn end do end do end do end function qmckl_compute_jastrow_champ_factor_single_een_gl_doc_f #+end_src # #+CALL: generate_c_header(table=qmckl_factor_een_args,rettyp=qmckl_exit_code),fname=get_value("Name")) #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_jastrow_champ_factor_single_een_gl_doc (const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, const int64_t dim_c_vector, const double* c_vector_full, const int64_t* lkpm_combined_index, const double* tmp_c, const double* dtmp_c, const double* delta_p, const double* delta_p_gl, const double* een_rescaled_n, const double* een_rescaled_single_n, const double* een_rescaled_n_gl, const double* een_rescaled_single_n_gl, double* const delta_een_gl ); qmckl_exit_code qmckl_compute_jastrow_champ_factor_single_een_gl (const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, const int64_t dim_c_vector, const double* c_vector_full, const int64_t* lkpm_combined_index, const double* tmp_c, const double* dtmp_c, const double* delta_p, const double* delta_p_gl, const double* een_rescaled_n, const double* een_rescaled_single_n, const double* een_rescaled_n_gl, const double* een_rescaled_single_n_gl, double* const delta_een_gl ); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_compute_jastrow_champ_factor_single_een_gl (const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, const int64_t dim_c_vector, const double* c_vector_full, const int64_t* lkpm_combined_index, const double* tmp_c, const double* dtmp_c, const double* delta_p, const double* delta_p_gl, const double* een_rescaled_n, const double* een_rescaled_single_n, const double* een_rescaled_n_gl, const double* een_rescaled_single_n_gl, double* const delta_een_gl ) { #ifdef HAVE_HPC return qmckl_compute_jastrow_champ_factor_single_een_gl_doc #else return qmckl_compute_jastrow_champ_factor_single_een_gl_doc #endif (context, num, walk_num, elec_num, nucl_num, cord_num, dim_c_vector, c_vector_full, lkpm_combined_index, tmp_c, dtmp_c, delta_p, delta_p_gl, een_rescaled_n, een_rescaled_single_n, een_rescaled_n_gl, een_rescaled_single_n_gl, delta_een_gl ); } #+end_src #+CALL: generate_c_interface(table=qmckl_factor_single_een_gl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_jastrow_champ_factor_single_een_gl_doc")) #+RESULTS: #+begin_src f90 :tangle (eval f) :comments org :exports none integer(c_int32_t) function qmckl_compute_jastrow_champ_factor_single_een_gl_doc & (context, & num, & walk_num, & elec_num, & nucl_num, & cord_num, & dim_c_vector, & c_vector_full, & lkpm_combined_index, & tmp_c, & dtmp_c, & delta_p, & delta_p_gl, & een_rescaled_n, & een_rescaled_single_n, & een_rescaled_n_gl, & een_rescaled_single_n_gl, & delta_een_gl) & 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 :: num 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 integer (c_int64_t) , intent(in) , value :: dim_c_vector real (c_double ) , intent(in) :: c_vector_full(nucl_num,dim_c_vector) integer (c_int64_t) , intent(in) :: lkpm_combined_index(dim_c_vector,4) real (c_double ) , intent(in) :: tmp_c(elec_num,nucl_num,0:cord_num,0:cord_num-1,walk_num) real (c_double ) , intent(in) :: dtmp_c(elec_num,4,nucl_num,0:cord_num,0:cord_num-1,walk_num) real (c_double ) , intent(in) :: delta_p(elec_num,nucl_num,0:cord_num,0:cord_num-1,walk_num) real (c_double ) , intent(in) :: delta_p_gl(elec_num,4,nucl_num,0:cord_num,0:cord_num-1,walk_num) real (c_double ) , intent(in) :: een_rescaled_n(elec_num,nucl_num,0:cord_num,walk_num) real (c_double ) , intent(in) :: een_rescaled_single_n(nucl_num,0:cord_num,walk_num) real (c_double ) , intent(in) :: een_rescaled_n_gl(elec_num,4,nucl_num,0:cord_num,walk_num) real (c_double ) , intent(in) :: een_rescaled_single_n_gl(4,nucl_num,0:cord_num,walk_num) real (c_double ) , intent(out) :: delta_een_gl(4, elec_num, walk_num) integer(c_int32_t), external :: qmckl_compute_jastrow_champ_factor_single_een_gl_doc_f info = qmckl_compute_jastrow_champ_factor_single_een_gl_doc_f & (context, & num, & walk_num, & elec_num, & nucl_num, & cord_num, & dim_c_vector, & c_vector_full, & lkpm_combined_index, & tmp_c, & dtmp_c, & delta_p, & delta_p_gl, & een_rescaled_n, & een_rescaled_single_n, & een_rescaled_n_gl, & een_rescaled_single_n_gl, & delta_een_gl) end function qmckl_compute_jastrow_champ_factor_single_een_gl_doc #+end_src ** test #+begin_src c :tangle (eval c_test) /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); rc = qmckl_set_point(context, 'N', elec_num, elec_coord, walk_num*elec_num*3); assert(rc == QMCKL_SUCCESS); double een_gl_old[walk_num][4][elec_num]; rc = qmckl_get_jastrow_champ_factor_een_gl(context, &een_gl_old[0][0][0], walk_num*elec_num*4); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_electron_coord(context, 'N', &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); rc = qmckl_set_single_point(context, 'N', 2, new_coords, 3); assert (rc == QMCKL_SUCCESS); double delta_een_gl[walk_num][elec_num][4]; rc = qmckl_get_jastrow_champ_single_een_gl(context, &delta_een_gl[0][0][0], walk_num*elec_num*4); assert (rc == QMCKL_SUCCESS); coords[0][2][0] = new_coords[0]; coords[0][2][1] = new_coords[1]; coords[0][2][2] = new_coords[2]; rc = qmckl_set_point(context, 'N', elec_num, &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); double een_gl_new[walk_num][4][elec_num]; rc = qmckl_get_jastrow_champ_factor_een_gl(context, &een_gl_new[0][0][0], walk_num*elec_num*4); assert (rc == QMCKL_SUCCESS); for (int nw = 0; nw < walk_num; nw++) { for (int i = 0; i < elec_num; i++) { for (int m = 0; m < 4; m++) { //printf("delta_een_gl[%d][%d][%d] = %f\n", nw, i, m, delta_een_gl[nw][i][m]); //printf("een_gl_[%d][%d][%d] = %f\n", nw, m,i, een_gl_new[nw][m][i]-een_gl_old[nw][m][i]); assert(fabs((een_gl_new[nw][m][i]- een_gl_old[nw][m][i]) - delta_een_gl[nw][i][m]) < 1.e-12); } } } #+end_src * ee distance rescaled single point gl ** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_ee_rescaled_single_gl(qmckl_context context, double* const distance_rescaled_gl); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_ee_rescaled_single_gl(qmckl_context context, double* const distance_rescaled_gl) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_exit_code rc; rc = qmckl_provide_ee_rescaled_single_gl(context); if (rc != QMCKL_SUCCESS) return rc; qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); size_t sze = 4 * ctx->electron.num * ctx->electron.walker.num; memcpy(distance_rescaled_gl, ctx->single_point.ee_rescaled_single_gl, 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_rescaled_single_gl(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_provide_ee_rescaled_single_gl(qmckl_context context) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); qmckl_exit_code rc = qmckl_provide_single_ee_distance(context); if(rc != QMCKL_SUCCESS) return rc; /* Compute if necessary */ if (ctx->single_point.date > ctx->single_point.ee_rescaled_single_gl_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { if (ctx->single_point.ee_rescaled_single_gl != NULL) { rc = qmckl_free(context, ctx->single_point.ee_rescaled_single_gl); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_ee_rescaled_single_gl", "Unable to free ctx->single_point.ee_rescaled_single_gl"); } ctx->single_point.ee_rescaled_single_gl = NULL; } } /* Allocate array */ if (ctx->single_point.ee_rescaled_single_gl == NULL) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = 4 * ctx->electron.num * ctx->electron.walker.num * sizeof(double); double* ee_rescaled_single_gl = (double*) qmckl_malloc(context, mem_info); if (ee_rescaled_single_gl == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_ee_rescaled_single_gl", NULL); } ctx->single_point.ee_rescaled_single_gl = ee_rescaled_single_gl; } qmckl_exit_code rc = qmckl_compute_ee_rescaled_single_gl(context, ctx->single_point.num, ctx->electron.num, ctx->jastrow_champ.rescale_factor_ee, ctx->electron.walker.num, ctx->single_point.single_ee_distance, ctx->electron.walker.point.coord.data, ctx->single_point.coord.data, ctx->single_point.ee_rescaled_single_gl); if (rc != QMCKL_SUCCESS) { return rc; } ctx->single_point.ee_rescaled_single_gl_date = ctx->single_point.date; } return QMCKL_SUCCESS; } #+end_src ** Compute :PROPERTIES: :Name: qmckl_compute_ee_rescaled_single_gl :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_ee_rescaled_single_gl_args | Variable | Type | In/Out | Description | |-------------------------+---------------------------------+--------+-------------------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~num~ | ~int64_t~ | in | Number of electrons | | ~elec_num~ | ~int64_t~ | in | Number of electrons | | ~rescale_factor_ee~ | ~double~ | in | Factor to rescale ee distances | | ~walk_num~ | ~int64_t~ | in | Number of walkers | | ~single_ee_distance~ | ~double[elec_num][walk_num]~ | in | Electron coordinates | | ~elec_coord~ | ~double[3][walk_num][elec_num]~ | in | Electron coordinates | | ~coord~ | ~double[3]~ | in | Electron coordinates | | ~ee_rescaled_single_gl~ | ~double[walk_num][elec_num][4]~ | out | Electron-electron rescaled distance derivatives | #+begin_src f90 :comments org :tangle (eval f) :noweb yes function qmckl_compute_ee_rescaled_single_gl_doc(context, num_in, & elec_num, rescale_factor_ee, walk_num, single_ee_distance, elec_coord, coord, ee_rescaled_single_gl) & bind(C) result(info) use qmckl implicit none integer(qmckl_context), intent(in), value :: context integer (c_int64_t) , intent(in) , value :: num_in integer (c_int64_t) , intent(in) , value :: elec_num real (c_double ) , intent(in) , value :: rescale_factor_ee integer (c_int64_t) , intent(in) , value :: walk_num real (c_double ) , intent(in) :: single_ee_distance(elec_num,walk_num) real (c_double ) , intent(in) :: elec_coord(elec_num,walk_num,3) real (c_double ) , intent(in) :: coord(3) real (c_double ) , intent(out) :: ee_rescaled_single_gl(4,elec_num,walk_num) integer(qmckl_exit_code) :: info integer*8 :: nw, i, ii, num double precision :: rij_inv, elel_dist_gl(4, elec_num), kappa_l num = num_in + 1 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 ee_rescaled_single_gl = 0.0d0 do nw = 1, walk_num ! prepare the actual een table do i = 1, elec_num rij_inv = 1.0d0 / single_ee_distance(i, nw) do ii = 1, 3 elel_dist_gl(ii, i) = (elec_coord(i,nw, ii) - coord(ii)) * rij_inv end do elel_dist_gl(4, i) = 2.0d0 * rij_inv end do do i = 1, elec_num kappa_l = -1 * rescale_factor_ee ee_rescaled_single_gl(1, i, nw) = elel_dist_gl(1, i) ee_rescaled_single_gl(2, i, nw) = elel_dist_gl(2, i) ee_rescaled_single_gl(3, i, nw) = elel_dist_gl(3, i) ee_rescaled_single_gl(4, i, nw) = elel_dist_gl(4, i) ee_rescaled_single_gl(4, i, nw) = ee_rescaled_single_gl(4, i, nw) & + ee_rescaled_single_gl(1, i, nw) * ee_rescaled_single_gl(1, i, nw) * kappa_l & + ee_rescaled_single_gl(2, i, nw) * ee_rescaled_single_gl(2, i, nw) * kappa_l & + ee_rescaled_single_gl(3, i, nw) * ee_rescaled_single_gl(3, i, nw) * kappa_l ee_rescaled_single_gl(1, i, nw) = ee_rescaled_single_gl(1, i, nw) * dexp(kappa_l * single_ee_distance(i,nw)) ee_rescaled_single_gl(2, i, nw) = ee_rescaled_single_gl(2, i, nw) * dexp(kappa_l * single_ee_distance(i,nw)) ee_rescaled_single_gl(3, i, nw) = ee_rescaled_single_gl(3, i, nw) * dexp(kappa_l * single_ee_distance(i,nw)) ee_rescaled_single_gl(4, i, nw) = ee_rescaled_single_gl(4, i, nw) * dexp(kappa_l * single_ee_distance(i,nw)) end do ee_rescaled_single_gl(1, num, nw) = 0.0d0 ee_rescaled_single_gl(2, num, nw) = 0.0d0 ee_rescaled_single_gl(3, num, nw) = 0.0d0 ee_rescaled_single_gl(4, num, nw) = 0.0d0 end do end function qmckl_compute_ee_rescaled_single_gl_doc #+end_src #+begin_src c :tangle (eval h_private_func) :comments org :exports none qmckl_exit_code qmckl_compute_ee_rescaled_single_gl_doc ( const qmckl_context context, const int64_t num, const int64_t elec_num, const double rescale_factor_ee, const int64_t walk_num, const double* single_ee_distance, const double* elec_coord, const double* coord, double* const ee_rescaled_single_gl ); qmckl_exit_code qmckl_compute_ee_rescaled_single_gl ( const qmckl_context context, const int64_t num, const int64_t elec_num, const double rescale_factor_ee, const int64_t walk_num, const double* single_ee_distance, const double* elec_coord, const double* coord, double* const ee_rescaled_single_gl ); #+end_src #+begin_src c :tangle (eval c) :comments org :exports none qmckl_exit_code qmckl_compute_ee_rescaled_single_gl ( const qmckl_context context, const int64_t num, const int64_t elec_num, const double rescale_factor_ee, const int64_t walk_num, const double* single_ee_distance, const double* elec_coord, const double* coord, double* const ee_rescaled_single_gl ) { #ifdef HAVE_HPC return qmckl_compute_ee_rescaled_single_gl_doc #else return qmckl_compute_ee_rescaled_single_gl_doc #endif (context, num, elec_num, rescale_factor_ee, walk_num, single_ee_distance, elec_coord, coord, ee_rescaled_single_gl); } #+end_src ** test #+begin_src c :tangle (eval c_test) /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); rc = qmckl_set_point(context, 'N', elec_num, elec_coord, walk_num*elec_num*3); assert(rc == QMCKL_SUCCESS); double ee_rescaled_gl[walk_num][elec_num][elec_num][4]; rc = qmckl_get_jastrow_champ_ee_distance_rescaled_gl(context, &ee_rescaled_gl[0][0][0][0]); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_electron_coord(context, 'N', &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); rc = qmckl_set_single_point(context, 'N', 2, new_coords, 3); assert (rc == QMCKL_SUCCESS); double single_ee_rescaled_gl[walk_num][elec_num][4]; rc = qmckl_get_ee_rescaled_single_gl(context, &single_ee_rescaled_gl[0][0][0]); assert (rc == QMCKL_SUCCESS); coords[0][2][0] = new_coords[0]; coords[0][2][1] = new_coords[1]; coords[0][2][2] = new_coords[2]; rc = qmckl_set_point(context, 'N', elec_num, &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_jastrow_champ_ee_distance_rescaled_gl(context, &ee_rescaled_gl[0][0][0][0]); assert (rc == QMCKL_SUCCESS); for (int nw = 0; nw < walk_num; nw++) { for (int i = 0; i < elec_num; i++) { for (int m = 0; m < 4; m++) { if (i == 2) continue; //printf("%f\n", ee_rescaled_gl[nw][2][i][m]); //printf("%f\n", single_ee_rescaled_gl[nw][i][m]); assert(fabs(ee_rescaled_gl[nw][2][i][m] - single_ee_rescaled_gl[nw][i][m]) < 1.e-12); } } } #+end_src * en distance rescaled single point gl ** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_en_rescaled_single_gl(qmckl_context context, double* distance_rescaled_gl); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_en_rescaled_single_gl(qmckl_context context, double* distance_rescaled_gl) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_exit_code rc; rc = qmckl_provide_en_rescaled_single_gl(context); if (rc != QMCKL_SUCCESS) return rc; qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); size_t sze = 4 * ctx->nucleus.num * ctx->electron.walker.num; memcpy(distance_rescaled_gl, ctx->single_point.en_rescaled_single_gl, sze * sizeof(double)); return QMCKL_SUCCESS; } #+end_src ** Provide #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_provide_en_rescaled_single_gl(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_provide_en_rescaled_single_gl(qmckl_context context) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); if (!(ctx->nucleus.provided)) { return QMCKL_NOT_PROVIDED; } qmckl_exit_code rc = qmckl_provide_single_en_distance(context); if(rc != QMCKL_SUCCESS) return rc; /* Compute if necessary */ if (ctx->single_point.date > ctx->single_point.en_rescaled_single_gl_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { if (ctx->single_point.en_rescaled_single_gl != NULL) { rc = qmckl_free(context, ctx->single_point.en_rescaled_single_gl); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_en_rescaled_single_gl", "Unable to free ctx->single_point.en_rescaled_single_gl"); } ctx->single_point.en_rescaled_single_gl = NULL; } } /* Allocate array */ if (ctx->single_point.en_rescaled_single_gl == NULL) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = 4 * ctx->nucleus.num * ctx->electron.walker.num * sizeof(double); double* en_rescaled_single_gl = (double*) qmckl_malloc(context, mem_info); if (en_rescaled_single_gl == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_en_rescaled_single_gl", NULL); } ctx->single_point.en_rescaled_single_gl = en_rescaled_single_gl; } qmckl_exit_code rc = qmckl_compute_en_rescaled_single_gl(context, ctx->nucleus.num, ctx->jastrow_champ.type_nucl_num, ctx->jastrow_champ.type_nucl_vector, ctx->jastrow_champ.rescale_factor_en, ctx->electron.walker.num, ctx->single_point.single_en_distance, ctx->single_point.coord.data, ctx->nucleus.coord.data, ctx->single_point.en_rescaled_single_gl); if (rc != QMCKL_SUCCESS) { return rc; } ctx->single_point.en_rescaled_single_gl_date = ctx->single_point.date; } return QMCKL_SUCCESS; } #+end_src ** Compute :PROPERTIES: :Name: qmckl_compute_en_rescaled_single_gl :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_en_rescaled_single_gl_args | Variable | Type | In/Out | Description | |-------------------------+---------------------------------+--------+---------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | | ~type_nucl_num~ | ~int64_t~ | in | Number of nucleus types | | ~type_nucl_vector~ | ~int64_t[nucl_num]~ | in | Array of nucleus types | | ~rescale_factor_en~ | ~double[nucl_num]~ | in | The factors for rescaled distances | | ~walk_num~ | ~int64_t~ | in | Number of walkers | | ~single_en_distance~ | ~double[walk_num][nucl_num]~ | in | Electron coordinates | | ~coord~ | ~double[3]~ | in | Electron coordinates | | ~nucl_coord~ | ~double[3][nucl_num]~ | in | Electron coordinates | | ~en_rescaled_single_gl~ | ~double[walk_num][nucl_num][4]~ | out | Electron-nucleus distance derivatives | #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_en_rescaled_single_gl_doc_f(context, nucl_num, & type_nucl_num, type_nucl_vector, rescale_factor_en, walk_num, & single_en_distance, coord, nucl_coord, en_rescaled_single_gl) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: nucl_num integer*8 , intent(in) :: type_nucl_num integer*8 , intent(in) :: type_nucl_vector(nucl_num) double precision , intent(in) :: rescale_factor_en(nucl_num) integer*8 , intent(in) :: walk_num double precision , intent(in) :: single_en_distance(nucl_num, walk_num) double precision , intent(in) :: coord(3) double precision , intent(in) :: nucl_coord(nucl_num,3) double precision , intent(out) :: en_rescaled_single_gl(4,nucl_num,walk_num) integer*8 :: nw, a, ii double precision :: ria_inv, elnuc_dist_gl(4, nucl_num), kappa_l info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif if (nucl_num <= 0) then info = QMCKL_INVALID_ARG_2 return endif if (walk_num <= 0) then info = QMCKL_INVALID_ARG_4 return endif en_rescaled_single_gl = 0.0d0 do nw = 1, walk_num ! prepare the actual een table do a = 1, nucl_num ria_inv = 1.0d0 / single_en_distance(a, nw) do ii = 1, 3 elnuc_dist_gl(ii, a) = (coord(ii) - nucl_coord(a, ii)) * ria_inv end do elnuc_dist_gl(4, a) = 2.0d0 * ria_inv end do do a = 1, nucl_num kappa_l = -1 * rescale_factor_en(type_nucl_vector(a)+1) en_rescaled_single_gl(1, a, nw) = elnuc_dist_gl(1, a) en_rescaled_single_gl(2, a, nw) = elnuc_dist_gl(2, a) en_rescaled_single_gl(3, a, nw) = elnuc_dist_gl(3, a) en_rescaled_single_gl(4, a, nw) = elnuc_dist_gl(4, a) en_rescaled_single_gl(4, a, nw) = en_rescaled_single_gl(4, a, nw) & + en_rescaled_single_gl(1, a, nw) * en_rescaled_single_gl(1, a, nw) * kappa_l & + en_rescaled_single_gl(2, a, nw) * en_rescaled_single_gl(2, a, nw) * kappa_l & + en_rescaled_single_gl(3, a, nw) * en_rescaled_single_gl(3, a, nw) * kappa_l en_rescaled_single_gl(1, a, nw) = en_rescaled_single_gl(1, a, nw) * dexp(kappa_l * single_en_distance(a,nw)) en_rescaled_single_gl(2, a, nw) = en_rescaled_single_gl(2, a, nw) * dexp(kappa_l * single_en_distance(a,nw)) en_rescaled_single_gl(3, a, nw) = en_rescaled_single_gl(3, a, nw) * dexp(kappa_l * single_en_distance(a,nw)) en_rescaled_single_gl(4, a, nw) = en_rescaled_single_gl(4, a, nw) * dexp(kappa_l * single_en_distance(a,nw)) end do end do end function qmckl_compute_en_rescaled_single_gl_doc_f #+end_src #+begin_src c :tangle (eval h_private_func) :comments org :exports none qmckl_exit_code qmckl_compute_en_rescaled_single_gl_doc ( const qmckl_context context, const int64_t nucl_num, const int64_t type_nucl_num, int64_t* const type_nucl_vector, const double* rescale_factor_en, const int64_t walk_num, const double* single_en_distance, const double* coord, const double* nucl_coord, double* const en_rescaled_single_gl ); qmckl_exit_code qmckl_compute_en_rescaled_single_gl ( const qmckl_context context, const int64_t nucl_num, const int64_t type_nucl_num, int64_t* const type_nucl_vector, const double* rescale_factor_en, const int64_t walk_num, const double* single_en_distance, const double* coord, const double* nucl_coord, double* const en_rescaled_single_gl ); #+end_src #+begin_src c :tangle (eval c) :comments org :exports none qmckl_exit_code qmckl_compute_en_rescaled_single_gl ( const qmckl_context context, const int64_t nucl_num, const int64_t type_nucl_num, int64_t* const type_nucl_vector, const double* rescale_factor_en, const int64_t walk_num, const double* single_en_distance, const double* coord, const double* nucl_coord, double* const en_rescaled_single_gl ) { #ifdef HAVE_HPC return qmckl_compute_en_rescaled_single_gl_doc #else return qmckl_compute_en_rescaled_single_gl_doc #endif (context, nucl_num, type_nucl_num, type_nucl_vector, rescale_factor_en, walk_num, single_en_distance, coord, nucl_coord, en_rescaled_single_gl ); } #+end_src #+CALL: generate_c_interface(table=qmckl_en_rescaled_single_gl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_en_rescaled_single_gl_doc") #+RESULTS: #+begin_src f90 :tangle (eval f) :comments org :exports none integer(c_int32_t) function qmckl_compute_en_rescaled_single_gl_doc & (context, & nucl_num, & type_nucl_num, & type_nucl_vector, & rescale_factor_en, & walk_num, & single_en_distance, & coord, & nucl_coord, & en_rescaled_single_gl) & 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 :: nucl_num integer (c_int64_t) , intent(in) , value :: type_nucl_num integer (c_int64_t) , intent(in) :: type_nucl_vector(nucl_num) real (c_double ) , intent(in) :: rescale_factor_en(nucl_num) integer (c_int64_t) , intent(in) , value :: walk_num real (c_double ) , intent(in) :: single_en_distance(nucl_num,walk_num) real (c_double ) , intent(in) :: coord(3) real (c_double ) , intent(in) :: nucl_coord(nucl_num,3) real (c_double ) , intent(out) :: en_rescaled_single_gl(4,nucl_num,walk_num) integer(c_int32_t), external :: qmckl_compute_en_rescaled_single_gl_doc_f info = qmckl_compute_en_rescaled_single_gl_doc_f & (context, & nucl_num, & type_nucl_num, & type_nucl_vector, & rescale_factor_en, & walk_num, & single_en_distance, & coord, & nucl_coord, & en_rescaled_single_gl) end function qmckl_compute_en_rescaled_single_gl_doc #+end_src ** test #+begin_src c :tangle (eval c_test) /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); rc = qmckl_set_point(context, 'N', elec_num, elec_coord, walk_num*elec_num*3); assert(rc == QMCKL_SUCCESS); double en_rescaled_gl[walk_num][nucl_num][elec_num][4]; rc = qmckl_get_electron_en_distance_rescaled_gl(context, &en_rescaled_gl[0][0][0][0]); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_electron_coord(context, 'N', &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); rc = qmckl_set_single_point(context, 'N', 2, new_coords, 3); assert (rc == QMCKL_SUCCESS); double single_en_rescaled_gl[walk_num][nucl_num][4]; rc = qmckl_get_en_rescaled_single_gl(context, &single_en_rescaled_gl[0][0][0]); assert (rc == QMCKL_SUCCESS); coords[0][2][0] = new_coords[0]; coords[0][2][1] = new_coords[1]; coords[0][2][2] = new_coords[2]; rc = qmckl_set_point(context, 'N', elec_num, &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_electron_en_distance_rescaled_gl(context, &en_rescaled_gl[0][0][0][0]); assert (rc == QMCKL_SUCCESS); for (int nw = 0; nw < walk_num; nw++) { for (int a = 0; a < nucl_num; a++) { for (int m = 0; m < 4; m++) { assert(fabs(en_rescaled_gl[nw][a][2][m] - single_en_rescaled_gl[nw][a][m]) < 1.e-12); } } } #+end_src * Delta ee gl ** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_jastrow_champ_single_ee_gl(qmckl_context context, double* const delta_ee_gl, const int64_t size_max); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_jastrow_champ_single_ee_gl(qmckl_context context, double* const delta_ee_gl, const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_exit_code rc; rc = qmckl_provide_jastrow_champ_single_ee_gl(context); if (rc != QMCKL_SUCCESS) return rc; qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int64_t sze = ctx->electron.walker.num * 4 * ctx->electron.num; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_jastrow_champ_single_ee_gl", "Array too small. Expected 4*walk_num*elec_num"); } memcpy(delta_ee_gl, ctx->single_point.delta_ee_gl, sze * sizeof(double)); return QMCKL_SUCCESS; } #+end_src #+begin_src f90 :tangle (eval fh_func) :comments org interface integer(qmckl_exit_code) function qmckl_get_jastrow_champ_single_ee_gl (context, & delta_ee_gl, size_max) bind(C) use, intrinsic :: iso_c_binding import implicit none integer (qmckl_context) , intent(in), value :: context integer(c_int64_t), intent(in), value :: size_max real(c_double), intent(out) :: delta_ee_gl(size_max) end function qmckl_get_jastrow_champ_single_ee_gl end interface #+end_src ** Provide :noexport: #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_provide_jastrow_champ_single_ee_gl(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_provide_jastrow_champ_single_ee_gl(qmckl_context context) { qmckl_exit_code rc; if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return qmckl_failwith( context, QMCKL_INVALID_CONTEXT, "qmckl_provide_jastrow_champ_single_ee_gl", NULL); } qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); if (!ctx->jastrow_champ.provided) { return qmckl_failwith( context, QMCKL_NOT_PROVIDED, "qmckl_provide_jastrow_champ_single_ee_gl", NULL); } /* Check if ee rescaled distance is provided */ rc = qmckl_provide_ee_distance_rescaled(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if ee rescaled distance deriv e is provided */ rc = qmckl_provide_ee_distance_rescaled_gl(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_ee_rescaled_single(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_ee_rescaled_single_gl(context); if(rc != QMCKL_SUCCESS) return rc; /* Compute if necessary */ if (ctx->single_point.date > ctx->single_point.delta_ee_gl_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { if (ctx->single_point.delta_ee_gl != NULL) { rc = qmckl_free(context, ctx->single_point.delta_ee_gl); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_jastrow_champ_single_ee_gl", "Unable to free ctx->single_point.delta_ee_gl"); } ctx->single_point.delta_ee_gl = NULL; } } /* Allocate array */ if (ctx->single_point.delta_ee_gl == NULL) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = ctx->electron.walker.num * 4 * ctx->electron.num * sizeof(double); double* delta_ee_gl = (double*) qmckl_malloc(context, mem_info); if (delta_ee_gl == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_jastrow_champ_single_ee_gl", NULL); } ctx->single_point.delta_ee_gl = delta_ee_gl; } rc = qmckl_compute_jastrow_champ_single_ee_gl(context, ctx->single_point.num, ctx->electron.walker.num, ctx->electron.num, ctx->electron.up_num, ctx->jastrow_champ.bord_num, ctx->jastrow_champ.b_vector, ctx->jastrow_champ.ee_distance_rescaled, ctx->jastrow_champ.ee_distance_rescaled_gl, ctx->single_point.ee_rescaled_single, ctx->single_point.ee_rescaled_single_gl, ctx->jastrow_champ.spin_independent, ctx->single_point.delta_ee_gl); if (rc != QMCKL_SUCCESS) { return rc; } ctx->single_point.delta_ee_gl_date = ctx->single_point.date; } return QMCKL_SUCCESS; } #+end_src ** Compute :PROPERTIES: :Name: qmckl_compute_jastrow_champ_single_ee_gl :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_single_ee_gl_args | Variable | Type | In/Out | Description | |---------------------------+-------------------------------------------+--------+-----------------------------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~num~ | ~int64_t~ | in | Number of walkers | | ~walk_num~ | ~int64_t~ | in | Number of walkers | | ~elec_num~ | ~int64_t~ | in | Number of electrons | | ~up_num~ | ~int64_t~ | in | Number of alpha electrons | | ~bord_num~ | ~int64_t~ | in | Number of coefficients | | ~b_vector~ | ~double[bord_num+1]~ | in | List of coefficients | | ~ee_distance_rescaled~ | ~double[walk_num][elec_num][elec_num]~ | in | Electron-electron distances | | ~ee_distance_rescaled_gl~ | ~double[walk_num][4][elec_num][elec_num]~ | in | Electron-electron distances | | ~ee_rescaled_single~ | ~double[walk_num][elec_num]~ | in | Electron-electron distances | | ~ee_rescaled_single_gl~ | ~double[walk_num][4][elec_num]~ | in | Electron-electron distances | | ~spin_independent~ | ~int32_t~ | in | If 1, same parameters for parallel and antiparallel spins | | ~delta_ee_gl~ | ~double[walk_num][elec_num][4]~ | out | Electron-electron distances | #+begin_src f90 :comments org :tangle (eval f) :noweb yes function qmckl_compute_jastrow_champ_single_ee_gl_doc( & context, num_in, walk_num, elec_num, up_num, bord_num, & b_vector, ee_distance_rescaled, ee_distance_rescaled_gl, & ee_rescaled_single, ee_rescaled_single_gl, & spin_independent, delta_ee_gl) & bind(C) result(info) use qmckl implicit none integer (qmckl_context), intent(in), value :: context integer (c_int64_t) , intent(in) , value :: num_in 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 :: up_num integer (c_int64_t) , intent(in) , value :: bord_num real (c_double ) , intent(in) :: b_vector(bord_num+1) real (c_double ) , intent(in) :: ee_distance_rescaled(elec_num,elec_num,walk_num) real (c_double ) , intent(in) :: ee_distance_rescaled_gl(4,elec_num,elec_num,walk_num) real (c_double ) , intent(in) :: ee_rescaled_single(elec_num,walk_num) real (c_double ) , intent(in) :: ee_rescaled_single_gl(4,elec_num,walk_num) integer (c_int32_t) , intent(in) , value :: spin_independent real (c_double ) , intent(out) :: delta_ee_gl(4,elec_num,walk_num) integer(qmckl_exit_code) :: info integer*8 :: i, j, k, nw, ii, num double precision :: x, x1, kf, x_old, x1_old double precision :: denom, invdenom, invdenom2, f double precision :: denom_old, invdenom_old, invdenom2_old, f_old double precision :: grad_c2, grad_c2_old double precision :: dx(4), dx_old(4) num = num_in + 1 info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif if (walk_num <= 0) then info = QMCKL_INVALID_ARG_3 return endif if (elec_num <= 0) then info = QMCKL_INVALID_ARG_4 return endif if (bord_num < 0) then info = QMCKL_INVALID_ARG_5 return endif if ((spin_independent < 0).or.(spin_independent > 1)) then info = QMCKL_INVALID_ARG_8 return endif do nw =1, walk_num delta_ee_gl(:,:,nw) = 0.0d0 do i = 1, elec_num if (i == num) cycle x = ee_rescaled_single(i,nw) x_old = ee_distance_rescaled(i,num,nw) denom = 1.0d0 + b_vector(2) * x invdenom = 1.0d0 / denom invdenom2 = invdenom * invdenom denom_old = 1.0d0 + b_vector(2) * x_old invdenom_old = 1.0d0 / denom_old invdenom2_old = invdenom_old * invdenom_old dx(1) = ee_rescaled_single_gl(1, i, nw) dx(2) = ee_rescaled_single_gl(2, i, nw) dx(3) = ee_rescaled_single_gl(3, i, nw) dx(4) = ee_rescaled_single_gl(4, i, nw) dx_old(1) = ee_distance_rescaled_gl(1, i, num, nw) dx_old(2) = ee_distance_rescaled_gl(2, i, num, nw) dx_old(3) = ee_distance_rescaled_gl(3, i, num, nw) dx_old(4) = ee_distance_rescaled_gl(4, i, num, nw) grad_c2 = dx(1)*dx(1) + dx(2)*dx(2) + dx(3)*dx(3) grad_c2_old = dx_old(1)*dx_old(1) + dx_old(2)*dx_old(2) + dx_old(3)*dx_old(3) if (spin_independent == 1) then f = b_vector(1) * invdenom2 f_old = b_vector(1) * invdenom2_old else if((i <= up_num .and. num <= up_num ) .or. (i > up_num .and. num > up_num)) then f = 0.5d0 * b_vector(1) * invdenom2 f_old = 0.5d0 * b_vector(1) * invdenom2_old else f = b_vector(1) * invdenom2 f_old = b_vector(1) * invdenom2_old end if end if delta_ee_gl(1,i,nw) = delta_ee_gl(1,i,nw) + f * dx(1) - f_old * dx_old(1) delta_ee_gl(2,i,nw) = delta_ee_gl(2,i,nw) + f * dx(2) - f_old * dx_old(2) delta_ee_gl(3,i,nw) = delta_ee_gl(3,i,nw) + f * dx(3) - f_old * dx_old(3) delta_ee_gl(4,i,nw) = delta_ee_gl(4,i,nw) & + f * (dx(4) - 2.d0 * b_vector(2) * grad_c2 * invdenom) & - f_old * (dx_old(4) - 2.d0 * b_vector(2) * grad_c2_old * invdenom_old) kf = 2.d0 x1 = x x1_old = x_old x = 1.d0 x_old = 1.d0 do k=2, bord_num f = b_vector(k+1) * kf * x f_old = b_vector(k+1) * kf * x_old delta_ee_gl(1,i,nw) = delta_ee_gl(1,i,nw) + f * x1 * dx(1) - f_old * x1_old * dx_old(1) delta_ee_gl(2,i,nw) = delta_ee_gl(2,i,nw) + f * x1 * dx(2) - f_old * x1_old * dx_old(2) delta_ee_gl(3,i,nw) = delta_ee_gl(3,i,nw) + f * x1 * dx(3) - f_old * x1_old * dx_old(3) delta_ee_gl(4,i,nw) = delta_ee_gl(4,i,nw) & + f * (x1 * dx(4) + (kf-1.d0) * grad_c2) & - f_old * (x1_old * dx_old(4) + (kf-1.d0) * grad_c2_old) x = x*x1 x_old = x_old*x1_old kf = kf + 1.d0 end do end do end do end function qmckl_compute_jastrow_champ_single_ee_gl_doc #+end_src # #+CALL: generate_c_header(table=qmckl_factor_ee_gl_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_jastrow_champ_single_ee_gl (const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t up_num, const int64_t bord_num, const double* b_vector, const double* ee_distance_rescaled, const double* ee_distance_rescaled_gl, const double* ee_rescaled_single, const double* ee_rescaled_single_gl, const int32_t spin_independent, double* const delta_ee_gl ); #+end_src #+begin_src c :tangle (eval h_private_func) :comments org qmckl_exit_code qmckl_compute_jastrow_champ_single_ee_gl_doc (const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t up_num, const int64_t bord_num, const double* b_vector, const double* ee_distance_rescaled, const double* ee_distance_rescaled_gl, const double* ee_rescaled_single, const double* ee_rescaled_single_gl, const int32_t spin_independent, double* const delta_ee_gl ); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes qmckl_exit_code qmckl_compute_jastrow_champ_single_ee_gl (const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t up_num, const int64_t bord_num, const double* b_vector, const double* ee_distance_rescaled, const double* ee_distance_rescaled_gl, const double* ee_rescaled_single, const double* ee_rescaled_single_gl, const int32_t spin_independent, double* const delta_ee_gl ) { #ifdef HAVE_HPC return qmckl_compute_jastrow_champ_single_ee_gl_doc #else return qmckl_compute_jastrow_champ_single_ee_gl_doc #endif (context, num, walk_num, elec_num, up_num, bord_num, b_vector, ee_distance_rescaled, ee_distance_rescaled_gl, ee_rescaled_single, ee_rescaled_single_gl, spin_independent, delta_ee_gl ); } #+end_src ** test #+begin_src c :tangle (eval c_test) /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); rc = qmckl_set_point(context, 'N', elec_num, elec_coord, walk_num*elec_num*3); assert(rc == QMCKL_SUCCESS); double ee_gl_old[walk_num][4][elec_num]; rc = qmckl_get_jastrow_champ_factor_ee_gl(context, &ee_gl_old[0][0][0], walk_num*elec_num*4); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_electron_coord(context, 'N', &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); rc = qmckl_set_single_point(context, 'N', 2, new_coords, 3); assert (rc == QMCKL_SUCCESS); double delta_ee_gl[walk_num][elec_num][4]; rc = qmckl_get_jastrow_champ_single_ee_gl(context, &delta_ee_gl[0][0][0], walk_num*elec_num*4); assert (rc == QMCKL_SUCCESS); coords[0][2][0] = new_coords[0]; coords[0][2][1] = new_coords[1]; coords[0][2][2] = new_coords[2]; rc = qmckl_set_point(context, 'N', elec_num, &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); double ee_gl_new[walk_num][4][elec_num]; rc = qmckl_get_jastrow_champ_factor_ee_gl(context, &ee_gl_new[0][0][0], walk_num*elec_num*4); assert (rc == QMCKL_SUCCESS); for (int nw = 0; nw < walk_num; nw++) { for (int i = 0; i < elec_num; i++) { for (int m = 0; m < 4; m++) { if (i == 2) continue; //printf("%f\n",(ee_gl_new[nw][m][i] - ee_gl_old[nw][m][i])); //printf("%f\n",delta_ee_gl[nw][i][m]); assert(fabs((ee_gl_new[nw][m][i] - ee_gl_old[nw][m][i]) - delta_ee_gl[nw][i][m]) < 1.e-12); } } } #+end_src * Delta en gl ** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_jastrow_champ_single_en_gl(qmckl_context context, double* const delta_en_gl, const int64_t size_max); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_jastrow_champ_single_en_gl(qmckl_context context, double* const delta_en_gl, const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_exit_code rc; rc = qmckl_provide_jastrow_champ_single_en_gl(context); if (rc != QMCKL_SUCCESS) return rc; qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int64_t sze = ctx->electron.walker.num * 4 * ctx->electron.num; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_jastrow_champ_single_en_gl", "Array too small. Expected 4*walker.num*elec_num"); } memcpy(delta_en_gl, ctx->single_point.delta_en_gl, sze*sizeof(double)); return QMCKL_SUCCESS; } #+end_src #+begin_src f90 :tangle (eval fh_func) :comments org interface integer(qmckl_exit_code) function qmckl_get_jastrow_champ_single_en_gl (context, & delta_en_gl, size_max) bind(C) use, intrinsic :: iso_c_binding import implicit none integer (qmckl_context) , intent(in), value :: context integer(c_int64_t), intent(in), value :: size_max real(c_double), intent(out) :: delta_en_gl(size_max) end function qmckl_get_jastrow_champ_single_en_gl end interface #+end_src ** Provide :noexport: #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_provide_jastrow_champ_single_en_gl(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_provide_jastrow_champ_single_en_gl(qmckl_context context) { qmckl_exit_code rc; if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return qmckl_failwith( context, QMCKL_INVALID_CONTEXT, "qmckl_provide_jastrow_champ_single_en_gl", NULL); } qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); if (!ctx->jastrow_champ.provided) { return qmckl_failwith( context, QMCKL_NOT_PROVIDED, "qmckl_provide_jastrow_champ_single_en_gl", NULL); } /* Check if en rescaled distance is provided */ rc = qmckl_provide_en_distance_rescaled(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if en rescaled distance derivatives is provided */ rc = qmckl_provide_en_distance_rescaled_gl(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if en rescaled distance is provided */ rc = qmckl_provide_en_rescaled_single(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if en rescaled distance derivatives is provided */ rc = qmckl_provide_en_rescaled_single_gl(context); if(rc != QMCKL_SUCCESS) return rc; /* Compute if necessary */ if (ctx->single_point.date > ctx->single_point.delta_en_gl_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { if (ctx->single_point.delta_en_gl != NULL) { rc = qmckl_free(context, ctx->single_point.delta_en_gl); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_jastrow_champ_single_en_gl", "Unable to free ctx->single_point.delta_en_gl"); } ctx->single_point.delta_en_gl = NULL; } } /* Allocate array */ if (ctx->single_point.delta_en_gl == NULL) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = ctx->electron.walker.num * 4 * ctx->electron.num * sizeof(double); double* delta_en_gl = (double*) qmckl_malloc(context, mem_info); if (delta_en_gl == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_jastrow_champ_single_en_gl", NULL); } ctx->single_point.delta_en_gl = delta_en_gl; } rc = qmckl_compute_jastrow_champ_single_en_gl(context, ctx->single_point.num, ctx->electron.walker.num, ctx->electron.num, ctx->nucleus.num, ctx->jastrow_champ.type_nucl_num, ctx->jastrow_champ.type_nucl_vector, ctx->jastrow_champ.aord_num, ctx->jastrow_champ.a_vector, ctx->jastrow_champ.en_distance_rescaled, ctx->jastrow_champ.en_distance_rescaled_gl, ctx->single_point.en_rescaled_single, ctx->single_point.en_rescaled_single_gl, ctx->single_point.delta_en_gl); if (rc != QMCKL_SUCCESS) { return rc; } ctx->single_point.delta_en_gl_date = ctx->single_point.date; } return QMCKL_SUCCESS; } #+end_src ** Compute :PROPERTIES: :Name: qmckl_compute_jastrow_champ_single_en_gl :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_single_en_gl_args | Variable | Type | In/Out | Description | |---------------------------+-------------------------------------------+--------+---------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~wnum~ | ~int64_t~ | in | Number of walkers | | ~walk_num~ | ~int64_t~ | in | Number of walkers | | ~elec_num~ | ~int64_t~ | in | Number of electrons | | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | | ~type_nucl_num~ | ~int64_t~ | in | Number of unique nuclei | | ~type_nucl_vector~ | ~int64_t[nucl_num]~ | in | IDs of unique nuclei | | ~aord_num~ | ~int64_t~ | in | Number of coefficients | | ~a_vector~ | ~double[type_nucl_num][aord_num+1]~ | in | List of coefficients | | ~en_distance_rescaled~ | ~double[walk_num][nucl_num][elec_num]~ | in | Electron-nucleus distances | | ~en_distance_rescaled_gl~ | ~double[walk_num][nucl_num][elec_num][4]~ | in | Electron-nucleus distance derivatives | | ~en_rescaled_single~ | ~double[walk_num][nucl_num]~ | in | Electron-nucleus distances | | ~en_rescaled_single_gl~ | ~double[walk_num][nucl_num][4]~ | in | Electron-nucleus distance derivatives | | ~delta_en_gl~ | ~double[walk_num][elec_num][4]~ | out | Electron-nucleus jastrow | #+begin_src f90 :comments org :tangle (eval f) :noweb yes function qmckl_compute_jastrow_champ_single_en_gl_doc( & context, num_in, walk_num, elec_num, nucl_num, type_nucl_num, & type_nucl_vector, aord_num, a_vector, & en_distance_rescaled, en_distance_rescaled_gl, en_rescaled_single, en_rescaled_single_gl, delta_en_gl) & bind(C) result(info) use qmckl implicit none integer (qmckl_context), intent(in), value :: context integer (c_int64_t) , intent(in) , value :: num_in 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 :: type_nucl_num integer (c_int64_t) , intent(in) :: type_nucl_vector(nucl_num) integer (c_int64_t) , intent(in) , value :: aord_num real (c_double ) , intent(in) :: a_vector(aord_num+1,type_nucl_num) real (c_double ) , intent(in) :: en_distance_rescaled(elec_num,nucl_num,walk_num) real (c_double ) , intent(in) :: en_distance_rescaled_gl(4, elec_num,nucl_num,walk_num) real (c_double ) , intent(in) :: en_rescaled_single(nucl_num,walk_num) real (c_double ) , intent(in) :: en_rescaled_single_gl(4, nucl_num,walk_num) real (c_double ) , intent(out) :: delta_en_gl(4,elec_num,walk_num) integer(qmckl_exit_code) :: info integer*8 :: i, a, k, nw, ii, num double precision :: x, x1, kf, x_old, x1_old double precision :: denom, invdenom, invdenom2, f double precision :: denom_old, invdenom_old, invdenom2_old, f_old double precision :: grad_c2, grad_c2_old double precision :: dx(4), dx_old(4) num = num_in + 1 info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif if (walk_num <= 0) then info = QMCKL_INVALID_ARG_3 return endif if (elec_num <= 0) then info = QMCKL_INVALID_ARG_4 return endif if (nucl_num <= 0) then info = QMCKL_INVALID_ARG_5 return endif if (aord_num < 0) then info = QMCKL_INVALID_ARG_8 return endif do nw =1, walk_num delta_en_gl(:,:,nw) = 0.0d0 do a = 1, nucl_num x_old = en_distance_rescaled(num,a,nw) x = en_rescaled_single(a,nw) denom = 1.0d0 + a_vector(2, type_nucl_vector(a)+1) * x invdenom = 1.0d0 / denom invdenom2 = invdenom*invdenom denom_old = 1.0d0 + a_vector(2, type_nucl_vector(a)+1) * x_old invdenom_old = 1.0d0 / denom_old invdenom2_old = invdenom_old*invdenom_old dx(1) = en_rescaled_single_gl(1,a,nw) dx(2) = en_rescaled_single_gl(2,a,nw) dx(3) = en_rescaled_single_gl(3,a,nw) dx(4) = en_rescaled_single_gl(4,a,nw) dx_old(1) = en_distance_rescaled_gl(1,num,a,nw) dx_old(2) = en_distance_rescaled_gl(2,num,a,nw) dx_old(3) = en_distance_rescaled_gl(3,num,a,nw) dx_old(4) = en_distance_rescaled_gl(4,num,a,nw) f = a_vector(1, type_nucl_vector(a)+1) * invdenom2 grad_c2 = dx(1)*dx(1) + dx(2)*dx(2) + dx(3)*dx(3) f_old = a_vector(1, type_nucl_vector(a)+1) * invdenom2_old grad_c2_old = dx_old(1)*dx_old(1) + dx_old(2)*dx_old(2) + dx_old(3)*dx_old(3) delta_en_gl(1,num,nw) = delta_en_gl(1,num,nw) + f * dx(1) - f_old * dx_old(1) delta_en_gl(2,num,nw) = delta_en_gl(2,num,nw) + f * dx(2) - f_old * dx_old(2) delta_en_gl(3,num,nw) = delta_en_gl(3,num,nw) + f * dx(3) - f_old * dx_old(3) delta_en_gl(4,num,nw) = delta_en_gl(4,num,nw) & + f * (dx(4) - 2.d0 * a_vector(2, type_nucl_vector(a)+1) * grad_c2 * invdenom) & - f_old * (dx_old(4) - 2.d0 * a_vector(2, type_nucl_vector(a)+1) * grad_c2_old * invdenom_old) kf = 2.d0 x1 = x x = 1.d0 x1_old = x_old x_old = 1.d0 do k=2, aord_num f = a_vector(k+1,type_nucl_vector(a)+1) * kf * x f_old = a_vector(k+1,type_nucl_vector(a)+1) * kf * x_old delta_en_gl(1,num,nw) = delta_en_gl(1,num,nw) + f * x1 * dx(1) - f_old * x1_old * dx_old(1) delta_en_gl(2,num,nw) = delta_en_gl(2,num,nw) + f * x1 * dx(2) - f_old * x1_old * dx_old(2) delta_en_gl(3,num,nw) = delta_en_gl(3,num,nw) + f * x1 * dx(3) - f_old * x1_old * dx_old(3) delta_en_gl(4,num,nw) = delta_en_gl(4,num,nw) & + f * (x1 * dx(4) + (kf-1.d0) * grad_c2) & - f_old * (x1_old * dx_old(4) + (kf-1.d0) * grad_c2_old) x = x*x1 x_old = x_old*x1_old kf = kf + 1.d0 end do end do end do end function qmckl_compute_jastrow_champ_single_en_gl_doc #+end_src # #+CALL: generate_c_header(table=qmckl_factor_en_gl_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_jastrow_champ_single_en_gl_doc ( const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t type_nucl_num, const int64_t* type_nucl_vector, const int64_t aord_num, const double* a_vector, const double* en_distance_rescaled, const double* en_distance_rescaled_gl, const double* en_rescaled_single, const double* en_rescaled_single_gl, double* const delta_en_gl ); qmckl_exit_code qmckl_compute_jastrow_champ_single_en_gl ( const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t type_nucl_num, const int64_t* type_nucl_vector, const int64_t aord_num, const double* a_vector, const double* en_distance_rescaled, const double* en_distance_rescaled_gl, const double* en_rescaled_single, const double* en_rescaled_single_gl, double* const delta_en_gl ); #+end_src #+begin_src c :tangle (eval c) :comments org :exports none qmckl_exit_code qmckl_compute_jastrow_champ_single_en_gl (const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t type_nucl_num, const int64_t* type_nucl_vector, const int64_t aord_num, const double* a_vector, const double* en_distance_rescaled, const double* en_distance_rescaled_gl, const double* en_rescaled_single, const double* en_rescaled_single_gl, double* const delta_en_gl ) { #ifdef HAVE_HPC return qmckl_compute_jastrow_champ_single_en_gl_doc #else return qmckl_compute_jastrow_champ_single_en_gl_doc #endif (context, num, walk_num, elec_num, nucl_num, type_nucl_num, type_nucl_vector, aord_num, a_vector, en_distance_rescaled, en_distance_rescaled_gl, en_rescaled_single, en_rescaled_single_gl, delta_en_gl ); } #+end_src ** test #+begin_src c :tangle (eval c_test) /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); rc = qmckl_set_point(context, 'N', elec_num, elec_coord, walk_num*elec_num*3); assert(rc == QMCKL_SUCCESS); double en_gl_old[walk_num][4][elec_num]; rc = qmckl_get_jastrow_champ_factor_en_gl(context, &en_gl_old[0][0][0], walk_num*elec_num*4); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_electron_coord(context, 'N', &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); rc = qmckl_set_single_point(context, 'N', 2, new_coords, 3); assert (rc == QMCKL_SUCCESS); double delta_en_gl[walk_num][elec_num][4]; rc = qmckl_get_jastrow_champ_single_en_gl(context, &delta_en_gl[0][0][0], walk_num*elec_num*4); assert (rc == QMCKL_SUCCESS); coords[0][2][0] = new_coords[0]; coords[0][2][1] = new_coords[1]; coords[0][2][2] = new_coords[2]; rc = qmckl_set_point(context, 'N', elec_num, &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); double en_gl_new[walk_num][4][elec_num]; rc = qmckl_get_jastrow_champ_factor_en_gl(context, &en_gl_new[0][0][0], walk_num*elec_num*4); assert (rc == QMCKL_SUCCESS); for (int nw = 0; nw < walk_num; nw++) { for (int i = 0; i < elec_num; i++) { for (int m = 0; m < 4; m++) { assert(fabs((en_gl_new[nw][m][i] - en_gl_old[nw][m][i]) - delta_en_gl[nw][i][m]) < 1.e-12); } } } #+end_src * Accept single electron move #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_jastrow_champ_single_accept(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_jastrow_champ_single_accept(qmckl_context context) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_exit_code rc; rc = qmckl_provide_jastrow_champ_single_en_gl(context); if (rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_jastrow_champ_single_ee_gl(context); if (rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_jastrow_champ_single_een_gl(context); if (rc != QMCKL_SUCCESS) return rc; qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); /* for (nw = 0; nw < ctx->electron.walker.num; nw++) { ctx->jastrow_champ.factor_ee[nw] = ctx->jastrow_champ.factor_ee[nw] + ctx->single_point.delta_ee[nw]; ctx->jastrow_champ.factor_en[nw] = ctx->jastrow_champ.factor_en[nw] + ctx->single_point.delta_en[nw]; ctx->jastrow_champ.factor_een[nw] = ctx->jastrow_champ.factor_een[nw] + ctx->single_point.delta_een[nw]; for (l = 0; l <= ctx->jastrow_champ.cord_num; l++){ for (a = 0; a < ctx->nucleus.num; a++){ for (i = 0; i < ctx->electron.num; i++) { for (m = 0; m < ctx->jastrow_champ.cord_num; m++){ ctx->jastrow_champ.tmp_c[nw][m][l][a][i] = ctx->jastrow_champ.tmp_c[nw][m][l][a][i] + ctx->single_point.delta_p[nw][m][l][a][i]; } } ctx->jastrow_champ.een_rescaled_n[nw][l][a][ctx->single_point.num] = ctx->single_point.een_rescaled_single_n[nw][l][a][ctx->single_point.num] + ctx->single_point.een_rescaled_delta_n[nw][l][a]; } for (i = 0; i < ctx->electron.num; i++) { ctx->jastrow_champ.een_rescaled_e[nw][l][i][ctx->single_point.num] = ctx->single_point.een_rescaled_single_e[nw][l][i][ctx->single_point.num] + ctx->single_point.een_rescaled_delta_e[nw][l][i]; ctx->jastrow_champ.een_rescaled_e[nw][l][ctx->single_point.num][i] = ctx->single_point.een_rescaled_single_e[nw][l][ctx->single_point.num][i] + ctx->single_point.een_rescaled_delta_e[nw][l][i]; } } for (k = 0; k < 4; k++) { for (i = 0; i < ctx->electron.num; i++) { ctx->jastrow_champ.factor_ee_gl[nw][k][i] = ctx->jastrow_champ.factor_ee_gl[nw][k][i] + ctx->single_point.delta_ee_gl[nw][k][i]; ctx->jastrow_champ.factor_en_gl[nw][k][i] = ctx->jastrow_champ.factor_en_gl[nw][k][i] + ctx->single_point.delta_en_gl[nw][k][i]; ctx->jastrow_champ.factor_een_gl[nw][k][i] = ctx->jastrow_champ.factor_een_gl[nw][k][i] + ctx->single_point.delta_een_gl[nw][k][i]; ctx->jastrow_champ.ee_distance_rescaled_gl[nw][i][ctx->single_point.num][k] = ctx->single_point.ee_rescaled_single_gl[nw][i][k]; ctx->jastrow_champ.ee_distance_rescaled_gl[nw][ctx->single_point.num][i][k] = ctx->single_point.ee_rescaled_single_gl[nw][i][k]; for (l = 0; l <= ctx->jastrow_champ.cord_num; l++){ ctx->jastrow_champ.een_rescaled_e_gl[nw][l][i][k][ctx->single_point.num] = ctx->jastrow_champ.een_rescaled_single_e_gl[nw][l][i][k]; ctx->jastrow_champ.een_rescaled_e_gl[nw][l][ctx->single_point.num][k][i] = ctx->jastrow_champ.een_rescaled_single_e_gl[nw][l][i][k]; for (m = 0; m < ctx->jastrow_champ.cord_num; m++){ for (a = 0; a < ctx->nucleus.num; a++){ ctx->jastrow_champ.dtmp_c[nw][i][k][a][l][m] = ctx->jastrow_champ.dtmp_c[nw][i][k][a][l][m] + ctx->single_point.delta_p_gl[nw][m][l][a][i][k] } } } for (a = 0; a < ctx->nucleus.num; a++) { ctx->jastrow_champ.en_distance_rescaled_gl[nw][a][ctx->single_point.num][k] = ctx->single_point.en_rescaled_single_gl[nw][a][k]; for (l = 0; l <= ctx->jastrow_champ.cord_num; l++) { ctx->jastorw_champ.een_rescaled_n_gl[nw][l][a][k][ctx->single_point.num] = ctx->single_point.een_rescaled_single_n_gl[nw][l][a][k]; } } } } for (a = 0; a < ctx->nucleus.num; a++) { ctx->jastrow_champ.en_distance_rescaled[nw][a][ctx->single_point.num] = ctx->single_point.en_rescaled_single[nw][a]; for (i = 0; i < ctx->electron.num; i++) { ctx->jastrow_champ.ee_distance_rescaled[nw][i][ctx->single_point.num] = ctx->single_point.ee_rescaled_single[nw][i]; ctx->jastrow_champ.ee_distance_rescaled[nw][ctx->single_point.num][i] = ctx->single_point.ee_rescaled_single[nw][i]; } } } ctx->jastrow_champ.ee_distance_rescaled_date = ctx->single_point.date; ctx->jastrow_champ.ee_distance_rescaled_gl_date = ctx->single_point.date; ctx->jastrow_champ.en_distance_rescaled_date = ctx->single_point.date; ctx->jastrow_champ.en_distance_rescaled_gl_date = ctx->single_point.date; ctx->jastrow_champ.een_rescaled_n_date = ctx->single_point.date; ctx->jastrow_champ.een_rescaled_e_date = ctx->single_point.date; ctx->jastrow_champ.een_rescaled_e_gl_date = ctx->single_point.date; ctx->jastrow_champ.een_rescaled_n_gl_date = ctx->single_point.date; ctx->jastrow_champ.factor_ee_date = ctx->single_point.date; ctx->jastrow_champ.factor_en_date = ctx->single_point.date; ctx->jastrow_champ.factor_een_date = ctx->single_point.date; ctx->jastrow_champ.factor_ee_gl_date = ctx->single_point.date; ctx->jastrow_champ.factor_en_gl_date = ctx->single_point.date; ctx->jastrow_champ.factor_een_gl_date = ctx->single_point.date; ctx->jastrow_champ.dtmp_c_date = ctx->single_point.date; ctx->jastrow_champ.tmp_c_date = ctx->single_point.date; for (k = 0; k < 3; k++) { ctx->point.coord[k] = ctx->single_point.coord[k]; } ctx->point.date = ctx->single_point.date; ctx->electron.walker.point.date = ctx->single_point.date; ctx->date = ctx->single_point.date; */ return QMCKL_SUCCESS; } #+end_src #+begin_src f90 :tangle (eval fh_func) :comments org interface integer(qmckl_exit_code) function qmckl_get_jastrow_champ_single_accept (context) bind(C) use, intrinsic :: iso_c_binding import implicit none integer (qmckl_context) , intent(in), value :: context end function qmckl_get_jastrow_champ_single_accept end interface #+end_src * End of files :noexport: #+begin_src c :tangle (eval h_private_type) #endif #+end_src #+begin_src c :tangle (eval h_private_func) #endif #+end_src ** Test #+begin_src c :tangle (eval c_test) rc = qmckl_context_destroy(context); assert (rc == QMCKL_SUCCESS); return 0; } #+end_src