#+TITLE: CHAMP Jastrow Factor Quad #+SETUPFILE: ../tools/theme.setup #+INCLUDE: ../tools/lib.org todo --optimize een_gl --optimized delta_p_gl * 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_QUAD_HPF #define QMCKL_JASTROW_CHAMP_QUAD_HPF #+end_src #+begin_src c :tangle (eval h_private_type) #ifndef QMCKL_JASTROW_CHAMP_QUAD_HPT #define QMCKL_JASTROW_CHAMP_QUAD_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" #include "qmckl_jastrow_champ_quad_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" #include "qmckl_jastrow_champ_quad_private_type.h" #include "qmckl_jastrow_champ_quad_private_func.h" #+end_src * Context ** Data structure #+begin_src c :comments org :tangle (eval h_private_type) typedef struct qmckl_jastrow_champ_quad_struct{ int64_t num; int64_t* indices; uint64_t date; qmckl_matrix coord; double * een_rescaled_quad_e; uint64_t een_rescaled_quad_e_date; uint64_t een_rescaled_quad_e_maxsize; double * een_rescaled_quad_n; uint64_t een_rescaled_quad_n_date; uint64_t een_rescaled_quad_n_maxsize; double* quad_ee_distance; uint64_t quad_ee_distance_date; uint64_t quad_ee_distance_maxsize; double* quad_en_distance; uint64_t quad_en_distance_date; uint64_t quad_en_distance_maxsize; double* delta_een; uint64_t delta_een_date; uint64_t delta_een_maxsize; double* delta_p; uint64_t delta_p_date; uint64_t delta_p_maxsize; double* ee_rescaled_quad; uint64_t ee_rescaled_quad_date; uint64_t ee_rescaled_quad_maxsize; double* en_rescaled_quad; uint64_t en_rescaled_quad_date; uint64_t en_rescaled_quad_maxsize; double* delta_en; uint64_t delta_en_date; uint64_t delta_en_maxsize; double* delta_ee; uint64_t delta_ee_date; uint64_t delta_ee_maxsize; double * een_rescaled_quad_e_gl; uint64_t een_rescaled_quad_e_gl_date; uint64_t een_rescaled_quad_e_gl_maxsize; double * een_rescaled_quad_n_gl; uint64_t een_rescaled_quad_n_gl_date; uint64_t een_rescaled_quad_n_gl_maxsize; double* delta_p_gl; uint64_t delta_p_gl_date; uint64_t delta_p_gl_maxsize; double* delta_een_gl; uint64_t delta_een_gl_date; uint64_t delta_een_gl_maxsize; double* delta_een_g; uint64_t delta_een_g_date; uint64_t delta_een_g_maxsize; double* ee_rescaled_quad_gl; uint64_t ee_rescaled_quad_gl_date; uint64_t ee_rescaled_quad_gl_maxsize; double* en_rescaled_quad_gl; uint64_t en_rescaled_quad_gl_date; uint64_t en_rescaled_quad_gl_maxsize; double* delta_en_gl; uint64_t delta_en_gl_date; uint64_t delta_en_gl_maxsize; double* delta_ee_gl; uint64_t delta_ee_gl_date; uint64_t delta_ee_gl_maxsize; double * forces_jastrow_quad_en; uint64_t forces_jastrow_quad_en_date; uint64_t forces_jastrow_quad_en_maxsize; double * forces_jastrow_quad_een; uint64_t forces_jastrow_quad_een_date; uint64_t forces_jastrow_quad_een_maxsize; double * forces_delta_p; uint64_t forces_delta_p_date; uint64_t forces_delta_p_maxsize; } qmckl_jastrow_champ_quad_struct; #+end_src ** Test :noexport: #+begin_src c :tangle (eval c_test) /* Reference input data */ int64_t walk_num = 1; 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 ; ielectron.walker.num > 1) { return qmckl_failwith( context, QMCKL_FAILURE, "qmckl_set_quad_points", "Not implemented for multiple walkers"); } if (size_max < 3*num) { return qmckl_failwith( context, QMCKL_INVALID_ARG_4, "qmckl_set_quad_points", "Array too small"); } qmckl_exit_code rc; if (ctx->quad_point.coord.data != NULL) { rc = qmckl_matrix_free(context, &(ctx->quad_point.coord)); assert (rc == QMCKL_SUCCESS); } ctx->quad_point.coord = qmckl_matrix_alloc(context, num, 3); if (ctx->quad_point.coord.data == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_set_quad_points", NULL); } if (ctx->quad_point.indices != NULL) { rc = qmckl_free(context, ctx->quad_point.indices); assert (rc == QMCKL_SUCCESS); } qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = num * sizeof(int64_t); ctx->quad_point.indices = qmckl_malloc(context, mem_info); if (ctx->quad_point.indices == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_set_quad_points", NULL); } ctx->quad_point.num = num; int64_t *b = ctx->quad_point.indices; for (int i = 0; i < num; i++){ b[i] = indices[i]; } if (transp == 'N') { double *a = ctx->quad_point.coord.data; for (int64_t i=0 ; i<3*num ; ++i) { a[i] = coord[i]; } } else { for (int64_t i=0 ; iquad_point.coord, i, 0) = coord[i*num + 0]; qmckl_mat(ctx->quad_point.coord, i, 1) = coord[i*num + 1]; qmckl_mat(ctx->quad_point.coord, i, 2) = coord[i*num + 2]; } } /* Increment the date of the quad point */ ctx->quad_point.date += 1UL; return QMCKL_SUCCESS; } #+end_src #+begin_src f90 :comments org :tangle (eval fh_func) :noweb yes interface integer(qmckl_exit_code) function qmckl_set_quad_points(context, & transp, num, indices, coord, size_max) bind(C) 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 integer (c_int64_t) , intent(in) :: indices(*) real (c_double ) , intent(in) :: coord(*) integer (c_int64_t) , intent(in) , value :: size_max end function end interface #+end_src ** Touch #+begin_src c :comments org :tangle (eval h_func) qmckl_exit_code qmckl_quad_touch (const qmckl_context context); #+end_src #+begin_src c :tangle (eval c) :exports none qmckl_exit_code qmckl_quad_touch(const qmckl_context context) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return qmckl_failwith( context, QMCKL_INVALID_CONTEXT, "qmckl_quad_touch", NULL); } qmckl_context_struct* const ctx = (qmckl_context_struct*) context; ctx->date += 1UL; ctx->point.date = ctx-> date; ctx->electron.walker.point.date = ctx-> date; ctx->quad_point.date = ctx-> date; return QMCKL_SUCCESS; } #+end_src * Electron-electron and electron-nucleus distances for quad point In order to calculate the $\delta J$, we need to have to updated distances for the quad electron. ** Electron-electron distances Electron-electron distance between the quad electron and all electrons for all walkers. Dimension is ~[walk_num][elec_num]~. *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_quad_electron_ee_distance(qmckl_context context, double* const distance, const int64_t size_max); #+end_src #+begin_src f90 :tangle (eval fh_func) :comments org :exports none interface integer(qmckl_exit_code) function qmckl_get_quad_electron_ee_distance(context, distance, size_max) & bind(C) use, intrinsic :: iso_c_binding import implicit none integer (c_int64_t) , intent(in) , value :: context real (c_double ) , intent(out) :: distance(*) integer (c_int64_t) , intent(in) , value :: size_max end function end interface #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_quad_electron_ee_distance(qmckl_context context, double* const distance, const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_exit_code rc; rc = qmckl_provide_quad_ee_distance(context); if (rc != QMCKL_SUCCESS) return rc; qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); if (distance == NULL) { return qmckl_failwith( context, QMCKL_INVALID_ARG_2, "qmckl_get_quad_electron_ee_distance", "distance is a NULL pointer"); } int64_t sze = ctx->electron.num * ctx->quad_point.num; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_quad_electron_ee_distance", "Array too small. Expected ctx->electron.num * ctx->quad_point.num"); } memcpy(distance, ctx->quad_point.quad_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_quad_ee_distance(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_provide_quad_ee_distance(qmckl_context context) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_exit_code rc; qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); /* Compute if necessary */ if (ctx->quad_point.date > ctx->quad_point.quad_ee_distance_date) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = ctx->electron.num * ctx->quad_point.num * sizeof(double); if (mem_info.size > ctx->quad_point.quad_ee_distance_maxsize) { if (ctx->quad_point.quad_ee_distance != NULL) { rc = qmckl_free(context, ctx->quad_point.quad_ee_distance); assert(rc == QMCKL_SUCCESS); ctx->quad_point.quad_ee_distance = NULL; } } /* Allocate array */ if (ctx->quad_point.quad_ee_distance == NULL) { double* quad_ee_distance = (double*) qmckl_malloc(context, mem_info); ctx->quad_point.quad_ee_distance_maxsize = mem_info.size; if (quad_ee_distance == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_quad_ee_distance", NULL); } ctx->quad_point.quad_ee_distance = quad_ee_distance; } rc = qmckl_compute_quad_ee_distance(context, ctx->quad_point.num, ctx->quad_point.indices, ctx->electron.num, ctx->electron.walker.num, ctx->electron.walker.point.coord.data, ctx->quad_point.coord.data, ctx->quad_point.quad_ee_distance); if (rc != QMCKL_SUCCESS) { return rc; } ctx->quad_point.quad_ee_distance_date = ctx->quad_point.date; } return QMCKL_SUCCESS; } #+end_src *** Compute :PROPERTIES: :Name: qmckl_compute_quad_ee_distance :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_quad_ee_distance_args |----------------------+---------------------------------+--------+-------------------------------------------------| | Variable | Type | In/Out | Description | |----------------------+---------------------------------+--------+-------------------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~num~ | ~int64_t~ | in | Index of quad electron | | ~indices~ | ~int64_t[num]~ | in | Indices of quad electron | | ~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 | | ~quad_coord~ | ~double[num][3]~ | in | quad electron coordinates | | ~quad_ee_distance~ | ~double[num][elec_num]~ | out | Electron-electron distances for quad electron | |----------------------+---------------------------------+--------+-------------------------------------------------| #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer(qmckl_exit_code) function qmckl_compute_quad_ee_distance(context, & num, indices, elec_num, walk_num, coord, quad_coord, quad_ee_distance) & result(info) bind(C) use, intrinsic :: iso_c_binding use qmckl implicit none integer(qmckl_context), intent(in) :: context integer (c_int64_t) , intent(in) , value :: elec_num, num integer (c_int64_t) , intent(in) :: indices(num) integer (c_int64_t) , intent(in) , value :: walk_num real (c_double ) , intent(in) :: coord(elec_num,walk_num,3) real (c_double ) , intent(in) :: quad_coord(3,num) real (c_double ) , intent(out) :: quad_ee_distance(elec_num,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,num info = qmckl_distance(context, 'T', 'N', elec_num, 1_8, & coord(1,1,1), elec_num*walk_num, & quad_coord(1,k), 3_8, & quad_ee_distance(1,k), elec_num) if (info /= QMCKL_SUCCESS) then exit endif quad_ee_distance(indices(k)+1,k) = 0.0d0 end do end function qmckl_compute_quad_ee_distance #+end_src #+begin_src c :tangle (eval h_private_func) :comments org :exports none qmckl_exit_code qmckl_compute_quad_ee_distance ( const qmckl_context context, const int64_t num, const int64_t* indices, const int64_t elec_num, const int64_t walk_num, const double* coord, const double* quad_coord, double* const quad_ee_distance ); #+end_src *** Test :noexport: #+begin_src c :tangle (eval c_test) printf("Quad e-e distance\n"); /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); double ee_distance[walk_num][elec_num][elec_num]; double quad_ee_distance[num][elec_num]; for (int elec = 0; elec < num; elec++){ rc = qmckl_set_electron_coord(context, 'N', walk_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], walk_num*elec_num*elec_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_quad_points(context, 'N', num, indices, &new_coords[0][0], 3*num); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_quad_electron_ee_distance(context,&quad_ee_distance[0][0],num*elec_num); assert (rc == QMCKL_SUCCESS); coords[0][indices[elec]][0] = new_coords[elec][0]; coords[0][indices[elec]][1] = new_coords[elec][1]; coords[0][indices[elec]][2] = new_coords[elec][2]; rc = qmckl_set_electron_coord(context, 'N', walk_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], walk_num*elec_num*elec_num); assert (rc == QMCKL_SUCCESS); for (int i = 0; i < elec_num; i++) { if (i == indices[elec]) continue; //printf("nw=%d elec=%ld i=%d\n",elec,indices[elec],i); //printf("ee_distance=%f quad_ee_distance=%f\n",ee_distance[0][indices[elec]][i],quad_ee_distance[elec][i]); assert(fabs((ee_distance[0][indices[elec]][i]-quad_ee_distance[elec][i])) < 1.e-12); assert(fabs((ee_distance[0][i][indices[elec]]-quad_ee_distance[elec][i])) < 1.e-12); } } printf("OK\n"); #+end_src ** Electron-nucleus distances *** Get Electron-nucleus distance between the quad electron and all nuclei for all walkers. Dimension is ~[num][nucl_num]~. #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_quad_electron_en_distance(qmckl_context context, double* distance, const int64_t size_max); #+end_src #+begin_src f90 :tangle (eval fh_func) :comments org :exports none interface integer(qmckl_exit_code) function qmckl_get_quad_electron_en_distance(context, distance, size_max) & bind(C) use, intrinsic :: iso_c_binding import implicit none integer (c_int64_t) , intent(in) , value :: context real (c_double ) , intent(out) :: distance(*) integer (c_int64_t) , intent(in) , value :: size_max end function end interface #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_quad_electron_en_distance(qmckl_context context, double* distance, const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_exit_code rc; rc = qmckl_provide_quad_en_distance(context); if (rc != QMCKL_SUCCESS) return rc; qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); if (distance == NULL) { return qmckl_failwith( context, QMCKL_INVALID_ARG_2, "qmckl_get_quad_electron_en_distance", "distance is a NULL pointer"); } int64_t sze = ctx->nucleus.num * ctx->quad_point.num; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_quad_electron_en_distance", "Array too small. Expected ctx->nucleus.num * ctx->quad_point.num"); } memcpy(distance, ctx->quad_point.quad_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_quad_en_distance(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_provide_quad_en_distance(qmckl_context context) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_exit_code rc; qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); if (!(ctx->nucleus.provided)) { return qmckl_failwith( context, QMCKL_NOT_PROVIDED, "qmckl_provide_quad_en_distance", NULL); } /* Compute if necessary */ if (ctx->quad_point.date > ctx->quad_point.quad_en_distance_date) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = ctx->nucleus.num * ctx->quad_point.num * sizeof(double); if (mem_info.size > ctx->quad_point.quad_en_distance_maxsize) { if (ctx->quad_point.quad_en_distance != NULL) { rc = qmckl_free(context, ctx->quad_point.quad_en_distance); assert (rc == QMCKL_SUCCESS); ctx->quad_point.quad_en_distance = NULL; } } /* Allocate array */ if (ctx->quad_point.quad_en_distance == NULL) { double* quad_en_distance = (double*) qmckl_malloc(context, mem_info); ctx->quad_point.quad_en_distance_maxsize = mem_info.size; if (quad_en_distance == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_quad_en_distance", NULL); } ctx->quad_point.quad_en_distance = quad_en_distance; } qmckl_exit_code rc = qmckl_compute_quad_en_distance(context, ctx->quad_point.num, ctx->quad_point.indices, ctx->nucleus.num, ctx->electron.walker.num, ctx->quad_point.coord.data, ctx->nucleus.coord.data, ctx->quad_point.quad_en_distance); if (rc != QMCKL_SUCCESS) { return rc; } ctx->quad_point.quad_en_distance_date = ctx->quad_point.date; } return QMCKL_SUCCESS; } #+end_src *** Compute :PROPERTIES: :Name: qmckl_compute_quad_en_distance :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_quad_en_distance_args |----------------------+-----------------------+--------+------------------------------------------------| | Variable | Type | In/Out | Description | |----------------------+-----------------------+--------+------------------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~num~ | ~int64_t~ | in | Number of quad electrons | | ~indices~ | ~int64_t[num]~ | in | Indices of quad electron | | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | | ~walk_num~ | ~int64_t~ | in | Number of walkers | | ~elec_coord~ | ~double[num][3]~ | in | Electron coordinates | | ~nucl_coord~ | ~double[3][nucl_num]~ | in | Nuclear coordinates | | ~quad_en_distance~ | ~double[num][nucl_num]~ | out | Electron-nucleus distances for quad-electron | |----------------------+-----------------------+--------+------------------------------------------------| #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_quad_en_distance(context, num, indices, nucl_num, walk_num, & elec_coord, nucl_coord, quad_en_distance) result(info) bind(C) use, intrinsic :: iso_c_binding use qmckl implicit none integer(qmckl_context), intent(in) :: context integer (c_int64_t) , intent(in) , value :: num integer (c_int64_t) , intent(in) :: indices(num) integer (c_int64_t) , intent(in) , value :: nucl_num, walk_num real (c_double ) , intent(in) :: elec_coord(3,num) real (c_double ) , intent(in) :: nucl_coord(nucl_num,3) real (c_double ) , intent(out) :: quad_en_distance(nucl_num, 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 do k=1,num info = qmckl_distance(context, 'T', 'N', nucl_num, 1_8, & nucl_coord(:,:), nucl_num, & elec_coord(:,k), 3_8, & quad_en_distance(:,k), nucl_num) end do end function qmckl_compute_quad_en_distance #+end_src #+begin_src c :tangle (eval h_private_func) :comments org :exports none qmckl_exit_code qmckl_compute_quad_en_distance ( const qmckl_context context, const int64_t num, const int64_t* indices, const int64_t nucl_num, const int64_t walk_num, const double* elec_coord, const double* nucl_coord, double* const quad_en_distance ); #+end_src *** Test :noexport: #+begin_src c :tangle (eval c_test) printf("Quad e-n distance\n"); /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); double en_distance[walk_num][elec_num][nucl_num]; double quad_en_distance[num][nucl_num]; for (int elec = 0; elec < num; elec++){ rc = qmckl_set_electron_coord(context, 'N', walk_num, elec_coord, walk_num*elec_num*3); assert(rc == QMCKL_SUCCESS); rc = qmckl_get_electron_en_distance(context, &en_distance[0][0][0], walk_num*elec_num*nucl_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_quad_points(context, 'N', num, indices, &new_coords[0][0], 3*num); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_quad_electron_en_distance(context,&quad_en_distance[0][0],num*nucl_num); assert (rc == QMCKL_SUCCESS); coords[0][indices[elec]][0] = new_coords[elec][0]; coords[0][indices[elec]][1] = new_coords[elec][1]; coords[0][indices[elec]][2] = new_coords[elec][2]; rc = qmckl_set_electron_coord(context, 'N', walk_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][0], walk_num*elec_num*nucl_num); assert (rc == QMCKL_SUCCESS); for (int a = 0; a < nucl_num; a++){ //printf("nw=%d elec=%ld a=%d\n",elec,indices[elec],a); //printf("en_distance=%f quad_en_distance=%f\n",en_distance[0][indices[elec]][a],quad_en_distance[elec][a]); assert(fabs((en_distance[0][indices[elec]][a]-quad_en_distance[elec][a])) < 1.e-12); } } printf("OK\n"); #+end_src * Electron-electron-nucleus Jastrow ** Electron-electron rescaled distances *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_een_rescaled_quad_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_quad_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_quad_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->quad_point.num * (ctx->jastrow_champ.cord_num + 1); if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_een_rescaled_quad_e", "Array too small. Expected ctx->electron.num * ctx->electron.num * ctx->quad_point.num * (ctx->jastrow_champ.cord_num + 1)"); } memcpy(distance_rescaled, ctx->quad_point.een_rescaled_quad_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_quad_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_quad_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_quad_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->quad_point.date > ctx->quad_point.een_rescaled_quad_e_date) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = ctx->electron.num * ctx->quad_point.num * (ctx->jastrow_champ.cord_num + 1) * sizeof(double); if (mem_info.size > ctx->quad_point.een_rescaled_quad_e_maxsize) { if (ctx->quad_point.een_rescaled_quad_e!= NULL) { rc = qmckl_free(context, ctx->quad_point.een_rescaled_quad_e); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_een_rescaled_quad_e", "Unable to free ctx->quad_point.een_rescaled_quad_e"); } ctx->quad_point.een_rescaled_quad_e = NULL; } } /* Allocate array */ if (ctx->quad_point.een_rescaled_quad_e == NULL) { double* een_rescaled_quad_e = (double*) qmckl_malloc(context, mem_info); ctx->quad_point.een_rescaled_quad_e_maxsize = mem_info.size; if (een_rescaled_quad_e == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_een_rescaled_quad_e", NULL); } ctx->quad_point.een_rescaled_quad_e = een_rescaled_quad_e; } rc = qmckl_compute_een_rescaled_quad_e(context, ctx->quad_point.num, ctx->quad_point.indices, ctx->electron.walker.num, ctx->electron.num, ctx->jastrow_champ.cord_num, ctx->jastrow_champ.rescale_factor_ee, ctx->quad_point.quad_ee_distance, ctx->jastrow_champ.een_rescaled_e, ctx->quad_point.een_rescaled_quad_e); if (rc != QMCKL_SUCCESS) { return rc; } ctx->quad_point.een_rescaled_quad_e_date = ctx->quad_point.date; } return QMCKL_SUCCESS; } #+end_src *** Compute :PROPERTIES: :Name: qmckl_compute_een_rescaled_quad_e :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_factor_een_rescaled_quad_e_args |-------------------------+----------------------------------------------------+--------+-------------------------------------------------------------| | Variable | Type | In/Out | Description | |-------------------------+----------------------------------------------------+--------+-------------------------------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~num~ | ~int64_t~ | in | Number of quad electron | | ~indices~ | ~int64_t[num]~ | in | Indices of quad electron | | ~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 | | ~quad_ee_distance~ | ~double[num][elec_num]~ | in | quad 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_quad_e~ | ~double[num][0:cord_num][elec_num]~ | out | quad electron-electron rescaled distances for each walker | |-------------------------+----------------------------------------------------+--------+-------------------------------------------------------------| #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_een_rescaled_quad_e_doc( & context, num, indices, walk_num, elec_num, cord_num, rescale_factor_ee, & quad_ee_distance, een_rescaled_e, een_rescaled_quad_e) & result(info) bind(C) use, intrinsic :: iso_c_binding use qmckl implicit none integer(qmckl_context), intent(in) :: context integer(c_int64_t) , intent(in), value :: num integer(c_int64_t) , intent(in) :: indices(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) :: quad_ee_distance(elec_num,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_quad_e(elec_num,0:cord_num,num) double precision,allocatable :: een_rescaled_quad_e_ij(:,:) double precision :: x integer*8 :: i, j, k, l, nw info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif if (walk_num <= 0) then info = QMCKL_INVALID_ARG_2 return endif if (elec_num <= 0) then info = QMCKL_INVALID_ARG_3 return endif if (cord_num < 0) then info = QMCKL_INVALID_ARG_4 return endif allocate(een_rescaled_quad_e_ij(elec_num, cord_num + 1)) ! Prepare table of exponentiated distances raised to appropriate power do nw = 1, num een_rescaled_quad_e_ij(:, 1) = 1.0d0 do j = 1, elec_num een_rescaled_quad_e_ij(j, 2) = dexp(-rescale_factor_ee * quad_ee_distance(j, nw)) end do do l = 2, cord_num do k = 1, elec_num een_rescaled_quad_e_ij(k, l + 1) = een_rescaled_quad_e_ij(k, l) * een_rescaled_quad_e_ij(k, 2) end do end do ! prepare the actual een table een_rescaled_quad_e(:,0,nw) = 1.0d0 do l = 1, cord_num do j = 1, elec_num x = een_rescaled_quad_e_ij(j, l + 1) een_rescaled_quad_e(j, l, nw) = x end do end do een_rescaled_quad_e(indices(nw)+1, :, nw) = 0.0d0 end do end function qmckl_compute_een_rescaled_quad_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_quad_e ( const qmckl_context context, const int64_t num, const int64_t* indices, const int64_t walk_num, const int64_t elec_num, const int64_t cord_num, const double rescale_factor_ee, const double* quad_ee_distance, const double* een_rescaled_e, double* const een_rescaled_quad_e ); #+end_src #+begin_src c :tangle (eval h_private_func) :comments org qmckl_exit_code qmckl_compute_een_rescaled_quad_e_doc ( const qmckl_context context, const int64_t num, const int64_t* indices, const int64_t walk_num, const int64_t elec_num, const int64_t cord_num, const double rescale_factor_ee, const double* quad_ee_distance, const double* een_rescaled_e, double* const een_rescaled_quad_e ); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes qmckl_exit_code qmckl_compute_een_rescaled_quad_e (const qmckl_context context, const int64_t num, const int64_t* indices, const int64_t walk_num, const int64_t elec_num, const int64_t cord_num, const double rescale_factor_ee, const double* quad_ee_distance, const double* een_rescaled_e, double* const een_rescaled_quad_e ) { #ifdef HAVE_HPC return qmckl_compute_een_rescaled_quad_e_doc #else return qmckl_compute_een_rescaled_quad_e_doc #endif (context, num, indices, walk_num, elec_num, cord_num, rescale_factor_ee, quad_ee_distance, een_rescaled_e, een_rescaled_quad_e); } #+end_src *** Test :noexport: #+begin_src c :tangle (eval c_test) printf("quad een rescaled e-e distance\n"); /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); double rescaled_een_ee_distance[walk_num][cord_num+1][elec_num][elec_num]; double quad_rescaled_een_ee_distance[num][cord_num+1][elec_num]; for (int elec = 0; elec < num; elec++){ rc = qmckl_set_electron_coord(context, 'N', walk_num, elec_coord, 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); rc = qmckl_get_electron_coord(context, 'N', &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); rc = qmckl_set_quad_points(context, 'N', num, indices, &new_coords[0][0], 3*num); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_een_rescaled_quad_e(context, &quad_rescaled_een_ee_distance[0][0][0], num*(cord_num+1)*elec_num); assert (rc == QMCKL_SUCCESS); coords[0][indices[elec]][0] = new_coords[elec][0]; coords[0][indices[elec]][1] = new_coords[elec][1]; coords[0][indices[elec]][2] = new_coords[elec][2]; rc = qmckl_set_electron_coord(context, 'N', walk_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 l = 0; l <= cord_num; l++){ for (int i = 0; i < elec_num; i++) { if (i == indices[elec]) continue; //printf("nw=%d elec=%ld i=%d\n",elec,indices[elec],i); //printf("rescaled_een_ee_distance=%f quad_rescaled_een_ee_distance=%f\n",rescaled_een_ee_distance[0][l][indices[elec]][i],quad_rescaled_een_ee_distance[elec][l][i]); assert(fabs((rescaled_een_ee_distance[0][l][indices[elec]][i]-quad_rescaled_een_ee_distance[elec][l][i])) < 1.e-12); assert(fabs((rescaled_een_ee_distance[0][l][i][indices[elec]]-quad_rescaled_een_ee_distance[elec][l][i])) < 1.e-12); } } } printf("OK\n"); #+end_src ** Electron-nucleus rescaled distances *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_een_rescaled_quad_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_quad_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_quad_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->quad_point.num * (ctx->jastrow_champ.cord_num + 1); if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "mckl_get_een_rescaled_quad_n", "Array too small. Expected ctx->nucleus.num * ctx->quad_point.num * (ctx->jastrow_champ.cord_num + 1)"); } memcpy(distance_rescaled, ctx->quad_point.een_rescaled_quad_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_quad_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_quad_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_quad_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->quad_point.date > ctx->quad_point.een_rescaled_quad_n_date) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = ctx->nucleus.num * ctx->quad_point.num * (ctx->jastrow_champ.cord_num + 1) * sizeof(double); if (mem_info.size > ctx->quad_point.een_rescaled_quad_n_maxsize) { if (ctx->quad_point.een_rescaled_quad_n != NULL) { rc = qmckl_free(context, ctx->quad_point.een_rescaled_quad_n); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_een_rescaled_quad_n", "Unable to free ctx->quad_point.een_rescaled_quad_n"); } ctx->quad_point.een_rescaled_quad_n = NULL; } } /* Allocate array */ if (ctx->quad_point.een_rescaled_quad_n == NULL) { double* een_rescaled_quad_n = (double*) qmckl_malloc(context, mem_info); ctx->quad_point.een_rescaled_quad_n_maxsize = mem_info.size; if (een_rescaled_quad_n == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_een_rescaled_quad_n", NULL); } ctx->quad_point.een_rescaled_quad_n = een_rescaled_quad_n; } rc = qmckl_compute_een_rescaled_quad_n(context, ctx->quad_point.num, ctx->quad_point.indices, 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->quad_point.quad_en_distance, ctx->jastrow_champ.een_rescaled_n, ctx->quad_point.een_rescaled_quad_n); if (rc != QMCKL_SUCCESS) { return rc; } ctx->quad_point.een_rescaled_quad_n_date = ctx->quad_point.date; } return QMCKL_SUCCESS; } #+end_src *** Compute :PROPERTIES: :Name: qmckl_compute_een_rescaled_quad_n :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_factor_een_rescaled_quad_n_args | Variable | Type | In/Out | Description | |-------------------------+----------------------------------------------------+--------+--------------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~num~ | ~int64_t~ | in | Number of quad electron | | ~indices~ | ~int64_t[num]~ | in | Indices of quad 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 | | ~quad_en_distance~ | ~double[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_quad_n~ | ~double[num][0:cord_num][nucl_num]~ | out | quad electron-nucleus rescaled distances | |-------------------------+----------------------------------------------------+--------+--------------------------------------------| #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_een_rescaled_quad_n( & context, num, indices, walk_num, elec_num, nucl_num, & type_nucl_num, type_nucl_vector, cord_num, rescale_factor_en, & quad_en_distance, een_rescaled_n, een_rescaled_quad_n) & result(info) bind(C) use, intrinsic :: iso_c_binding use qmckl implicit none integer(qmckl_context), intent(in) :: context integer(c_int64_t) , intent(in), value :: num integer(c_int64_t) , intent(in) :: indices(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(type_nucl_num) real(c_double) , intent(in) :: quad_en_distance(nucl_num,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_quad_n(nucl_num,0:cord_num,num) double precision :: x integer*8 :: i, a, k, l, nw info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif if (walk_num <= 0) then info = QMCKL_INVALID_ARG_2 return endif if (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, num ! prepare the actual een table een_rescaled_quad_n(:, 0, nw) = 1.0d0 do a = 1, nucl_num een_rescaled_quad_n(a, 1, nw) = dexp(-rescale_factor_en(type_nucl_vector(a)+1) * quad_en_distance(a, nw)) end do do l = 2, cord_num do a = 1, nucl_num een_rescaled_quad_n(a, l, nw) = een_rescaled_quad_n(a, l - 1, nw) * een_rescaled_quad_n(a, 1, nw) end do end do end do end function qmckl_compute_een_rescaled_quad_n #+end_src #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_een_rescaled_quad_n ( const qmckl_context context, const int64_t num, const int64_t* indices, 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* quad_en_distance, const double* een_rescaled_n, double* const een_rescaled_quad_n ); #+end_src *** Test :noexport: #+begin_src c :tangle (eval c_test) printf("Quad een rescaled e-n distance\n"); /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); double rescaled_een_en_distance[walk_num][cord_num+1][nucl_num][elec_num]; double quad_rescaled_een_en_distance[num][cord_num+1][nucl_num]; for (int elec = 0; elec < num; elec++){ rc = qmckl_set_electron_coord(context, 'N', walk_num, elec_coord, 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); rc = qmckl_get_electron_coord(context, 'N', &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); rc = qmckl_set_quad_points(context, 'N', num, indices, &new_coords[0][0], 3*num); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_een_rescaled_quad_n(context, &quad_rescaled_een_en_distance[0][0][0], num*(cord_num+1)*nucl_num); assert (rc == QMCKL_SUCCESS); coords[0][indices[elec]][0] = new_coords[elec][0]; coords[0][indices[elec]][1] = new_coords[elec][1]; coords[0][indices[elec]][2] = new_coords[elec][2]; rc = qmckl_set_electron_coord(context, 'N', walk_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 l = 0; l <= cord_num; l++){ for (int a = 0; a < nucl_num; a++) { //printf("nw=%d elec=%ld a=%d\n",elec,indices[elec],a); //printf("rescaled_een_en_distance=%f quad_rescaled_een_en_distance=%f\n",rescaled_een_en_distance[0][l][a][indices[elec]],quad_rescaled_een_en_distance[elec][l][a]); assert(fabs((rescaled_een_en_distance[0][l][a][indices[elec]]-quad_rescaled_een_en_distance[elec][l][a])) < 1.e-12); } } } printf("OK\n"); #+end_src ** $\delta P$ matrix *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_jastrow_champ_quad_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_quad_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_quad_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->quad_point.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->quad_point.delta_p, 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_jastrow_champ_quad_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_quad_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_quad_e(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_een_rescaled_quad_n(context); if(rc != QMCKL_SUCCESS) return rc; /* Compute if necessary */ if (ctx->quad_point.date > ctx->quad_point.delta_p_date) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = ctx->quad_point.num * ctx->jastrow_champ.cord_num * (ctx->jastrow_champ.cord_num + 1) * ctx->nucleus.num * ctx->electron.num * sizeof(double); if (mem_info.size > ctx->quad_point.delta_p_maxsize) { if (ctx->quad_point.delta_p != NULL) { rc = qmckl_free(context, ctx->quad_point.delta_p); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_jastrow_champ_delta_p", "Unable to free ctx->quad_point.delta_p"); } ctx->quad_point.delta_p = NULL; } } /* Allocate array */ if (ctx->quad_point.delta_p == NULL) { double* delta_p = (double*) qmckl_malloc(context, mem_info); ctx->quad_point.delta_p_maxsize = mem_info.size; if (delta_p == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_jastrow_champ_delta_p", NULL); } ctx->quad_point.delta_p = delta_p; } rc = qmckl_compute_jastrow_champ_quad_delta_p_doc(context, ctx->quad_point.num, ctx->quad_point.indices, 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->quad_point.een_rescaled_quad_n, ctx->quad_point.een_rescaled_quad_e, ctx->quad_point.delta_p); if (rc != QMCKL_SUCCESS) { return rc; } ctx->quad_point.delta_p_date = ctx->quad_point.date; } return QMCKL_SUCCESS; } #+end_src *** Compute :PROPERTIES: :Name: qmckl_compute_jastrow_champ_quad_delta_p_doc :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_factor_quad_delta_p_args | Variable | Type | In/Out | Description | |-------------------------+------------------------------------------------------------------+--------+---------------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~num~ | ~int64_t~ | in | quad point index | | ~indices~ | ~int64_t[num]~ | in | Indices of quad points | | ~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_quad_n~ | ~double[num][0:cord_num][nucl_num]~ | in | Electron-nucleus quad rescaled distances | | ~een_rescaled_quad_e~ | ~double[num][0:cord_num][elec_num]~ | in | Electron-electron quad rescaled distances | | ~delta_p~ | ~double[num][0:cord_num-1][0:cord_num][nucl_num][elec_num]~ | out | quad point matrix P | |-------------------------+------------------------------------------------------------------+--------+---------------------------------------------| #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_jastrow_champ_quad_delta_p_doc( & context, num, indices, walk_num, elec_num, nucl_num, cord_num, & een_rescaled_n, een_rescaled_e, een_rescaled_quad_n, een_rescaled_quad_e, delta_p) & result(info) bind(C) use, intrinsic :: iso_c_binding use qmckl implicit none integer(qmckl_context), intent(in) :: context integer(c_int64_t), intent(in), value :: num, walk_num, elec_num, cord_num, nucl_num integer(c_int64_t) , intent(in) :: indices(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_quad_n(nucl_num, 0:cord_num, num) real(c_double) , intent(in) :: een_rescaled_quad_e(elec_num, 0:cord_num, num) real(c_double) , intent(out) :: delta_p(elec_num, nucl_num,0:cord_num, 0:cord_num-1, num) double precision :: een_rescaled_delta_e(elec_num) integer*8 :: i, a, c, j, l, k, p, m, n, nw, idx double precision :: dn, dn2 integer*8 :: LDA, LDB, LDC 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 do nw=1, num idx = indices(nw)+1 do i=0, cord_num-1 een_rescaled_delta_e(:) = een_rescaled_quad_e(:,i,nw) - een_rescaled_e(:,idx,i,1) do c=0,cord_num do a=1,nucl_num dn = een_rescaled_quad_n(a,c,nw) - een_rescaled_n(idx,a,c,1) dn2 = een_rescaled_quad_n(a,c,nw) do j=1,elec_num delta_p(j,a,c,i,nw) = een_rescaled_e(j,idx,i,1)*dn + een_rescaled_delta_e(j) * dn2 enddo end do end do info = qmckl_dgemm(context, 'T', 'N', 1_8, nucl_num * (cord_num+1_8), elec_num, 1.0d0, & een_rescaled_delta_e,elec_num, & een_rescaled_n(1,1,0,1),elec_num, & 1.0d0, & delta_p(idx,1,0,i,nw),elec_num) enddo end do end function qmckl_compute_jastrow_champ_quad_delta_p_doc #+end_src #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_jastrow_champ_quad_delta_p_doc (const qmckl_context context, const int64_t num, const int64_t* indices, 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_quad_n, const double* een_rescaled_quad_e, double* const delta_p ); qmckl_exit_code qmckl_compute_jastrow_champ_quad_delta_p (const qmckl_context context, const int64_t num, const int64_t* indices, 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_quad_n, const double* een_rescaled_quad_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_quad_delta_p (const qmckl_context context, const int64_t num, const int64_t* indices, 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_quad_n, const double* een_rescaled_quad_e, double* const delta_p ) { #ifdef HAVE_HPC return qmckl_compute_jastrow_champ_quad_delta_p_doc #else return qmckl_compute_jastrow_champ_quad_delta_p_doc #endif (context, num, indices, walk_num, elec_num, nucl_num, cord_num, een_rescaled_n, een_rescaled_e, een_rescaled_quad_n, een_rescaled_quad_e, delta_p ); } #+end_src *** Test :noexport: #+begin_src c :tangle (eval c_test) printf("Delta p\n"); /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); double p_old[walk_num][cord_num][cord_num+1][nucl_num][elec_num]; double delta_p[num][cord_num][cord_num+1][nucl_num][elec_num]; double p_new[walk_num][cord_num][cord_num+1][nucl_num][elec_num]; for (int elec = 0; elec < num; elec++){ rc = qmckl_set_electron_coord(context, 'N', walk_num, elec_coord, walk_num*elec_num*3); assert(rc == QMCKL_SUCCESS); rc = qmckl_get_jastrow_champ_tmp_c(context, &p_old[0][0][0][0][0], walk_num*cord_num*(cord_num+1)*nucl_num*elec_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_quad_points(context, 'N', num, indices, &new_coords[0][0], 3*num); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_jastrow_champ_quad_delta_p(context, &delta_p[0][0][0][0][0], num*cord_num*(cord_num+1)*nucl_num*elec_num); assert (rc == QMCKL_SUCCESS); coords[0][indices[elec]][0] = new_coords[elec][0]; coords[0][indices[elec]][1] = new_coords[elec][1]; coords[0][indices[elec]][2] = new_coords[elec][2]; rc = qmckl_set_electron_coord(context, 'N', walk_num, &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_jastrow_champ_tmp_c(context, &p_new[0][0][0][0][0], walk_num*cord_num*(cord_num+1)*nucl_num*elec_num); assert (rc == QMCKL_SUCCESS); 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++){ //printf("nw=%d elec=%ld l=%d m=%d a=%d i=%d\n",elec,indices[elec],l,m,a,i); //printf("p_new=%f p_old=%f delta_p=%f\n",p_new[0][l][m][a][i],p_old[0][l][m][a][i],delta_p[elec][l][m][a][i]); assert(fabs(((p_new[0][l][m][a][i]-p_old[0][l][m][a][i])-delta_p[elec][l][m][a][i])) < 1.e-12); } } } } } printf("OK\n"); #+end_src ** Electron-electron-nucleus Jastrow value *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_jastrow_champ_quad_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_quad_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_quad_een(context); if (rc != QMCKL_SUCCESS) return rc; qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int64_t sze = ctx->quad_point.num; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_jastrow_champ_quad_eenn", "Array too small. Expected ctx->quad_point.num"); } memcpy(delta_een, ctx->quad_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_quad_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_quad_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_quad_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; if (ctx->jastrow_champ.cord_num > 0) { rc = qmckl_provide_jastrow_champ_quad_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->quad_point.date > ctx->quad_point.delta_een_date) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = ctx->quad_point.num * sizeof(double); if (mem_info.size > ctx->quad_point.delta_een_maxsize) { if (ctx->quad_point.delta_een != NULL) { rc = qmckl_free(context, ctx->quad_point.delta_een); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_jastrow_champ_quad_een", "Unable to free ctx->quad_point.delta_een"); } ctx->quad_point.delta_een = NULL; } } /* Allocate array */ if (ctx->quad_point.delta_een == NULL) { double* delta_een = (double*) qmckl_malloc(context, mem_info); ctx->quad_point.delta_een_maxsize = mem_info.size; if (delta_een == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_jastrow_champ_quad_een", NULL); } ctx->quad_point.delta_een = delta_een; } rc = qmckl_compute_jastrow_champ_factor_quad_een_doc(context, ctx->quad_point.num, ctx->quad_point.indices, 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->quad_point.delta_p, ctx->jastrow_champ.een_rescaled_n, ctx->jastrow_champ.een_rescaled_e, ctx->quad_point.een_rescaled_quad_n, ctx->quad_point.een_rescaled_quad_e, ctx->quad_point.delta_een); if (rc != QMCKL_SUCCESS) { return rc; } ctx->quad_point.delta_een_date = ctx->quad_point.date; } return QMCKL_SUCCESS; } #+end_src *** Compute :PROPERTIES: :Name: qmckl_compute_jastrow_champ_factor_quad_een_doc :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_factor_quad_een_args | Variable | Type | In/Out | Description | |-------------------------+------------------------------------------------------------------+--------+---------------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~num~ | ~int64_t~ | in | quad point number | | ~indices~ | ~int64_t[num]~ | in | Indices of quad points | | ~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 | P matrix | | ~delta_p~ | ~double[num][0:cord_num-1][0:cord_num][nucl_num][elec_num]~ | in | quad electron P matrix | | ~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_quad_n~ | ~double[num][0:cord_num][nucl_num]~ | in | Electron-nucleus quad rescaled distances | | ~een_rescaled_quad_e~ | ~double[wnum][0:cord_num][elec_num]~ | in | Electron-electron quad rescaled distances | | ~delta_een~ | ~double[num]~ | out | Electron-nucleus jastrow | |-------------------------+------------------------------------------------------------------+--------+---------------------------------------------| #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_jastrow_champ_factor_quad_een_doc( & context, num, indices, 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_quad_n, & een_rescaled_quad_e, delta_een) & result(info) bind(C) use, intrinsic :: iso_c_binding use qmckl implicit none integer(qmckl_context), intent(in) :: context integer(c_int64_t) , intent(in), value :: num, walk_num, elec_num, cord_num, nucl_num, dim_c_vector integer(c_int64_t) , intent(in) :: indices(num) integer(c_int64_t) , intent(in) :: lkpm_combined_index(dim_c_vector,4) real(c_double) , intent(in) :: c_vector_full(nucl_num, dim_c_vector) 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, 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_quad_n(nucl_num, 0:cord_num, num) real(c_double) , intent(in) :: een_rescaled_quad_e(elec_num, 0:cord_num, num) real(c_double) , intent(out) :: delta_een(num) double precision :: een_rescaled_delta_n(nucl_num, 0:cord_num) integer*8 :: i, a, j, l, k, p, m, n, nw, idx double precision :: accu, accu2, cn integer*8 :: LDA, LDB, LDC 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 if (cord_num == 0) return do nw =1, num idx= indices(nw)+1 een_rescaled_delta_n(:,:) = een_rescaled_quad_n(:,:,nw) - een_rescaled_n(idx,:,:,1) 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,1) * delta_p(j,a,m+l,k,nw) end do accu = accu + een_rescaled_delta_n(a,m) * (tmp_c(idx,a,m+l,k,1) + delta_p(idx,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_quad_een_doc #+end_src #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_jastrow_champ_factor_quad_een_doc (const qmckl_context context, const int64_t num, const int64_t* indices, 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_quad_n, const double* een_rescaled_quad_e, double* const delta_een ); qmckl_exit_code qmckl_compute_jastrow_champ_factor_quad_een (const qmckl_context context, const int64_t num, const int64_t* indices, 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_quad_n, const double* een_rescaled_quad_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_quad_een (const qmckl_context context, const int64_t num, const int64_t* indices, 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_quad_n, const double* een_rescaled_quad_e, double* const delta_een ) { #ifdef HAVE_HPC return qmckl_compute_jastrow_champ_factor_quad_een_doc #else return qmckl_compute_jastrow_champ_factor_quad_een_doc #endif (context, num, indices, 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_quad_n, een_rescaled_quad_e, delta_een ); } #+end_src *** Test :noexport: #+begin_src c :tangle (eval c_test) printf("Delta een\n"); /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); double jastrow_een_old[walk_num]; double delta_een[num]; double jastrow_een_new[walk_num]; for (int elec = 0; elec < num; elec++){ rc = qmckl_set_electron_coord(context, 'N', walk_num, elec_coord, walk_num*elec_num*3); assert(rc == QMCKL_SUCCESS); 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_quad_points(context, 'N', num, indices, &new_coords[0][0], 3*num); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_jastrow_champ_quad_een(context, &delta_een[0], num); assert (rc == QMCKL_SUCCESS); coords[0][indices[elec]][0] = new_coords[elec][0]; coords[0][indices[elec]][1] = new_coords[elec][1]; coords[0][indices[elec]][2] = new_coords[elec][2]; rc = qmckl_set_electron_coord(context, 'N', walk_num, &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_jastrow_champ_factor_een(context, &jastrow_een_new[0], walk_num); assert (rc == QMCKL_SUCCESS); //printf("jastrow_een_old %f\n", jastrow_een_old[0]); //printf("jastrow_een_new %f\n", jastrow_een_new[0]); //printf("delta_een %f\n", delta_een[elec]); assert(fabs((jastrow_een_new[0]-jastrow_een_old[0])-delta_een[elec]) < 1.e-12); } printf("OK\n"); #+end_src ** Electron-nucleus rescaled distance derivative *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_een_rescaled_quad_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_quad_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_quad_n_gl(context); if (rc != QMCKL_SUCCESS) return rc; qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int64_t sze = 3 * ctx->nucleus.num * ctx->quad_point.num * (ctx->jastrow_champ.cord_num + 1); if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_jastrow_champ_quad_een_gl", "Array too small. Expected 3 * ctx->nucleus.num * ctx->quad_point.num * (ctx->jastrow_champ.cord_num + 1)"); } memcpy(distance_rescaled, ctx->quad_point.een_rescaled_quad_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_quad_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_quad_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 en distance is provided */ qmckl_exit_code rc = qmckl_provide_quad_en_distance(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if ee distance is provided */ rc = qmckl_provide_een_rescaled_quad_n(context); if(rc != QMCKL_SUCCESS) return rc; /* Compute if necessary */ if (ctx->quad_point.date > ctx->quad_point.een_rescaled_quad_n_gl_date) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = 3 * ctx->nucleus.num * ctx->quad_point.num * (ctx->jastrow_champ.cord_num + 1) * sizeof(double); if (mem_info.size > ctx->quad_point.een_rescaled_quad_n_gl_maxsize) { if (ctx->quad_point.een_rescaled_quad_n_gl != NULL) { rc = qmckl_free(context, ctx->quad_point.een_rescaled_quad_n_gl); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_een_rescaled_quad_n_gl", "Unable to free ctx->quad_pont.een_rescaled_quad_n_gl"); } ctx->quad_point.een_rescaled_quad_n_gl = NULL; } } /* Allocate array */ if (ctx->quad_point.een_rescaled_quad_n_gl == NULL) { double* een_rescaled_quad_n_gl = (double*) qmckl_malloc(context, mem_info); ctx->quad_point.een_rescaled_quad_n_gl_maxsize = mem_info.size; if (een_rescaled_quad_n_gl == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_een_rescaled_quad_n_gl", NULL); } ctx->quad_point.een_rescaled_quad_n_gl = een_rescaled_quad_n_gl; } rc = qmckl_compute_een_rescaled_quad_n_gl(context, ctx->quad_point.num, ctx->quad_point.indices, 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->quad_point.coord.data, ctx->nucleus.coord.data, ctx->quad_point.quad_en_distance, ctx->quad_point.een_rescaled_quad_n, ctx->quad_point.een_rescaled_quad_n_gl); if (rc != QMCKL_SUCCESS) { return rc; } ctx->quad_point.een_rescaled_quad_n_gl_date = ctx->quad_point.date; } return QMCKL_SUCCESS; } #+end_src *** Compute :PROPERTIES: :Name: qmckl_compute_een_rescaled_quad_n_gl :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_compute_een_rescaled_quad_n_gl_args |----------------------------+---------------------------------------------+--------+-------------------------------------------------------| | Variable | Type | In/Out | Description | |----------------------------+---------------------------------------------+--------+-------------------------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~num~ | ~int64_t~ | in | Number of quad points | | ~indices~ | ~int64_t[num]~ | in | Indices of quad points | | ~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[num][3]~ | in | Electron coordinates | | ~coord_n~ | ~double[3][nucl_num]~ | in | Nuclear coordinates | | ~quad_en_distance~ | ~double[num][nucl_num]~ | in | Electron-nucleus quad distances | | ~een_rescaled_quad_n~ | ~double[num][0:cord_num][nucl_num]~ | in | Electron-nucleus rescaled quad distances | | ~een_rescaled_quad_n_gl~ | ~double[num][0:cord_num][nucl_num][3]~ | out | Electron-nucleus rescaled quad distances derivative | |----------------------------+---------------------------------------------+--------+-------------------------------------------------------| #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_een_rescaled_quad_n_gl( & context, num, indices, walk_num, nucl_num, type_nucl_num, type_nucl_vector, & cord_num, rescale_factor_en, coord_ee, coord_n, quad_en_distance, & een_rescaled_quad_n, een_rescaled_quad_n_gl) & result(info) bind(C) use, intrinsic :: iso_c_binding use qmckl implicit none integer(qmckl_context), intent(in), value :: context integer(c_int64_t) , intent(in), value :: num integer(c_int64_t) , intent(in) :: indices(num) 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(type_nucl_num) real(c_double) , intent(in) :: coord_ee(3,num) real(c_double) , intent(in) :: coord_n(nucl_num,3) real(c_double) , intent(in) :: quad_en_distance(nucl_num,num) real(c_double) , intent(in) :: een_rescaled_quad_n(nucl_num,0:cord_num,num) real(c_double) , intent(out) :: een_rescaled_quad_n_gl(3,nucl_num,0:cord_num,num) double precision,allocatable :: elnuc_dist_gl(:,:) double precision :: x, ria_inv, kappa_l integer*8 :: i, a, k, l, nw, ii allocate(elnuc_dist_gl(3, 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_quad_n_gl = 0.0d0 do nw = 1, num ! prepare the actual een table do a = 1, nucl_num ria_inv = 1.0d0 / quad_en_distance(a, nw) do ii = 1, 3 elnuc_dist_gl(ii, a) = (coord_ee(ii,nw) - coord_n(a, ii)) * ria_inv end do 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_quad_n_gl(1, a, l, nw) = kappa_l * elnuc_dist_gl(1, a) een_rescaled_quad_n_gl(2, a, l, nw) = kappa_l * elnuc_dist_gl(2, a) een_rescaled_quad_n_gl(3, a, l, nw) = kappa_l * elnuc_dist_gl(3, a) een_rescaled_quad_n_gl(1, a, l, nw) = een_rescaled_quad_n_gl(1, a, l, nw) * & een_rescaled_quad_n(a, l, nw) een_rescaled_quad_n_gl(2, a, l, nw) = een_rescaled_quad_n_gl(2, a, l, nw) * & een_rescaled_quad_n(a, l, nw) een_rescaled_quad_n_gl(3, a, l, nw) = een_rescaled_quad_n_gl(3, a, l, nw) * & een_rescaled_quad_n(a, l, nw) end do end do end do end function qmckl_compute_een_rescaled_quad_n_gl #+end_src #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_een_rescaled_quad_n_gl ( const qmckl_context context, const int64_t num, const int64_t* indices, 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* quad_en_distance, const double* een_rescaled_quad_n, double* const een_rescaled_quad_n_gl ); #+end_src *** Test :noexport: #+begin_src c :tangle (eval c_test) printf("Een rescaled quad n gl\n"); /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); double een_rescaled_en_gl[walk_num][cord_num+1][nucl_num][4][elec_num]; double een_rescaled_quad_n_gl[num][cord_num+1][nucl_num][3]; for (int elec = 0; elec < num; elec++){ rc = qmckl_set_electron_coord(context, 'N', walk_num, elec_coord, 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); rc = qmckl_get_electron_coord(context, 'N', &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); rc = qmckl_set_quad_points(context, 'N', num, indices, &new_coords[0][0], 3*num); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_een_rescaled_quad_n_gl(context, &een_rescaled_quad_n_gl[0][0][0][0], num*(cord_num+1)*nucl_num*3); assert (rc == QMCKL_SUCCESS); coords[0][indices[elec]][0] = new_coords[elec][0]; coords[0][indices[elec]][1] = new_coords[elec][1]; coords[0][indices[elec]][2] = new_coords[elec][2]; rc = qmckl_set_electron_coord(context, 'N', walk_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 a = 0; a < nucl_num; a++) { for (int m = 0; m < 3; m++) { //printf("nw %d l %d a %d m %d\n", indices[elec], l, a, m); //printf(" %f %f\n", een_rescaled_en_gl[0][l][a][m][indices[elec]], een_rescaled_quad_n_gl[elec][l][a][m]); assert(fabs(een_rescaled_en_gl[0][l][a][m][indices[elec]] - een_rescaled_quad_n_gl[elec][l][a][m]) < 1.e-12); } } } } printf("OK\n"); #+end_src ** Electron-electron rescaled distances derivative *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_een_rescaled_quad_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_quad_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_quad_e_gl(context); if (rc != QMCKL_SUCCESS) return rc; qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int64_t sze = 3 * ctx->electron.num * ctx->quad_point.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 3 * ctx->electron.num * ctx->quad_point.num * (ctx->jastrow_champ.cord_num + 1)"); } memcpy(distance_rescaled, ctx->quad_point.een_rescaled_quad_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_quad_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_quad_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_quad_e(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_quad_ee_distance(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_quad_en_distance(context); if(rc != QMCKL_SUCCESS) return rc; /* Compute if necessary */ if (ctx->quad_point.date > ctx->quad_point.een_rescaled_quad_e_gl_date) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = 3 * ctx->electron.num * ctx->quad_point.num * (ctx->jastrow_champ.cord_num + 1) * sizeof(double); if (mem_info.size > ctx->quad_point.een_rescaled_quad_e_gl_maxsize) { if (ctx->quad_point.een_rescaled_quad_e_gl != NULL) { rc = qmckl_free(context, ctx->quad_point.een_rescaled_quad_e_gl); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_een_rescaled_e_gl", "Unable to free ctx->quad_point.een_rescaled_quad_e_gl"); } ctx->quad_point.een_rescaled_quad_e_gl = NULL; } } /* Allocate array */ if (ctx->quad_point.een_rescaled_quad_e_gl == NULL) { double* een_rescaled_quad_e_gl = (double*) qmckl_malloc(context, mem_info); ctx->quad_point.een_rescaled_quad_e_gl_maxsize = mem_info.size; if (een_rescaled_quad_e_gl == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_een_rescaled_quad_e_gl", NULL); } ctx->quad_point.een_rescaled_quad_e_gl = een_rescaled_quad_e_gl; } rc = qmckl_compute_een_rescaled_quad_e_gl(context, ctx->quad_point.num, ctx->quad_point.indices, ctx->electron.walker.num, ctx->electron.num, ctx->jastrow_champ.cord_num, ctx->jastrow_champ.rescale_factor_ee, ctx->quad_point.coord.data, ctx->electron.walker.point.coord.data, ctx->quad_point.quad_ee_distance, ctx->quad_point.een_rescaled_quad_e, ctx->quad_point.een_rescaled_quad_e_gl); if (rc != QMCKL_SUCCESS) { return rc; } ctx->quad_point.een_rescaled_quad_e_gl_date = ctx->quad_point.date; } return QMCKL_SUCCESS; } #+end_src *** Compute :PROPERTIES: :Name: qmckl_compute_een_rescaled_quad_e_gl :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_een_rescaled_quad_e_gl_args |----------------------------+---------------------------------------------+--------+--------------------------------------------------------| | Variable | Type | In/Out | Description | |----------------------------+---------------------------------------------+--------+--------------------------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~num~ | ~int64_t~ | in | Index of quad electron | | ~indices~ | ~int64_t[num]~ | in | Indices of quad points | | ~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[num][3]~ | in | quad electron coordinates | | ~coord_ee~ | ~double[3][walk_num][elec_num]~ | in | Electron coordinates | | ~quad_ee_distance~ | ~double[num][elec_num]~ | in | Electron-electron quad distances | | ~een_rescaled_quad_e~ | ~double[num][0:cord_num][elec_num]~ | in | Electron-electron rescaled quad distances | | ~een_rescaled_quad_e_gl~ | ~double[num][0:cord_num][elec_num][3]~ | out | Electron-electron rescaled quad distances derivative | |----------------------------+---------------------------------------------+--------+--------------------------------------------------------| #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_een_rescaled_quad_e_gl_doc( & context, num, indices, walk_num, elec_num, cord_num, rescale_factor_ee, & coord, coord_ee, quad_ee_distance, een_rescaled_quad_e, een_rescaled_quad_e_gl) & result(info) bind(C) use, intrinsic :: iso_c_binding use qmckl implicit none integer(qmckl_context), intent(in) :: context integer(c_int64_t) , intent(in), value :: num integer(c_int64_t) , intent(in) :: indices(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,num) real(c_double) , intent(in) :: coord_ee(elec_num,walk_num,3) real(c_double) , intent(in) :: quad_ee_distance(elec_num,num) real(c_double) , intent(in) :: een_rescaled_quad_e(elec_num,0:cord_num,num) real(c_double) , intent(out) :: een_rescaled_quad_e_gl(3,elec_num,0:cord_num,num) double precision,allocatable :: elec_dist_gl(:,:) double precision :: x, rij_inv, kappa_l integer*8 :: i, j, k, l, nw, ii 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_quad_e_gl = 0.0d0 ! Prepare table of exponentiated distances raised to appropriate power do nw = 1, num do i = 1, elec_num if (i == indices(nw)+1) cycle rij_inv = 1.0d0 / quad_ee_distance(i, nw) do ii = 1, 3 elec_dist_gl(ii, i) = (coord(ii, nw) - coord_ee(i, 1, ii)) * rij_inv end do end do elec_dist_gl(:, indices(nw)+1) = 0.0d0 do l = 1, cord_num kappa_l = - dble(l) * rescale_factor_ee do i = 1, elec_num een_rescaled_quad_e_gl(1, i, l, nw) = kappa_l * elec_dist_gl(1, i) een_rescaled_quad_e_gl(2, i, l, nw) = kappa_l * elec_dist_gl(2, i) een_rescaled_quad_e_gl(3, i, l, nw) = kappa_l * elec_dist_gl(3, i) een_rescaled_quad_e_gl(1,i,l,nw) = een_rescaled_quad_e_gl(1,i,l,nw) * een_rescaled_quad_e(i,l,nw) een_rescaled_quad_e_gl(2,i,l,nw) = een_rescaled_quad_e_gl(2,i,l,nw) * een_rescaled_quad_e(i,l,nw) een_rescaled_quad_e_gl(3,i,l,nw) = een_rescaled_quad_e_gl(3,i,l,nw) * een_rescaled_quad_e(i,l,nw) end do end do end do end function qmckl_compute_een_rescaled_quad_e_gl_doc #+end_src #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_een_rescaled_quad_e_gl ( const qmckl_context context, const int64_t num, const int64_t* indices, 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* quad_ee_distance, const double* een_rescaled_quad_e, double* const een_rescaled_quad_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_quad_e_gl_doc ( const qmckl_context context, const int64_t num, const int64_t* indices, 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* quad_ee_distance, const double* een_rescaled_quad_e, double* const een_rescaled_quad_e_gl ); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes qmckl_exit_code qmckl_compute_een_rescaled_quad_e_gl ( const qmckl_context context, const int64_t num, const int64_t* indices, 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* quad_ee_distance, const double* een_rescaled_quad_e, double* const een_rescaled_quad_e_gl ) { #ifdef HAVE_HPC return qmckl_compute_een_rescaled_quad_e_gl_doc #else return qmckl_compute_een_rescaled_quad_e_gl_doc #endif (context, num, indices, walk_num, elec_num, cord_num, rescale_factor_ee, coord, coord_ee, quad_ee_distance, een_rescaled_quad_e, een_rescaled_quad_e_gl ); } #+end_src *** Test :noexport: #+begin_src c :tangle (eval c_test) printf("Een rescaled quad e gl\n"); /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); double een_rescaled_ee_gl[walk_num][cord_num+1][elec_num][4][elec_num]; double een_rescaled_quad_e_gl[num][cord_num+1][elec_num][3]; for (int elec = 0; elec < num; elec++){ rc = qmckl_set_electron_coord(context, 'N', walk_num, elec_coord, walk_num*elec_num*3); 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); rc = qmckl_get_electron_coord(context, 'N', &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); rc = qmckl_set_quad_points(context, 'N', num, indices, &new_coords[0][0], 3*num); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_een_rescaled_quad_e_gl(context, &een_rescaled_quad_e_gl[0][0][0][0], num*(cord_num+1)*elec_num*3); assert (rc == QMCKL_SUCCESS); coords[0][indices[elec]][0] = new_coords[elec][0]; coords[0][indices[elec]][1] = new_coords[elec][1]; coords[0][indices[elec]][2] = new_coords[elec][2]; rc = qmckl_set_electron_coord(context, 'N', walk_num, &coords[0][0][0], walk_num*elec_num*3); 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 i = 0; i < elec_num; i++) { for (int m = 0; m < 3; m++) { //printf("een_rescaled_ee_gl[0][l][i][m][indices[elec]] %i %i %i %f \n", l, m ,i, een_rescaled_ee_gl[0][l][i][m][indices[elec]]); //printf("een_rescaled_ee_gl[0][l][indices[elec]][m][i] %i %i %i %f \n", l, m ,i, een_rescaled_ee_gl[0][l][indices[elec]][m][i]); //printf("een_rescaled_quad_e_gl[elec][l][i][m] %i %i %i %f\n", l, m, i,een_rescaled_quad_e_gl[elec][l][i][m]); assert(fabs(een_rescaled_ee_gl[0][l][i][m][indices[elec]] - een_rescaled_quad_e_gl[elec][l][i][m]) < 1.e-12); assert(fabs(een_rescaled_ee_gl[0][l][indices[elec]][m][i] + een_rescaled_quad_e_gl[elec][l][i][m]) < 1.e-12); } } } } printf("OK\n"); #+end_src ** $\delta P$ matrix gradients and Laplacian *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_jastrow_champ_quad_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_quad_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_quad_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->quad_point.num * ctx->jastrow_champ.cord_num * (ctx->jastrow_champ.cord_num + 1) * ctx->nucleus.num * 3; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_jastrow_champ_quad_delta_p_gl", "Array too small."); } memcpy(delta_p_gl, ctx->quad_point.delta_p_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_jastrow_champ_quad_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_quad_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_quad_e(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_een_rescaled_n(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_een_rescaled_e(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_een_rescaled_n_gl(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_een_rescaled_e_gl(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_een_rescaled_quad_n(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_een_rescaled_quad_e_gl(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_een_rescaled_quad_n_gl(context); if(rc != QMCKL_SUCCESS) return rc; /* Compute if necessary */ if (ctx->quad_point.date > ctx->quad_point.delta_p_gl_date) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = 3 * ctx->quad_point.num * ctx->jastrow_champ.cord_num * (ctx->jastrow_champ.cord_num + 1) * ctx->nucleus.num * sizeof(double); if (mem_info.size > ctx->quad_point.delta_p_gl_maxsize) { if (ctx->quad_point.delta_p_gl != NULL) { rc = qmckl_free(context, ctx->quad_point.delta_p_gl); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_jastrow_champ_delta_p_gl", "Unable to free ctx->quad_point.delta_p_gl"); } ctx->quad_point.delta_p_gl = NULL; } } /* Allocate array */ if (ctx->quad_point.delta_p_gl == NULL) { double* delta_p_gl = (double*) qmckl_malloc(context, mem_info); ctx->quad_point.delta_p_gl_maxsize = mem_info.size; if (delta_p_gl == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_jastrow_champ_delta_p_gl", NULL); } ctx->quad_point.delta_p_gl = delta_p_gl; } rc = qmckl_compute_jastrow_champ_quad_delta_p_gl(context, ctx->quad_point.num, ctx->quad_point.indices, 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->quad_point.een_rescaled_quad_n, ctx->quad_point.een_rescaled_quad_e, ctx->jastrow_champ.een_rescaled_n_gl, ctx->jastrow_champ.een_rescaled_e_gl, ctx->quad_point.een_rescaled_quad_n_gl, ctx->quad_point.een_rescaled_quad_e_gl, ctx->quad_point.delta_p_gl); if (rc != QMCKL_SUCCESS) { return rc; } ctx->quad_point.delta_p_gl_date = ctx->quad_point.date; } return QMCKL_SUCCESS; } #+end_src *** Compute :PROPERTIES: :Name: qmckl_compute_jastrow_champ_quad_delta_p_gl_doc :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_factor_quad_delta_p_gl_args | Variable | Type | In/Out | Description | |----------------------------+---------------------------------------------------------------------+--------+---------------------------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~num~ | ~int64_t~ | in | Index of quad electron | | ~indices~ | ~int64_t[num]~ | in | Indices of quad points | | ~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_quad_n~ | ~double[num][0:cord_num][nucl_num]~ | in | Electron-nucleus quad rescaled distances | | ~een_rescaled_quad_e~ | ~double[num][0:cord_num][elec_num]~ | in | Electron-electron quad rescaled distances | | ~een_rescaled_n_gl~ | ~double[walk_num][0:cord_num][nucl_num][4][elec_num]~ | in | Electron-nucleus rescaled distances derivatives | | ~een_rescaled_e_gl~ | ~double[walk_num][0:cord_num][elec_num][4][elec_num]~ | in | Electron-electron rescaled distances derivatives | | ~een_rescaled_quad_n_gl~ | ~double[num][0:cord_num][nucl_num][3]~ | in | Electron-nucleus quad rescaled distances derivatives | | ~een_rescaled_quad_e_gl~ | ~double[num][0:cord_num][elec_num][3]~ | in | Electron-electron quad rescaled distances derivatives | | ~delta_p_gl~ | ~double[num][0:cord_num-1][0:cord_num][3][nucl_num]~ | out | Delta P matrix gradient and Laplacian | |----------------------------+---------------------------------------------------------------------+--------+---------------------------------------------------------| #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer(qmckl_exit_code) function qmckl_compute_jastrow_champ_quad_delta_p_gl_doc( & context, num, indices, walk_num, elec_num, nucl_num, cord_num, & een_rescaled_n, een_rescaled_e, een_rescaled_quad_n, een_rescaled_quad_e, & een_rescaled_n_gl, een_rescaled_e_gl, een_rescaled_quad_n_gl, een_rescaled_quad_e_gl, delta_p_gl) & result(info) bind(C) use, intrinsic :: iso_c_binding use qmckl implicit none integer(qmckl_context), intent(in) :: context integer(c_int64_t) , intent(in), value :: num, walk_num, elec_num, cord_num, nucl_num integer(c_int64_t) , intent(in) :: indices(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_quad_n(nucl_num, 0:cord_num, num) real(c_double) , intent(in) :: een_rescaled_quad_e(elec_num, 0:cord_num, 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_quad_n_gl(3, nucl_num, 0:cord_num, num) real(c_double) , intent(in) :: een_rescaled_quad_e_gl(3,elec_num, 0:cord_num, num) real(c_double) , intent(out) :: delta_p_gl(nucl_num,3,0:cord_num, 0:cord_num-1, num) double precision :: delta_e_gl(elec_num,3) integer*8 :: i, a, j, l, k, p, m, n, nw, idx double precision :: tmp, accu integer*8 :: LDA, LDB, LDC 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) then delta_p_gl = 0.d0 return endif do nw=1, num idx = indices(nw)+1 do m=0, cord_num-1 do k = 1, 3 do j = 1, elec_num delta_e_gl(j,k) = een_rescaled_quad_e_gl(k,j,m,nw) - een_rescaled_e_gl(idx, k, j, m, 1) end do delta_e_gl(idx,k) = 0.0d0 end do do l=0, cord_num do k = 1, 3 do a = 1, nucl_num accu = 0.0d0 do j = 1, elec_num accu = accu + delta_e_gl(j,k) * een_rescaled_n(j,a,l,1) end do delta_p_gl(a,k,l,m,nw) = accu end do end do end do end do end do end function qmckl_compute_jastrow_champ_quad_delta_p_gl_doc #+end_src #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_jastrow_champ_quad_delta_p_gl_doc (const qmckl_context context, const int64_t num, const int64_t* indices, 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_quad_n, const double* een_rescaled_quad_e, const double* een_rescaled_n_gl, const double* een_rescaled_e_gl, const double* een_rescaled_quad_n_gl, const double* een_rescaled_quad_e_gl, double* const delta_p_gl ); qmckl_exit_code qmckl_compute_jastrow_champ_quad_delta_p_gl (const qmckl_context context, const int64_t num, const int64_t* indices, 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_quad_n, const double* een_rescaled_quad_e, const double* een_rescaled_n_gl, const double* een_rescaled_e_gl, const double* een_rescaled_quad_n_gl, const double* een_rescaled_quad_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_quad_delta_p_gl (const qmckl_context context, const int64_t num, const int64_t* indices, 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_quad_n, const double* een_rescaled_quad_e, const double* een_rescaled_n_gl, const double* een_rescaled_e_gl, const double* een_rescaled_quad_n_gl, const double* een_rescaled_quad_e_gl, double* const delta_p_gl ) { #ifdef HAVE_HPC return qmckl_compute_jastrow_champ_quad_delta_p_gl_doc #else return qmckl_compute_jastrow_champ_quad_delta_p_gl_doc #endif (context, num, indices, walk_num, elec_num, nucl_num, cord_num, een_rescaled_n, een_rescaled_e, een_rescaled_quad_n, een_rescaled_quad_e, een_rescaled_n_gl, een_rescaled_e_gl, een_rescaled_quad_n_gl, een_rescaled_quad_e_gl, delta_p_gl); } #+end_src *** Test :noexport: #+begin_src c :tangle (eval c_test) printf("Delta P gl\n"); /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); double p_gl_old[walk_num][cord_num][cord_num+1][nucl_num][4][elec_num]; double delta_p_gl[num][cord_num][cord_num+1][3][nucl_num]; double p_gl_new[walk_num][cord_num][cord_num+1][nucl_num][4][elec_num]; for (int elec = 0; elec < num; elec++){ rc = qmckl_set_electron_coord(context, 'N', walk_num, elec_coord, walk_num*elec_num*3); assert(rc == QMCKL_SUCCESS); rc = qmckl_get_jastrow_champ_dtmp_c(context, &p_gl_old[0][0][0][0][0][0], walk_num*cord_num*(cord_num+1)*nucl_num*4*elec_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_quad_points(context, 'N', num, indices, &new_coords[0][0], 3*num); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_jastrow_champ_quad_delta_p_gl(context, &delta_p_gl[0][0][0][0][0], 3*num*cord_num*(cord_num+1)*nucl_num); assert (rc == QMCKL_SUCCESS); coords[0][indices[elec]][0] = new_coords[elec][0]; coords[0][indices[elec]][1] = new_coords[elec][1]; coords[0][indices[elec]][2] = new_coords[elec][2]; rc = qmckl_set_electron_coord(context, 'N', walk_num, &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_jastrow_champ_dtmp_c(context, &p_gl_new[0][0][0][0][0][0], walk_num*cord_num*(cord_num+1)*nucl_num*4*elec_num); assert (rc == QMCKL_SUCCESS); 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 k = 0; k < 3; k++){ if (fabs(((p_gl_new[0][l][m][a][k][indices[elec]]-p_gl_old[0][l][m][a][k][indices[elec]])-delta_p_gl[elec][l][m][k][a])) > 1.e-12) { printf("p_gl[%d][%d][%d][%d][%d][%d] = %f\n", 0, l, m, a, k, indices[elec], p_gl_new[0][l][m][a][k][indices[elec]] - p_gl_old[0][l][m][a][k][indices[elec]]); printf("delta_p_gl[%d][%d][%d][%d][%d] = %f\n", elec, l, m, a, k, delta_p_gl[elec][l][m][k][a]); } assert(fabs(((p_gl_new[0][l][m][a][k][indices[elec]]-p_gl_old[0][l][m][a][k][indices[elec]])-delta_p_gl[elec][l][m][k][a])) < 1.e-12); } } } } } printf("OK\n"); #+end_src ** Electron-electron-nucleus Jastrow gradients and Laplacian *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_jastrow_champ_quad_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_quad_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_quad_een_gl(context); if (rc != QMCKL_SUCCESS) return rc; qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); if (delta_een_gl == NULL) { return qmckl_failwith( context, QMCKL_INVALID_ARG_2, "qmckl_get_jastrow_champ_quad_een_gl", "Array is NULL."); } int64_t sze = 3 * ctx->quad_point.num; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_jastrow_champ_quad_een_gl", "Array too small. Expected 3 * ctx->quad_point.num"); } memcpy(delta_een_gl, ctx->quad_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_quad_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_quad_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_quad_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; if (ctx->jastrow_champ.cord_num > 0) { rc = qmckl_provide_jastrow_champ_quad_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_quad_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->quad_point.date > ctx->quad_point.delta_een_gl_date) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = 3 * ctx->quad_point.num * sizeof(double); if (mem_info.size > ctx->quad_point.delta_een_gl_maxsize) { if (ctx->quad_point.delta_een_gl != NULL) { rc = qmckl_free(context, ctx->quad_point.delta_een_gl); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_jastrow_champ_quad_een_gl", "Unable to free ctx->quad_point.delta_een_gl"); } ctx->quad_point.delta_een_gl = NULL; } } /* Allocate array */ if (ctx->quad_point.delta_een_gl == NULL) { double* delta_een_gl = (double*) qmckl_malloc(context, mem_info); ctx->quad_point.delta_een_gl_maxsize = mem_info.size; if (delta_een_gl == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_jastrow_champ_quad_een_gl", NULL); } ctx->quad_point.delta_een_gl = delta_een_gl; } rc = qmckl_compute_jastrow_champ_factor_quad_een_gl_doc(context, ctx->quad_point.num, ctx->quad_point.indices, 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->quad_point.delta_p, ctx->quad_point.delta_p_gl, ctx->jastrow_champ.een_rescaled_n, ctx->quad_point.een_rescaled_quad_n, ctx->jastrow_champ.een_rescaled_n_gl, ctx->quad_point.een_rescaled_quad_n_gl, ctx->quad_point.delta_een_gl); if (rc != QMCKL_SUCCESS) { return rc; } ctx->quad_point.delta_een_gl_date = ctx->quad_point.date; } return QMCKL_SUCCESS; } #+end_src *** Compute :PROPERTIES: :Name: qmckl_compute_jastrow_champ_factor_quad_een_gl_doc :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_factor_quad_een_gl_args | Variable | Type | In/Out | Description | |----------------------------+---------------------------------------------------------------------+--------+----------------------------------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~num~ | ~int64_t~ | in | Index of quad electron | | ~indices~ | ~int64_t[num]~ | in | Indices of quad points | | ~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 | P matrix | | ~dtmp_c~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][4][elec_num]~ | in | P matrix derivative | | ~delta_p~ | ~double[num][0:cord_num-1][0:cord_num][nucl_num][elec_num]~ | in | Delta P matrix | | ~delta_p_gl~ | ~double[num][0:cord_num-1][0:cord_num][3][nucl_num]~ | in | Delta P matrix derivative | | ~een_rescaled_n~ | ~double[walk_num][0:cord_num][nucl_num][elec_num]~ | in | Electron-nucleus rescaled distances | | ~een_rescaled_quad_n~ | ~double[num][0:cord_num][nucl_num]~ | in | Electron-nucleus quad rescaled distances | | ~een_rescaled_n_gl~ | ~double[walk_num][0:cord_num][nucl_num][4][elec_num]~ | in | Electron-nucleus rescaled distances derivatives | | ~een_rescaled_quad_n_gl~ | ~double[num][0:cord_num][nucl_num][3]~ | in | Electron-nucleus quad rescaled distances derivatives | | ~delta_een_gl~ | ~double[num][3][elec_num]~ | out | Delta electron-electron-nucleus jastrow gradient and Laplacian | |----------------------------+---------------------------------------------------------------------+--------+----------------------------------------------------------------| #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer(qmckl_exit_code) function qmckl_compute_jastrow_champ_factor_quad_een_gl_doc( & context, num, indices, 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_quad_n, & een_rescaled_n_gl, een_rescaled_quad_n_gl, delta_een_gl) & result(info) bind(C) use, intrinsic :: iso_c_binding use qmckl implicit none integer(qmckl_context), intent(in) :: context integer(c_int64_t) , intent(in), value :: num, walk_num, elec_num, cord_num, nucl_num, dim_c_vector integer(c_int64_t) , intent(in) :: indices(num) integer(c_int64_t) , intent(in) :: lkpm_combined_index(dim_c_vector,4) real(c_double) , intent(in) :: c_vector_full(nucl_num, dim_c_vector) 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, num) real(c_double) , intent(in) :: delta_p_gl(nucl_num, 3, 0:cord_num, 0:cord_num-1, 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_quad_n(nucl_num, 0:cord_num, 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_quad_n_gl(3, nucl_num, 0:cord_num, num) real(c_double) , intent(out) :: delta_een_gl(3, num) integer*8 :: i, a, j, l, k, p, m, n, nw, kk, idx double precision :: accu, accu2, cn integer*8 :: LDA, LDB, LDC double precision :: een_rescaled_delta_n_gl(3, nucl_num, 0:cord_num) double precision :: een_rescaled_delta_n(nucl_num, 0:cord_num) 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 if (cord_num == 0) return do nw =1, num idx = indices(nw)+1 een_rescaled_delta_n(:,:) = een_rescaled_quad_n(:,:,nw) - een_rescaled_n(idx, :, :, 1) een_rescaled_delta_n_gl(:,:,:) = een_rescaled_quad_n_gl(:,:,:,nw) - een_rescaled_n_gl(idx, :,:,:,1) 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 kk = 1, 3 do a = 1, nucl_num cn = c_vector_full(a, n) if(cn == 0.d0) cycle delta_een_gl(kk,nw) = delta_een_gl(kk,nw) + ( & delta_p_gl(a,kk,m ,k,nw) * een_rescaled_n(idx,a,m+l,1) + & delta_p_gl(a,kk,m+l,k,nw) * een_rescaled_n(idx,a,m ,1) + & delta_p(idx,a,m ,k,nw) * een_rescaled_n_gl(idx,kk,a,m+l,1) + & delta_p(idx,a,m+l,k,nw) * een_rescaled_n_gl(idx,kk,a,m ,1) ) * cn delta_een_gl(kk,nw) = delta_een_gl(kk,nw) + ( & (dtmp_c(idx,kk,a,m ,k,1) + delta_p_gl(a,kk,m ,k,nw)) * een_rescaled_delta_n(a,m+l) + & (dtmp_c(idx,kk,a,m+l,k,1) + delta_p_gl(a,kk,m+l,k,nw)) * een_rescaled_delta_n(a,m) + & (tmp_c(idx,a,m ,k,1) + delta_p(idx,a,m ,k,nw)) * een_rescaled_delta_n_gl(kk,a,m+l) + & (tmp_c(idx,a,m+l,k,1) + delta_p(indices(nw)+1,a,m+l,k,nw)) * een_rescaled_delta_n_gl(kk,a,m) )* cn end do end do end do end do end function qmckl_compute_jastrow_champ_factor_quad_een_gl_doc #+end_src #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_jastrow_champ_factor_quad_een_gl_doc (const qmckl_context context, const int64_t num, const int64_t* indices, 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_quad_n, const double* een_rescaled_n_gl, const double* een_rescaled_quad_n_gl, double* const delta_een_gl ); qmckl_exit_code qmckl_compute_jastrow_champ_factor_quad_een_gl (const qmckl_context context, const int64_t num, const int64_t* indices, 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_quad_n, const double* een_rescaled_n_gl, const double* een_rescaled_quad_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_quad_een_gl (const qmckl_context context, const int64_t num, const int64_t* indices, 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_quad_n, const double* een_rescaled_n_gl, const double* een_rescaled_quad_n_gl, double* const delta_een_gl ) { #ifdef HAVE_HPC return qmckl_compute_jastrow_champ_factor_quad_een_gl_doc #else return qmckl_compute_jastrow_champ_factor_quad_een_gl_doc #endif (context, num, indices, 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_quad_n, een_rescaled_n_gl, een_rescaled_quad_n_gl, delta_een_gl ); } #+end_src *** Test :noexport: #+begin_src c :tangle (eval c_test) printf("Delta een gl\n"); /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); double een_gl_old[walk_num][4][elec_num]; double delta_een_gl[num][3]; double een_gl_new[walk_num][4][elec_num]; for (int elec = 0; elec < num; elec++){ rc = qmckl_set_electron_coord(context, 'N', walk_num, elec_coord, walk_num*elec_num*3); assert(rc == QMCKL_SUCCESS); rc = qmckl_get_jastrow_champ_factor_een_gl(context, &een_gl_old[0][0][0], walk_num*4*elec_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_quad_points(context, 'N', num, indices, &new_coords[0][0], 3*num); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_jastrow_champ_quad_een_gl(context, &delta_een_gl[0][0], num*3); assert (rc == QMCKL_SUCCESS); coords[0][indices[elec]][0] = new_coords[elec][0]; coords[0][indices[elec]][1] = new_coords[elec][1]; coords[0][indices[elec]][2] = new_coords[elec][2]; rc = qmckl_set_electron_coord(context, 'N', walk_num, &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_jastrow_champ_factor_een_gl(context, &een_gl_new[0][0][0], walk_num*4*elec_num); assert (rc == QMCKL_SUCCESS); for (int m = 0; m < 3; m++) { printf("delta_een_gl[%d][%d] = %f\n", elec, m, delta_een_gl[elec][m]); printf("een_gl_[%d][%d][%d] = %f\n", 0, m, elec, een_gl_new[0][m][indices[elec]]-een_gl_old[0][m][indices[elec]]); assert(fabs((een_gl_new[0][m][indices[elec]]- een_gl_old[0][m][indices[elec]]) - delta_een_gl[elec][m]) < 1.e-12); } } printf("OK\n"); #+end_src * Electron-electron Jastrow ** Electron-electron rescaled distance *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_ee_rescaled_quad(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_quad(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_quad(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->quad_point.num; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "todo", "Array too small. Expected ctx->electron.num * ctx->quad_point.num "); } memcpy(distance_rescaled, ctx->quad_point.ee_rescaled_quad, 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_quad(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_provide_ee_rescaled_quad(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_quad_ee_distance(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_ee_distance(context); if(rc != QMCKL_SUCCESS) return rc; /* Compute if necessary */ if (ctx->quad_point.date > ctx->quad_point.ee_rescaled_quad_date) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = ctx->electron.num * ctx->quad_point.num * sizeof(double); if (mem_info.size > ctx->quad_point.ee_rescaled_quad_maxsize) { if (ctx->quad_point.ee_rescaled_quad!= NULL) { rc = qmckl_free(context, ctx->quad_point.ee_rescaled_quad); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_ee_rescaled_quad", "Unable to free ctx->quad_point.ee_rescaled_quad"); } ctx->quad_point.ee_rescaled_quad = NULL; } } /* Allocate array */ if (ctx->quad_point.ee_rescaled_quad == NULL) { double* ee_rescaled_quad = (double*) qmckl_malloc(context, mem_info); ctx->quad_point.ee_rescaled_quad_maxsize = mem_info.size; if (ee_rescaled_quad == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_ee_rescaled_quad", NULL); } ctx->quad_point.ee_rescaled_quad = ee_rescaled_quad; } rc = qmckl_compute_ee_rescaled_quad(context, ctx->quad_point.num, ctx->electron.num, ctx->jastrow_champ.rescale_factor_ee, ctx->electron.walker.num, ctx->quad_point.quad_ee_distance, ctx->quad_point.ee_rescaled_quad); if (rc != QMCKL_SUCCESS) { return rc; } ctx->quad_point.ee_rescaled_quad_date = ctx->quad_point.date; } return QMCKL_SUCCESS; } #+end_src *** Compute :PROPERTIES: :Name: qmckl_compute_ee_rescaled_quad :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_ee_rescaled_quad_args | Variable | Type | In/Out | Description | |----------------------+------------------------------+--------+--------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~num~ | ~int64_t~ | in | Index of quad electron | | ~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 | | ~quad_ee_distance~ | ~double[num][elec_num]~ | in | quad electron-electron distances | | ~ee_rescaled_quad~ | ~double[num][elec_num]~ | out | Electron-electron rescaled distances | |----------------------+------------------------------+--------+--------------------------------------| #+begin_src f90 :comments org :tangle (eval f) :noweb yes function qmckl_compute_ee_rescaled_quad_doc(context, num,& elec_num, rescale_factor_ee, walk_num, & quad_ee_distance, ee_rescaled_quad) & bind(C) result(info) use qmckl implicit none integer(qmckl_context), intent(in), value :: context integer(c_int64_t) , intent(in) , value :: num 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) :: quad_ee_distance(elec_num,num) real (c_double ) , intent(out) :: ee_rescaled_quad(elec_num,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,num do i=1,elec_num ee_rescaled_quad(i,k) = (1.0d0 - dexp(-rescale_factor_ee * quad_ee_distance(i,k))) * inverse_rescale_factor_ee enddo end do end function qmckl_compute_ee_rescaled_quad_doc #+end_src #+begin_src c :tangle (eval h_private_func) :comments org :exports none qmckl_exit_code qmckl_compute_ee_rescaled_quad_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* quad_ee_distance, double* const ee_rescaled_quad ); qmckl_exit_code qmckl_compute_ee_rescaled_quad( 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* quad_ee_distance, double* const ee_rescaled_quad ); #+end_src #+begin_src c :tangle (eval c) :comments org :exports none qmckl_exit_code qmckl_compute_ee_rescaled_quad ( 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* quad_ee_distance, double* const ee_rescaled_quad ) { #ifdef HAVE_HPC return qmckl_compute_ee_rescaled_quad_doc #else return qmckl_compute_ee_rescaled_quad_doc #endif (context, num, elec_num, rescale_factor_ee, walk_num, quad_ee_distance, ee_rescaled_quad); } #+end_src *** Test :noexport: #+begin_src c :tangle (eval c_test) printf("ee rescaled quad\n"); /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); double ee_rescaled[walk_num][elec_num][elec_num]; double quad_ee_rescaled[num][elec_num]; for (int elec = 0; elec < num; elec++){ rc = qmckl_set_electron_coord(context, 'N', walk_num, elec_coord, walk_num*elec_num*3); assert(rc == QMCKL_SUCCESS); rc = qmckl_get_jastrow_champ_ee_distance_rescaled(context, &ee_rescaled[0][0][0], walk_num*elec_num*elec_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_quad_points(context, 'N', num, indices, &new_coords[0][0], 3*num); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_ee_rescaled_quad(context, &quad_ee_rescaled[0][0], num*elec_num); assert (rc == QMCKL_SUCCESS); coords[0][indices[elec]][0] = new_coords[elec][0]; coords[0][indices[elec]][1] = new_coords[elec][1]; coords[0][indices[elec]][2] = new_coords[elec][2]; rc = qmckl_set_electron_coord(context, 'N', walk_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], walk_num*elec_num*elec_num); assert (rc == QMCKL_SUCCESS); for (int i = 0; i < elec_num; i++){ //printf("nw %d i %d %f %f\n", nw, i, ee_rescaled[nw][2][i], quad_ee_rescaled[nw][i]); assert(fabs(ee_rescaled[0][indices[elec]][i]-quad_ee_rescaled[elec][i]) < 1.e-12); assert(fabs(ee_rescaled[0][i][indices[elec]]-quad_ee_rescaled[elec][i]) < 1.e-12); } } printf("OK\n"); #+end_src ** Electron-electron Jastrow value *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_jastrow_champ_quad_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_quad_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_quad_ee", NULL); } qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); rc = qmckl_provide_jastrow_champ_quad_ee(context); if (rc != QMCKL_SUCCESS) return rc; int64_t sze=ctx->quad_point.num; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_jastrow_champ_quad_ee", "Array too small. Expected quad_point.num"); } memcpy(delta_ee, ctx->quad_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_quad_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_quad_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_quad_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_quad_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_quad_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_quad_ee", NULL); } rc = qmckl_provide_ee_distance_rescaled(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_ee_rescaled_quad(context); if(rc != QMCKL_SUCCESS) return rc; /* Compute if necessary */ if (ctx->quad_point.date > ctx->quad_point.delta_ee_date) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = ctx->quad_point.num * sizeof(double); if (mem_info.size > ctx->quad_point.delta_ee_maxsize) { if (ctx->quad_point.delta_ee != NULL) { rc = qmckl_free(context, ctx->quad_point.delta_ee); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_jastrow_champ_quad_ee", "Unable to free ctx->quad_point.delta_ee"); } ctx->quad_point.delta_ee = NULL; } } /* Allocate array */ if (ctx->quad_point.delta_ee == NULL) { double* delta_ee = (double*) qmckl_malloc(context, mem_info); ctx->quad_point.delta_ee_maxsize = mem_info.size; if (delta_ee == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_jastrow_champ_quad_ee", NULL); } ctx->quad_point.delta_ee = delta_ee; } rc = qmckl_compute_jastrow_champ_quad_ee(context, ctx->quad_point.num, ctx->quad_point.indices, 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->quad_point.ee_rescaled_quad, ctx->jastrow_champ.spin_independent, ctx->quad_point.delta_ee); if (rc != QMCKL_SUCCESS) { return rc; } ctx->quad_point.delta_ee_date = ctx->quad_point.date; } return QMCKL_SUCCESS; } #+end_src *** Compute :PROPERTIES: :Name: qmckl_compute_jastrow_champ_quad_ee :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_quad_ee_args | Variable | Type | In/Out | Description | |------------------------+----------------------------------------+--------+---------------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~num~ | ~int64_t~ | in | Index of quad point | | ~num~ | ~int64_t~ | in | Index of quad point | | ~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 rescaled distances | | ~ee_rescaled_quad~ | ~double[num][elec_num]~ | in | Electron-electron rescaled quad distances | | ~delta_ee~ | ~double[num]~ | out | quad electron-electron Jastrow | |------------------------+----------------------------------------+--------+---------------------------------------------| #+begin_src f90 :comments org :tangle (eval f) :noweb yes function qmckl_compute_jastrow_champ_quad_ee_doc(context, & num, indices, walk_num, elec_num, up_num, bord_num, b_vector, & ee_distance_rescaled, ee_rescaled_quad, 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 integer (c_int64_t) , intent(in) :: indices(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 :: 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_quad(elec_num,num) integer (c_int32_t) , intent(in), value :: spin_independent real (c_double ) , intent(out) :: delta_ee(num) integer(qmckl_exit_code) :: info integer*8 :: i, j, k, nw double precision :: x, xk, y, yk 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, num delta_ee(nw) = 0.0d0 do i=1,elec_num if (i.ne.(indices(nw)+1)) then x = ee_distance_rescaled(i,indices(nw)+1,1) y = ee_rescaled_quad(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 .and. (indices(nw)+1) <= up_num ) .or. (i > up_num .and. (indices(nw)+1) > 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_quad_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_quad_ee_doc (const qmckl_context context, const int64_t num, const int64_t* indices, 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 #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_jastrow_champ_quad_ee (const qmckl_context context, const int64_t num, const int64_t* indices, 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_quad, 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_quad_ee (const qmckl_context context, const int64_t num, const int64_t* indices, 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_quad, const int32_t spin_independent, double* const delta_ee ) { #ifdef HAVE_HPC return qmckl_compute_jastrow_champ_quad_ee_doc #else return qmckl_compute_jastrow_champ_quad_ee_doc #endif (context, num, indices, walk_num, elec_num, up_num, bord_num, b_vector, ee_distance_rescaled, ee_rescaled_quad, spin_independent, delta_ee); } #+end_src *** Test :noexport: #+begin_src c :tangle (eval c_test) printf("Delta ee\n"); /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); double jastrow_ee_old[walk_num]; double delta_ee[num]; double jastrow_ee_new[walk_num]; for (int elec = 0; elec < num; elec++){ rc = qmckl_set_electron_coord(context, 'N', walk_num, elec_coord, walk_num*elec_num*3); assert(rc == QMCKL_SUCCESS); 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_quad_points(context, 'N', num, indices, &new_coords[0][0], 3*num); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_jastrow_champ_quad_ee(context, &delta_ee[0], num); assert (rc == QMCKL_SUCCESS); coords[0][indices[elec]][0] = new_coords[elec][0]; coords[0][indices[elec]][1] = new_coords[elec][1]; coords[0][indices[elec]][2] = new_coords[elec][2]; rc = qmckl_set_electron_coord(context, 'N', walk_num, &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_jastrow_champ_factor_ee(context, &jastrow_ee_new[0], walk_num); assert (rc == QMCKL_SUCCESS); //printf("%f %f %f %3.14f\n", jastrow_ee_new[0], jastrow_ee_old[0], delta_ee[elec], fabs((jastrow_ee_new[0] - jastrow_ee_old[0]) - delta_ee[elec])); assert(fabs((jastrow_ee_new[0] - jastrow_ee_old[0]) - delta_ee[elec]) < 1.e-12); } printf("OK\n"); #+end_src ** Electron-electron rescaled distances derivatives *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_ee_rescaled_quad_gl(qmckl_context context, double* const distance_rescaled_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_ee_rescaled_quad_gl(qmckl_context context, double* const distance_rescaled_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_ee_rescaled_quad_gl(context); if (rc != QMCKL_SUCCESS) return rc; qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); if (distance_rescaled_gl == NULL) { return qmckl_failwith( context, QMCKL_INVALID_ARG_2, "qmckl_get_ee_rescaled_quad_gl", "Array is NULL"); } int64_t sze = 3 * ctx->electron.num * ctx->quad_point.num; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_ee_rescaled_quad_gl", "Array too small. Expected 3 * ctx->electron.num * ctx->quad_point.num"); } memcpy(distance_rescaled_gl, ctx->quad_point.ee_rescaled_quad_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_quad_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_quad_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_quad_ee_distance(context); if(rc != QMCKL_SUCCESS) return rc; /* Compute if necessary */ if (ctx->quad_point.date > ctx->quad_point.ee_rescaled_quad_gl_date) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = 3 * ctx->electron.num * ctx->quad_point.num * sizeof(double); if (mem_info.size > ctx->quad_point.ee_rescaled_quad_gl_maxsize) { if (ctx->quad_point.ee_rescaled_quad_gl != NULL) { rc = qmckl_free(context, ctx->quad_point.ee_rescaled_quad_gl); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_ee_rescaled_quad_gl", "Unable to free ctx->quad_point.ee_rescaled_quad_gl"); } ctx->quad_point.ee_rescaled_quad_gl = NULL; } } /* Allocate array */ if (ctx->quad_point.ee_rescaled_quad_gl == NULL) { double* ee_rescaled_quad_gl = (double*) qmckl_malloc(context, mem_info); ctx->quad_point.ee_rescaled_quad_gl_maxsize = mem_info.size; if (ee_rescaled_quad_gl == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_ee_rescaled_quad_gl", NULL); } ctx->quad_point.ee_rescaled_quad_gl = ee_rescaled_quad_gl; } qmckl_exit_code rc = qmckl_compute_ee_rescaled_quad_gl(context, ctx->quad_point.num, ctx->quad_point.indices, ctx->electron.num, ctx->jastrow_champ.rescale_factor_ee, ctx->electron.walker.num, ctx->quad_point.quad_ee_distance, ctx->electron.walker.point.coord.data, ctx->quad_point.coord.data, ctx->quad_point.ee_rescaled_quad_gl); if (rc != QMCKL_SUCCESS) { return rc; } ctx->quad_point.ee_rescaled_quad_gl_date = ctx->quad_point.date; } return QMCKL_SUCCESS; } #+end_src *** Compute :PROPERTIES: :Name: qmckl_compute_ee_rescaled_quad_gl :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_ee_rescaled_quad_gl_args | Variable | Type | In/Out | Description | |-------------------------+---------------------------------+--------+--------------------------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~num~ | ~int64_t~ | in | Index of quad electron | | ~indices~ | ~int64_t[num]~ | in | Index of quad electron | | ~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 | | ~quad_ee_distance~ | ~double[elec_num][num]~ | in | quad electron-electron distances | | ~elec_coord~ | ~double[3][walk_num][elec_num]~ | in | Electron coordinates | | ~coord~ | ~double[num][3]~ | in | quad electron coordinates | | ~ee_rescaled_quad_gl~ | ~double[num][elec_num][3]~ | out | Electron-electron rescaled quad distance derivatives | |-------------------------+---------------------------------+--------+--------------------------------------------------------| #+begin_src f90 :comments org :tangle (eval f) :noweb yes function qmckl_compute_ee_rescaled_quad_gl_doc(context, num, indices, & elec_num, rescale_factor_ee, walk_num, quad_ee_distance, elec_coord, coord, ee_rescaled_quad_gl) & bind(C) result(info) use qmckl implicit none integer(qmckl_context), intent(in), value :: context integer (c_int64_t) , intent(in) , value :: num integer (c_int64_t) , intent(in) :: indices(num) 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) :: quad_ee_distance(elec_num,num) real (c_double ) , intent(in) :: elec_coord(elec_num,walk_num,3) real (c_double ) , intent(in) :: coord(3,num) real (c_double ) , intent(out) :: ee_rescaled_quad_gl(3,elec_num,num) integer(qmckl_exit_code) :: info integer*8 :: nw, i, ii double precision :: rij_inv, elel_dist_gl(3, elec_num), kappa_l 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_quad_gl = 0.0d0 do nw = 1, num ! prepare the actual een table do i = 1, elec_num rij_inv = 1.0d0 / quad_ee_distance(i, nw) do ii = 1, 3 elel_dist_gl(ii, i) = (elec_coord(i,1, ii) - coord(ii,nw)) * rij_inv end do end do do i = 1, elec_num kappa_l = -1 * rescale_factor_ee ee_rescaled_quad_gl(1, i, nw) = elel_dist_gl(1, i) ee_rescaled_quad_gl(2, i, nw) = elel_dist_gl(2, i) ee_rescaled_quad_gl(3, i, nw) = elel_dist_gl(3, i) ee_rescaled_quad_gl(4, i, nw) = elel_dist_gl(4, i) ee_rescaled_quad_gl(1, i, nw) = ee_rescaled_quad_gl(1, i, nw) * dexp(kappa_l * quad_ee_distance(i,nw)) ee_rescaled_quad_gl(2, i, nw) = ee_rescaled_quad_gl(2, i, nw) * dexp(kappa_l * quad_ee_distance(i,nw)) ee_rescaled_quad_gl(3, i, nw) = ee_rescaled_quad_gl(3, i, nw) * dexp(kappa_l * quad_ee_distance(i,nw)) end do ee_rescaled_quad_gl(1, indices(nw)+1, nw) = 0.0d0 ee_rescaled_quad_gl(2, indices(nw)+1, nw) = 0.0d0 ee_rescaled_quad_gl(3, indices(nw)+1, nw) = 0.0d0 end do end function qmckl_compute_ee_rescaled_quad_gl_doc #+end_src #+begin_src c :tangle (eval h_private_func) :comments org :exports none qmckl_exit_code qmckl_compute_ee_rescaled_quad_gl_doc ( const qmckl_context context, const int64_t num, const int64_t* indices, const int64_t elec_num, const double rescale_factor_ee, const int64_t walk_num, const double* quad_ee_distance, const double* elec_coord, const double* coord, double* const ee_rescaled_quad_gl ); qmckl_exit_code qmckl_compute_ee_rescaled_quad_gl ( const qmckl_context context, const int64_t num, const int64_t* indices, const int64_t elec_num, const double rescale_factor_ee, const int64_t walk_num, const double* quad_ee_distance, const double* elec_coord, const double* coord, double* const ee_rescaled_quad_gl ); #+end_src #+begin_src c :tangle (eval c) :comments org :exports none qmckl_exit_code qmckl_compute_ee_rescaled_quad_gl ( const qmckl_context context, const int64_t num, const int64_t* indices, const int64_t elec_num, const double rescale_factor_ee, const int64_t walk_num, const double* quad_ee_distance, const double* elec_coord, const double* coord, double* const ee_rescaled_quad_gl ) { #ifdef HAVE_HPC return qmckl_compute_ee_rescaled_quad_gl_doc #else return qmckl_compute_ee_rescaled_quad_gl_doc #endif (context, num, indices, elec_num, rescale_factor_ee, walk_num, quad_ee_distance, elec_coord, coord, ee_rescaled_quad_gl); } #+end_src *** Test :noexport: #+begin_src c :tangle (eval c_test) printf("ee rescaled quad gl\n"); /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); double ee_rescaled_gl[walk_num][elec_num][elec_num][4]; double quad_ee_rescaled_gl[num][elec_num][3]; for (int elec = 0; elec < num; elec++){ rc = qmckl_set_electron_coord(context, 'N', walk_num, elec_coord, 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], walk_num*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_quad_points(context, 'N', num, indices, &new_coords[0][0], 3*num); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_ee_rescaled_quad_gl(context, &quad_ee_rescaled_gl[0][0][0], num*elec_num*3); assert (rc == QMCKL_SUCCESS); coords[0][indices[elec]][0] = new_coords[elec][0]; coords[0][indices[elec]][1] = new_coords[elec][1]; coords[0][indices[elec]][2] = new_coords[elec][2]; rc = qmckl_set_electron_coord(context, 'N', walk_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], walk_num*elec_num*elec_num*4); assert (rc == QMCKL_SUCCESS); for (int i = 0; i < elec_num; i++) { for (int m = 0; m < 3; m++) { if (i == indices[elec]) continue; //printf("%f\n", ee_rescaled_gl[0][indices[elec]][i][m]); //printf("%f\n", quad_ee_rescaled_gl[elec][i][m]); assert(fabs(ee_rescaled_gl[0][indices[elec]][i][m] - quad_ee_rescaled_gl[elec][i][m]) < 1.e-12); assert(fabs(ee_rescaled_gl[0][i][indices[elec]][m] + quad_ee_rescaled_gl[elec][i][m]) < 1.e-12); } } } printf("OK\n"); #+end_src ** Electron-electron Jastrow gradients and Laplacian *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_jastrow_champ_quad_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_quad_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_quad_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->quad_point.num * 3 * ctx->electron.num; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_jastrow_champ_quad_ee_gl", "Array too small. Expected 3*num*elec_num"); } memcpy(delta_ee_gl, ctx->quad_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_quad_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_quad_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_quad_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_quad_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_quad_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_quad_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_quad(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_ee_rescaled_quad_gl(context); if(rc != QMCKL_SUCCESS) return rc; /* Compute if necessary */ if (ctx->quad_point.date > ctx->quad_point.delta_ee_gl_date) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = ctx->quad_point.num * 3 * ctx->electron.num * sizeof(double); if (mem_info.size > ctx->quad_point.delta_ee_gl_maxsize) { if (ctx->quad_point.delta_ee_gl != NULL) { rc = qmckl_free(context, ctx->quad_point.delta_ee_gl); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_jastrow_champ_quad_ee_gl", "Unable to free ctx->quad_point.delta_ee_gl"); } ctx->quad_point.delta_ee_gl = NULL; } } /* Allocate array */ if (ctx->quad_point.delta_ee_gl == NULL) { double* delta_ee_gl = (double*) qmckl_malloc(context, mem_info); ctx->quad_point.delta_ee_gl_maxsize = mem_info.size; if (delta_ee_gl == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_jastrow_champ_quad_ee_gl", NULL); } ctx->quad_point.delta_ee_gl = delta_ee_gl; } rc = qmckl_compute_jastrow_champ_quad_ee_gl(context, ctx->quad_point.num, ctx->quad_point.indices, 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->quad_point.ee_rescaled_quad, ctx->quad_point.ee_rescaled_quad_gl, ctx->jastrow_champ.spin_independent, ctx->quad_point.delta_ee_gl); if (rc != QMCKL_SUCCESS) { return rc; } ctx->quad_point.delta_ee_gl_date = ctx->quad_point.date; } return QMCKL_SUCCESS; } #+end_src *** Compute :PROPERTIES: :Name: qmckl_compute_jastrow_champ_quad_ee_gl :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_quad_ee_gl_args | Variable | Type | In/Out | Description | |---------------------------+-------------------------------------------+--------+-----------------------------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~num~ | ~int64_t~ | in | Index of quad electron | | ~indices~ | ~int64_t[num]~ | in | Index of quad electron | | ~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 rescaled distances | | ~ee_distance_rescaled_gl~ | ~double[walk_num][4][elec_num][elec_num]~ | in | Electron-electron rescaled distances derivatives | | ~ee_rescaled_quad~ | ~double[num][elec_num]~ | in | Electron-electron rescaled quad distances | | ~ee_rescaled_quad_gl~ | ~double[num][3][elec_num]~ | in | Electron-electron rescaled quad distances derivatives | | ~spin_independent~ | ~int32_t~ | in | If 1, same parameters for parallel and antiparallel spins | | ~delta_ee_gl~ | ~double[num][elec_num][3]~ | out | quad electron-electron jastrow gradients and Laplacian | |---------------------------+-------------------------------------------+--------+-----------------------------------------------------------| #+begin_src f90 :comments org :tangle (eval f) :noweb yes function qmckl_compute_jastrow_champ_quad_ee_gl_doc( & context, num, indices, walk_num, elec_num, up_num, bord_num, & b_vector, ee_distance_rescaled, ee_distance_rescaled_gl, & ee_rescaled_quad, ee_rescaled_quad_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 integer (c_int64_t) , intent(in) :: indices(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 :: 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_quad(elec_num,num) real (c_double ) , intent(in) :: ee_rescaled_quad_gl(3,elec_num,num) integer (c_int32_t) , intent(in) , value :: spin_independent real (c_double ) , intent(out) :: delta_ee_gl(3,elec_num,num) integer(qmckl_exit_code) :: info integer*8 :: i, j, k, nw, ii 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(3), dx_old(3) 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, num delta_ee_gl(:,:,nw) = 0.0d0 do i = 1, elec_num if (i == (indices(nw)+1)) cycle x = ee_rescaled_quad(i,nw) x_old = ee_distance_rescaled(i,indices(nw)+1,1) 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_quad_gl(1, i, nw) dx(2) = ee_rescaled_quad_gl(2, i, nw) dx(3) = ee_rescaled_quad_gl(3, i, nw) dx_old(1) = ee_distance_rescaled_gl(1, i, indices(nw)+1, 1) dx_old(2) = ee_distance_rescaled_gl(2, i, indices(nw)+1, 1) dx_old(3) = ee_distance_rescaled_gl(3, i, indices(nw)+1, 1) 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. (indices(nw)+1) <= up_num ) .or. (i > up_num .and. (indices(nw)+1) > 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(1,indices(nw)+1,nw) = delta_ee_gl(1,indices(nw)+1,nw) - f * dx(1) + f_old * dx_old(1) delta_ee_gl(2,indices(nw)+1,nw) = delta_ee_gl(2,indices(nw)+1,nw) - f * dx(2) + f_old * dx_old(2) delta_ee_gl(3,indices(nw)+1,nw) = delta_ee_gl(3,indices(nw)+1,nw) - f * dx(3) + f_old * dx_old(3) 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(1,indices(nw)+1,nw) = delta_ee_gl(1,indices(nw)+1,nw) - f * x1 * dx(1) + f_old * x1_old * dx_old(1) delta_ee_gl(2,indices(nw)+1,nw) = delta_ee_gl(2,indices(nw)+1,nw) - f * x1 * dx(2) + f_old * x1_old * dx_old(2) delta_ee_gl(3,indices(nw)+1,nw) = delta_ee_gl(3,indices(nw)+1,nw) - f * x1 * dx(3) + f_old * x1_old * dx_old(3) 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_quad_ee_gl_doc #+end_src #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_jastrow_champ_quad_ee_gl (const qmckl_context context, const int64_t num, const int64_t* indices, 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_quad, const double* ee_rescaled_quad_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_quad_ee_gl_doc (const qmckl_context context, const int64_t num, const int64_t* indices, 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_quad, const double* ee_rescaled_quad_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_quad_ee_gl (const qmckl_context context, const int64_t num, const int64_t* indices, 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_quad, const double* ee_rescaled_quad_gl, const int32_t spin_independent, double* const delta_ee_gl ) { #ifdef HAVE_HPC return qmckl_compute_jastrow_champ_quad_ee_gl_doc #else return qmckl_compute_jastrow_champ_quad_ee_gl_doc #endif (context, num, indices, walk_num, elec_num, up_num, bord_num, b_vector, ee_distance_rescaled, ee_distance_rescaled_gl, ee_rescaled_quad, ee_rescaled_quad_gl, spin_independent, delta_ee_gl ); } #+end_src *** Test :noexport: #+begin_src c :tangle (eval c_test) printf("Delta ee gl\n"); /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); double ee_gl_old[walk_num][4][elec_num]; double delta_ee_gl[num][elec_num][3]; double ee_gl_new[walk_num][4][elec_num]; for (int elec = 0; elec < num; elec++){ rc = qmckl_set_electron_coord(context, 'N', walk_num, elec_coord, walk_num*elec_num*3); assert(rc == QMCKL_SUCCESS); 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_quad_points(context, 'N', num, indices, &new_coords[0][0], 3*num); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_jastrow_champ_quad_ee_gl(context, &delta_ee_gl[0][0][0], num*elec_num*4); assert (rc == QMCKL_SUCCESS); coords[0][indices[elec]][0] = new_coords[elec][0]; coords[0][indices[elec]][1] = new_coords[elec][1]; coords[0][indices[elec]][2] = new_coords[elec][2]; rc = qmckl_set_electron_coord(context, 'N', walk_num, &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); 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 i = 0; i < elec_num; i++) { for (int m = 0; m < 3; m++) { //printf("%f\n",(ee_gl_new[0][m][i] - ee_gl_old[0][m][i])); //printf("%f\n",delta_ee_gl[elec][i][m]); assert(fabs((ee_gl_new[0][m][i] - ee_gl_old[0][m][i]) - delta_ee_gl[elec][i][m]) < 1.e-12); } } } printf("OK\n"); #+end_src * Electron-nucleus Jastrow ** Electron-nucleus rescaled distance *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_en_rescaled_quad(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_quad(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_quad(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->quad_point.num; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "todo", "Array too small. Expected ctx->nucleus.num * ctx->quad_point.num "); } memcpy(distance_rescaled, ctx->quad_point.en_rescaled_quad, 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_quad(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_provide_en_rescaled_quad(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_quad_en_distance(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_en_distance(context); if(rc != QMCKL_SUCCESS) return rc; /* Compute if necessary */ if (ctx->quad_point.date > ctx->quad_point.en_rescaled_quad_date) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = ctx->nucleus.num * ctx->quad_point.num * sizeof(double); if (mem_info.size > ctx->quad_point.en_rescaled_quad_maxsize) { if (ctx->quad_point.en_rescaled_quad!= NULL) { rc = qmckl_free(context, ctx->quad_point.en_rescaled_quad); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_en_rescaled_quad", "Unable to free ctx->quad_point.en_rescaled_quad"); } ctx->quad_point.en_rescaled_quad = NULL; } } /* Allocate array */ if (ctx->quad_point.en_rescaled_quad == NULL) { double* en_rescaled_quad = (double*) qmckl_malloc(context, mem_info); ctx->quad_point.en_rescaled_quad_maxsize = mem_info.size; if (en_rescaled_quad == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_en_rescaled_quad", NULL); } ctx->quad_point.en_rescaled_quad = en_rescaled_quad; } rc = qmckl_compute_en_rescaled_quad(context, ctx->quad_point.num, 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->quad_point.quad_en_distance, ctx->quad_point.en_rescaled_quad); if (rc != QMCKL_SUCCESS) { return rc; } ctx->quad_point.en_rescaled_quad_date = ctx->quad_point.date; } return QMCKL_SUCCESS; } #+end_src *** Compute :PROPERTIES: :Name: qmckl_compute_en_rescaled_quad :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_en_rescaled_quad_args | Variable | Type | In/Out | Description | |----------------------+------------------------------+--------+-------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | | ~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 | | ~quad_en_distance~ | ~double[num][nucl_num]~ | in | quad electron-nucleus distances | | ~en_rescaled_quad~ | ~double[num][nucl_num]~ | out | Electron-nucleus rescaled distances | |----------------------+------------------------------+--------+-------------------------------------| #+begin_src f90 :comments org :tangle (eval f) :noweb yes function qmckl_compute_en_rescaled_quad_doc(context, num,& nucl_num, type_nucl_num, type_nucl_vector, rescale_factor_en, & walk_num, quad_en_distance, en_rescaled_quad) & bind(C) result(info) use qmckl implicit none integer (qmckl_context), intent(in), value :: context integer (c_int64_t) , intent(in) , value :: 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) 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) :: quad_en_distance(nucl_num,num) real (c_double ) , intent(out) :: en_rescaled_quad(nucl_num,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 k=1,num do i=1, nucl_num en_rescaled_quad(i,k) = (1.0d0 - dexp(-rescale_factor_en(type_nucl_vector(i)+1) * & quad_en_distance(i,k))) / rescale_factor_en(type_nucl_vector(i)+1) end do end do end function qmckl_compute_en_rescaled_quad_doc #+end_src #+begin_src c :tangle (eval h_private_func) :comments org :exports none qmckl_exit_code qmckl_compute_en_rescaled_quad_doc ( const qmckl_context context, const int64_t num, 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* quad_en_distance, double* const en_rescaled_quad ); qmckl_exit_code qmckl_compute_en_rescaled_quad ( const qmckl_context context, const int64_t num, 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* quad_en_distance, double* const en_rescaled_quad ); #+end_src #+begin_src c :tangle (eval c) :comments org :exports none qmckl_exit_code qmckl_compute_en_rescaled_quad( const qmckl_context context, const int64_t num, 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* quad_en_distance, double* const en_rescaled_quad ) { #ifdef HAVE_HPC return qmckl_compute_en_rescaled_quad_doc #else return qmckl_compute_en_rescaled_quad_doc #endif (context, num, nucl_num, type_nucl_num, type_nucl_vector, rescale_factor_en, walk_num, quad_en_distance, en_rescaled_quad ); } #+end_src *** Test :noexport: #+begin_src c :tangle (eval c_test) printf("En rescaled quad\n"); /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); double en_rescaled[walk_num][nucl_num][elec_num]; double quad_en_rescaled[num][nucl_num]; for (int elec = 0; elec < num; elec++){ rc = qmckl_set_electron_coord(context, 'N', walk_num, elec_coord, walk_num*elec_num*3); assert(rc == QMCKL_SUCCESS); rc = qmckl_get_jastrow_champ_en_distance_rescaled(context, &en_rescaled[0][0][0], walk_num*nucl_num*elec_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_quad_points(context, 'N', num, indices, &new_coords[0][0], 3*num); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_en_rescaled_quad(context, &quad_en_rescaled[0][0], num*nucl_num); assert (rc == QMCKL_SUCCESS); coords[0][indices[elec]][0] = new_coords[elec][0]; coords[0][indices[elec]][1] = new_coords[elec][1]; coords[0][indices[elec]][2] = new_coords[elec][2]; rc = qmckl_set_electron_coord(context, 'N', walk_num, &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_jastrow_champ_en_distance_rescaled(context, &en_rescaled[0][0][0], walk_num*nucl_num*elec_num); assert (rc == QMCKL_SUCCESS); for (int a = 0; a < nucl_num; a++){ assert(fabs(en_rescaled[0][a][indices[elec]]-quad_en_rescaled[elec][a]) < 1.e-12); } } printf("OK\n"); #+end_src ** quad electron-nucleus Jastrow value *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_jastrow_champ_quad_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_quad_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_quad_en", NULL); } qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); qmckl_exit_code rc; rc = qmckl_provide_jastrow_champ_quad_en(context); if (rc != QMCKL_SUCCESS) return rc; int64_t sze=ctx->quad_point.num; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_jastrow_champ_quad_en", "Array too small. Expected quad_point.num"); } memcpy(delta_en, ctx->quad_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_quad_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_quad_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_quad_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_quad_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_quad_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_quad_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_quad(context); if(rc != QMCKL_SUCCESS) return rc; /* Compute if necessary */ if (ctx->quad_point.date > ctx->quad_point.delta_en_date) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = ctx->quad_point.num * sizeof(double); if (mem_info.size > ctx->quad_point.delta_en_maxsize) { if (ctx->quad_point.delta_en != NULL) { rc = qmckl_free(context, ctx->quad_point.delta_en); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_jastrow_champ_quad_en", "Unable to free ctx->quad_point.delta_en"); } ctx->quad_point.delta_en = NULL; } } /* Allocate array */ if (ctx->quad_point.delta_en == NULL) { double* delta_en = (double*) qmckl_malloc(context, mem_info); ctx->quad_point.delta_en_maxsize = mem_info.size; if (delta_en == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_jastrow_champ_quad_en", NULL); } ctx->quad_point.delta_en = delta_en; } rc = qmckl_compute_jastrow_champ_quad_en(context, ctx->quad_point.num, ctx->quad_point.indices, 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->quad_point.en_rescaled_quad, ctx->quad_point.delta_en); if (rc != QMCKL_SUCCESS) { return rc; } ctx->quad_point.delta_en_date = ctx->quad_point.date; } return QMCKL_SUCCESS; } #+end_src *** Compute :PROPERTIES: :Name: qmckl_compute_jastrow_champ_quad_en_doc :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_quad_en_args |------------------------+----------------------------------------+--------+--------------------------------------------| | Variable | Type | In/Out | Description | |------------------------+----------------------------------------+--------+--------------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~num~ | ~int64_t~ | in | Index of quad point | | ~indices~ | ~int64_t[num]~ | in | Indices of electrons | | ~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 rescaled distances | | ~en_rescaled_quad ~ | ~double[num][nucl_num]~ | in | Electron-nucleus rescaled quad distances | | ~delta_en~ | ~double[num]~ | out | quad electron-nucleus jastrow | |------------------------+----------------------------------------+--------+--------------------------------------------| #+begin_src f90 :comments org :tangle (eval f) :noweb yes function qmckl_compute_jastrow_champ_quad_en_doc( & context, num, indices, walk_num, elec_num, nucl_num, type_nucl_num, & type_nucl_vector, aord_num, a_vector, & en_distance_rescaled, en_rescaled_quad, 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 integer (c_int64_t) , intent(in) :: indices(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 :: 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_quad(nucl_num,num) real (c_double ) , intent(out) :: delta_en(num) integer(qmckl_exit_code) :: info integer*8 :: i, a, p, nw double precision :: x, power_ser, y 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, num delta_en(nw) = 0.0d0 do a = 1, nucl_num x = en_distance_rescaled(indices(nw)+1, a, 1) y = en_rescaled_quad(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(indices(nw)+1, a, 1) y = y * en_rescaled_quad(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_quad_en_doc #+end_src #+CALL: generate_c_header(table=qmckl_quad_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_quad_en ( const qmckl_context context, const int64_t num, const int64_t* indices, 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_quad, double* const delta_en ); qmckl_exit_code qmckl_compute_jastrow_champ_quad_en_doc ( const qmckl_context context, const int64_t num, const int64_t* indices, 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_quad, double* const delta_en ); #+end_src #+begin_src c :tangle (eval c) :comments org qmckl_exit_code qmckl_compute_jastrow_champ_quad_en ( const qmckl_context context, const int64_t num, const int64_t* indices, 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_quad, double* const delta_en ) { #ifdef HAVE_HPC return qmckl_compute_jastrow_champ_quad_en_doc #else return qmckl_compute_jastrow_champ_quad_en_doc #endif (context, num, indices, walk_num, elec_num, nucl_num, type_nucl_num, type_nucl_vector, aord_num, a_vector, en_distance_rescaled, en_rescaled_quad, delta_en ); } #+end_src *** Test :noexport: #+begin_src c :tangle (eval c_test) printf("Delta en\n"); /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); double jastrow_en_old[walk_num]; double delta_en[num]; double jastrow_en_new[walk_num]; for (int elec = 0; elec < num; elec++){ rc = qmckl_set_electron_coord(context, 'N', walk_num, elec_coord, walk_num*elec_num*3); assert(rc == QMCKL_SUCCESS); 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_quad_points(context, 'N', num, indices, &new_coords[0][0], 3*num); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_jastrow_champ_quad_en(context, &delta_en[0], num); assert (rc == QMCKL_SUCCESS); coords[0][indices[elec]][0] = new_coords[elec][0]; coords[0][indices[elec]][1] = new_coords[elec][1]; coords[0][indices[elec]][2] = new_coords[elec][2]; rc = qmckl_set_electron_coord(context, 'N', walk_num, &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_jastrow_champ_factor_en(context, &jastrow_en_new[0], walk_num); assert (rc == QMCKL_SUCCESS); //printf("electron %d walk %d \n", elec, 0); //printf("jastrow_en_new %f\n", jastrow_en_new[0]); //printf("jastrow_en_old %f\n", jastrow_en_old[0]); //printf("delta_en %f\n", delta_en[elec]); //printf("diff %f\n", jastrow_en_new[0] - jastrow_en_old[0] - delta_en[elec]); assert(fabs((jastrow_en_new[0] - jastrow_en_old[0]) - delta_en[elec]) < 1.e-12); } printf("OK\n"); #+end_src ** Electron-nucleus rescaled distances derivatives *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_en_rescaled_quad_gl(qmckl_context context, double* distance_rescaled_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_en_rescaled_quad_gl(qmckl_context context, double* distance_rescaled_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_en_rescaled_quad_gl(context); if (rc != QMCKL_SUCCESS) return rc; qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); size_t sze = 3 * ctx->nucleus.num * ctx->quad_point.num; memcpy(distance_rescaled_gl, ctx->quad_point.en_rescaled_quad_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_en_rescaled_quad_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_quad_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_quad_en_distance(context); if(rc != QMCKL_SUCCESS) return rc; /* Compute if necessary */ if (ctx->quad_point.date > ctx->quad_point.en_rescaled_quad_gl_date) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = 3 * ctx->nucleus.num * ctx->quad_point.num * sizeof(double); if (mem_info.size > ctx->quad_point.en_rescaled_quad_gl_maxsize) { if (ctx->quad_point.en_rescaled_quad_gl != NULL) { rc = qmckl_free(context, ctx->quad_point.en_rescaled_quad_gl); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_en_rescaled_quad_gl", "Unable to free ctx->quad_point.en_rescaled_quad_gl"); } ctx->quad_point.en_rescaled_quad_gl = NULL; } } /* Allocate array */ if (ctx->quad_point.en_rescaled_quad_gl == NULL) { double* en_rescaled_quad_gl = (double*) qmckl_malloc(context, mem_info); ctx->quad_point.en_rescaled_quad_gl_maxsize = mem_info.size; if (en_rescaled_quad_gl == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_en_rescaled_quad_gl", NULL); } ctx->quad_point.en_rescaled_quad_gl = en_rescaled_quad_gl; } qmckl_exit_code rc = qmckl_compute_en_rescaled_quad_gl(context, ctx->quad_point.num, 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->quad_point.quad_en_distance, ctx->quad_point.coord.data, ctx->nucleus.coord.data, ctx->quad_point.en_rescaled_quad_gl); if (rc != QMCKL_SUCCESS) { return rc; } ctx->quad_point.en_rescaled_quad_gl_date = ctx->quad_point.date; } return QMCKL_SUCCESS; } #+end_src *** Compute :PROPERTIES: :Name: qmckl_compute_en_rescaled_quad_gl :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_en_rescaled_quad_gl_args | Variable | Type | In/Out | Description | |-------------------------+---------------------------------+--------+-------------------------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~num~ | ~int64_t~ | in | Number of quad points | | ~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 | | ~quad_en_distance~ | ~double[num][nucl_num]~ | in | quad electorn distances | | ~coord~ | ~double[num][3]~ | in | quad electron coordinates | | ~nucl_coord~ | ~double[3][nucl_num]~ | in | Nucleus coordinates | | ~en_rescaled_quad_gl~ | ~double[num][nucl_num][3]~ | out | Electron-nucleus rescaled quad distance derivatives | |-------------------------+---------------------------------+--------+-------------------------------------------------------| #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_en_rescaled_quad_gl_doc_f(context, num, nucl_num, & type_nucl_num, type_nucl_vector, rescale_factor_en, walk_num, & quad_en_distance, coord, nucl_coord, en_rescaled_quad_gl) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: num 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) :: quad_en_distance(nucl_num, num) double precision , intent(in) :: coord(3,num) double precision , intent(in) :: nucl_coord(nucl_num,3) double precision , intent(out) :: en_rescaled_quad_gl(3,nucl_num,num) integer*8 :: nw, a, ii double precision :: ria_inv, elnuc_dist_gl(3, 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_quad_gl = 0.0d0 do nw = 1, num ! prepare the actual een table do a = 1, nucl_num ria_inv = 1.0d0 / quad_en_distance(a, nw) do ii = 1, 3 elnuc_dist_gl(ii, a) = (coord(ii,nw) - nucl_coord(a, ii)) * ria_inv end do end do do a = 1, nucl_num kappa_l = -1 * rescale_factor_en(type_nucl_vector(a)+1) en_rescaled_quad_gl(1, a, nw) = elnuc_dist_gl(1, a) en_rescaled_quad_gl(2, a, nw) = elnuc_dist_gl(2, a) en_rescaled_quad_gl(3, a, nw) = elnuc_dist_gl(3, a) en_rescaled_quad_gl(1, a, nw) = en_rescaled_quad_gl(1, a, nw) * dexp(kappa_l * quad_en_distance(a,nw)) en_rescaled_quad_gl(2, a, nw) = en_rescaled_quad_gl(2, a, nw) * dexp(kappa_l * quad_en_distance(a,nw)) en_rescaled_quad_gl(3, a, nw) = en_rescaled_quad_gl(3, a, nw) * dexp(kappa_l * quad_en_distance(a,nw)) end do end do end function qmckl_compute_en_rescaled_quad_gl_doc_f #+end_src #+begin_src c :tangle (eval h_private_func) :comments org :exports none qmckl_exit_code qmckl_compute_en_rescaled_quad_gl_doc ( const qmckl_context context, const int64_t num, 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* quad_en_distance, const double* coord, const double* nucl_coord, double* const en_rescaled_quad_gl ); qmckl_exit_code qmckl_compute_en_rescaled_quad_gl ( const qmckl_context context, const int64_t num, 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* quad_en_distance, const double* coord, const double* nucl_coord, double* const en_rescaled_quad_gl ); #+end_src #+begin_src c :tangle (eval c) :comments org :exports none qmckl_exit_code qmckl_compute_en_rescaled_quad_gl ( const qmckl_context context, const int64_t num, 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* quad_en_distance, const double* coord, const double* nucl_coord, double* const en_rescaled_quad_gl ) { #ifdef HAVE_HPC return qmckl_compute_en_rescaled_quad_gl_doc #else return qmckl_compute_en_rescaled_quad_gl_doc #endif (context, num, nucl_num, type_nucl_num, type_nucl_vector, rescale_factor_en, walk_num, quad_en_distance, coord, nucl_coord, en_rescaled_quad_gl ); } #+end_src #+CALL: generate_c_interface(table=qmckl_en_rescaled_quad_gl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_en_rescaled_quad_gl_doc") #+RESULTS: #+begin_src f90 :tangle (eval f) :comments org :exports none integer(c_int32_t) function qmckl_compute_en_rescaled_quad_gl_doc & (context, & num, & nucl_num, & type_nucl_num, & type_nucl_vector, & rescale_factor_en, & walk_num, & quad_en_distance, & coord, & nucl_coord, & en_rescaled_quad_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 :: 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) :: quad_en_distance(nucl_num,num) real (c_double ) , intent(in) :: coord(3,num) real (c_double ) , intent(in) :: nucl_coord(nucl_num,3) real (c_double ) , intent(out) :: en_rescaled_quad_gl(3,nucl_num,num) integer(c_int32_t), external :: qmckl_compute_en_rescaled_quad_gl_doc_f info = qmckl_compute_en_rescaled_quad_gl_doc_f & (context, & num, & nucl_num, & type_nucl_num, & type_nucl_vector, & rescale_factor_en, & walk_num, & quad_en_distance, & coord, & nucl_coord, & en_rescaled_quad_gl) end function qmckl_compute_en_rescaled_quad_gl_doc #+end_src *** Test :noexport: #+begin_src c :tangle (eval c_test) printf("En rescaled quad gl\n"); /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); double en_rescaled_gl[walk_num][nucl_num][elec_num][4]; double quad_en_rescaled_gl[num][nucl_num][3]; for (int elec = 0; elec < num; elec++){ rc = qmckl_set_electron_coord(context, 'N', walk_num, elec_coord, walk_num*elec_num*3); assert(rc == QMCKL_SUCCESS); rc = qmckl_get_jastrow_champ_en_distance_rescaled_gl(context, &en_rescaled_gl[0][0][0][0], walk_num*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_quad_points(context, 'N', num, indices, &new_coords[0][0], 3*num); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_en_rescaled_quad_gl(context, &quad_en_rescaled_gl[0][0][0], num*nucl_num*3); assert (rc == QMCKL_SUCCESS); coords[0][indices[elec]][0] = new_coords[elec][0]; coords[0][indices[elec]][1] = new_coords[elec][1]; coords[0][indices[elec]][2] = new_coords[elec][2]; rc = qmckl_set_electron_coord(context, 'N', walk_num, &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_jastrow_champ_en_distance_rescaled_gl(context, &en_rescaled_gl[0][0][0][0], walk_num*nucl_num*elec_num*4); assert (rc == QMCKL_SUCCESS); for (int a = 0; a < nucl_num; a++) { for (int m = 0; m < 3; m++) { assert(fabs(en_rescaled_gl[0][a][indices[elec]][m] - quad_en_rescaled_gl[elec][a][m]) < 1.e-12); } } } printf("OK\n"); #+end_src ** Electron-nucleus Jastrow gradients and Laplacian *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_jastrow_champ_quad_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_quad_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_quad_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->quad_point.num * 3 * ctx->electron.num; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_jastrow_champ_quad_en_gl", "Array too small. Expected 3*quad_point.num*elec_num"); } memcpy(delta_en_gl, ctx->quad_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_quad_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_quad_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_quad_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_quad_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_quad_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_quad_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_quad(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if en rescaled distance derivatives is provided */ rc = qmckl_provide_en_rescaled_quad_gl(context); if(rc != QMCKL_SUCCESS) return rc; /* Compute if necessary */ if (ctx->quad_point.date > ctx->quad_point.delta_en_gl_date) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = ctx->quad_point.num * 3 * ctx->electron.num * sizeof(double); if (mem_info.size > ctx->quad_point.delta_en_gl_maxsize) { if (ctx->quad_point.delta_en_gl != NULL) { rc = qmckl_free(context, ctx->quad_point.delta_en_gl); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_jastrow_champ_quad_en_gl", "Unable to free ctx->quad_point.delta_en_gl"); } ctx->quad_point.delta_en_gl = NULL; } } /* Allocate array */ if (ctx->quad_point.delta_en_gl == NULL) { double* delta_en_gl = (double*) qmckl_malloc(context, mem_info); ctx->quad_point.delta_en_gl_maxsize = mem_info.size; if (delta_en_gl == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_jastrow_champ_quad_en_gl", NULL); } ctx->quad_point.delta_en_gl = delta_en_gl; } rc = qmckl_compute_jastrow_champ_quad_en_gl(context, ctx->quad_point.num, ctx->quad_point.indices, 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->quad_point.en_rescaled_quad, ctx->quad_point.en_rescaled_quad_gl, ctx->quad_point.delta_en_gl); if (rc != QMCKL_SUCCESS) { return rc; } ctx->quad_point.delta_en_gl_date = ctx->quad_point.date; } return QMCKL_SUCCESS; } #+end_src *** Compute :PROPERTIES: :Name: qmckl_compute_jastrow_champ_quad_en_gl :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_quad_en_gl_args | Variable | Type | In/Out | Description | |---------------------------+-------------------------------------------+--------+---------------------------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~num~ | ~int64_t~ | in | Index of quad electron | | ~indices~ | ~int64_t[num]~ | in | Indices of quad electrons | | ~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 rescaled distances | | ~en_distance_rescaled_gl~ | ~double[walk_num][nucl_num][elec_num][4]~ | in | Electron-nucleus rescaled distance derivatives | | ~en_rescaled_quad~ | ~double[num][nucl_num]~ | in | Electron-nucleus rescaled quad distances | | ~en_rescaled_quad_gl~ | ~double[num][nucl_num][3]~ | in | Electron-nucleus rescaled quad distance derivatives | | ~delta_en_gl~ | ~double[num][elec_num][3]~ | out | quad electron-nucleus Jastrow gradients and Laplacian | |---------------------------+-------------------------------------------+--------+---------------------------------------------------------| #+begin_src f90 :comments org :tangle (eval f) :noweb yes function qmckl_compute_jastrow_champ_quad_en_gl_doc( & context, num, indices, 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_quad, en_rescaled_quad_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 integer (c_int64_t) , intent(in) :: indices(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 :: 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_quad(nucl_num,num) real (c_double ) , intent(in) :: en_rescaled_quad_gl(3, nucl_num,num) real (c_double ) , intent(out) :: delta_en_gl(3,elec_num,num) integer(qmckl_exit_code) :: info integer*8 :: i, a, k, nw, ii 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(3), dx_old(3) 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, num delta_en_gl(:,:,nw) = 0.0d0 do a = 1, nucl_num x_old = en_distance_rescaled(indices(nw)+1,a,1) x = en_rescaled_quad(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_quad_gl(1,a,nw) dx(2) = en_rescaled_quad_gl(2,a,nw) dx(3) = en_rescaled_quad_gl(3,a,nw) dx_old(1) = en_distance_rescaled_gl(1,indices(nw)+1,a,1) dx_old(2) = en_distance_rescaled_gl(2,indices(nw)+1,a,1) dx_old(3) = en_distance_rescaled_gl(3,indices(nw)+1,a,1) 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,indices(nw)+1,nw) = delta_en_gl(1,indices(nw)+1,nw) + f * dx(1) - f_old * dx_old(1) delta_en_gl(2,indices(nw)+1,nw) = delta_en_gl(2,indices(nw)+1,nw) + f * dx(2) - f_old * dx_old(2) delta_en_gl(3,indices(nw)+1,nw) = delta_en_gl(3,indices(nw)+1,nw) + f * dx(3) - f_old * dx_old(3) 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,indices(nw)+1,nw) = delta_en_gl(1,indices(nw)+1,nw) + f * x1 * dx(1) - f_old * x1_old * dx_old(1) delta_en_gl(2,indices(nw)+1,nw) = delta_en_gl(2,indices(nw)+1,nw) + f * x1 * dx(2) - f_old * x1_old * dx_old(2) delta_en_gl(3,indices(nw)+1,nw) = delta_en_gl(3,indices(nw)+1,nw) + f * x1 * dx(3) - f_old * x1_old * dx_old(3) 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_quad_en_gl_doc #+end_src #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_jastrow_champ_quad_en_gl_doc ( const qmckl_context context, const int64_t num, const int64_t* indices, 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_quad, const double* en_rescaled_quad_gl, double* const delta_en_gl ); qmckl_exit_code qmckl_compute_jastrow_champ_quad_en_gl ( const qmckl_context context, const int64_t num, const int64_t* indices, 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_quad, const double* en_rescaled_quad_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_quad_en_gl (const qmckl_context context, const int64_t num, const int64_t* indices, 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_quad, const double* en_rescaled_quad_gl, double* const delta_en_gl ) { #ifdef HAVE_HPC return qmckl_compute_jastrow_champ_quad_en_gl_doc #else return qmckl_compute_jastrow_champ_quad_en_gl_doc #endif (context, num, indices, 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_quad, en_rescaled_quad_gl, delta_en_gl ); } #+end_src *** Test :noexport: #+begin_src c :tangle (eval c_test) printf("Delta en gl\n"); /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); double en_gl_old[walk_num][4][elec_num]; double delta_en_gl[num][elec_num][3]; double en_gl_new[walk_num][4][elec_num]; for (int elec = 0; elec < num; elec++){ rc = qmckl_set_electron_coord(context, 'N', walk_num, elec_coord, walk_num*elec_num*3); assert(rc == QMCKL_SUCCESS); 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_quad_points(context, 'N', num, indices, &new_coords[0][0], 3*num); assert (rc == QMCKL_SUCCESS); rc = qmckl_get_jastrow_champ_quad_en_gl(context, &delta_en_gl[0][0][0], num*elec_num*4); assert (rc == QMCKL_SUCCESS); coords[0][indices[elec]][0] = new_coords[elec][0]; coords[0][indices[elec]][1] = new_coords[elec][1]; coords[0][indices[elec]][2] = new_coords[elec][2]; rc = qmckl_set_electron_coord(context, 'N', walk_num, &coords[0][0][0], walk_num*elec_num*3); assert (rc == QMCKL_SUCCESS); 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 i = 0; i < elec_num; i++) { for (int m = 0; m < 3; m++) { assert(fabs((en_gl_new[0][m][i] - en_gl_old[0][m][i]) - delta_en_gl[elec][i][m]) < 1.e-12); } } } printf("OK\n"); #+end_src * Force of quad en Jastrow ** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_forces_jastrow_quad_en(qmckl_context context, double* const forces_jastrow_quad_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_forces_jastrow_quad_en(qmckl_context context, double* const forces_jastrow_quad_en, const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_exit_code rc; rc = qmckl_provide_forces_jastrow_quad_en(context); if (rc != QMCKL_SUCCESS) return rc; qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); const int64_t sze = ctx->quad_point.num * 3 * ctx->nucleus.num; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_forces_jastrow_quad_en", "input array too small"); } memcpy(forces_jastrow_quad_en, ctx->quad_point.forces_jastrow_quad_en, sze * sizeof(double)); return QMCKL_SUCCESS; } #+end_src #+begin_src f90 :tangle (eval fh_func) :comments org :exports none interface integer(qmckl_exit_code) function qmckl_get_forces_jastrow_quad_en (context, & forces_jastrow_quad_en, size_max) bind(C) use, intrinsic :: iso_c_binding import implicit none integer (c_int64_t) , intent(in) , value :: context integer (c_int64_t) , intent(in) , value :: size_max real(c_double), intent(out) :: forces_jastrow_quad_en(size_max) end function qmckl_get_forces_jastrow_quad_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_forces_jastrow_quad_en(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :export none qmckl_exit_code qmckl_provide_forces_jastrow_quad_en(qmckl_context context) { qmckl_exit_code rc = QMCKL_SUCCESS; if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return qmckl_failwith( context, QMCKL_INVALID_CONTEXT, "qmckl_provide_forces_jastrow_quad_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_forces_jastrow_quad_en", 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_quad(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if en rescaled distance derivatives is provided */ rc = qmckl_provide_en_rescaled_quad_gl(context); if(rc != QMCKL_SUCCESS) return rc; /* Compute if necessary */ if (ctx->quad_point.date > ctx->quad_point.forces_jastrow_quad_en_date) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = ctx->quad_point.num * 3 * ctx->nucleus.num * sizeof(double); if (mem_info.size > ctx->quad_point.forces_jastrow_quad_en_maxsize) { if (ctx->quad_point.forces_jastrow_quad_en != NULL) { rc = qmckl_free(context, ctx->quad_point.forces_jastrow_quad_en); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_forces_jastrow_quad_en", "Unable to free ctx->quad_point.forces_jastrow_quad_en"); } ctx->quad_point.forces_jastrow_quad_en = NULL; } } /* Allocate array */ if (ctx->quad_point.forces_jastrow_quad_en == NULL) { double* forces_jastrow_quad_en = (double*) qmckl_malloc(context, mem_info); ctx->quad_point.forces_jastrow_quad_en_maxsize = mem_info.size; if (forces_jastrow_quad_en == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_forces_jastrow_quad_en", NULL); } ctx->quad_point.forces_jastrow_quad_en = forces_jastrow_quad_en; } rc = qmckl_compute_forces_jastrow_quad_en(context, ctx->quad_point.num, ctx->quad_point.indices, 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->quad_point.en_rescaled_quad, ctx->quad_point.en_rescaled_quad_gl, ctx->quad_point.forces_jastrow_quad_en); if (rc != QMCKL_SUCCESS) { return rc; } ctx->quad_point.forces_jastrow_quad_en_date = ctx->quad_point.date; } return QMCKL_SUCCESS; } #+end_src ** Compute :PROPERTIES: :Name: qmckl_compute_forces_jastrow_quad_en :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_forces_quad_en_args | Variable | Type | In/Out | Description | |---------------------------+-------------------------------------------+--------+---------------------------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~num~ | ~int64_t~ | in | Index of quad electron | | ~indices~ | ~int64_t[num]~ | in | Indices of quad electrons | | ~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 rescaled distances | | ~en_distance_rescaled_gl~ | ~double[walk_num][nucl_num][elec_num][4]~ | in | Electron-nucleus rescaled distance derivatives | | ~en_rescaled_quad~ | ~double[num][nucl_num]~ | in | Electron-nucleus rescaled quad distances | | ~en_rescaled_quad_gl~ | ~double[num][nucl_num][3]~ | in | Electron-nucleus rescaled quad distance derivatives | | ~forces_jastrow_quad_en~ | ~double[num][nucl_num][3]~ | out | quad electron-nucleus forces | |---------------------------+-------------------------------------------+--------+---------------------------------------------------------| #+begin_src f90 :comments org :tangle (eval f) :noweb yes function qmckl_compute_forces_jastrow_quad_en_doc( & context, num, indices, 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_quad, en_rescaled_quad_gl, forces_jastrow_quad_en) & bind(C) result(info) use qmckl implicit none integer (qmckl_context), intent(in), value :: context integer (c_int64_t) , intent(in) , value :: num integer (c_int64_t) , intent(in) :: indices(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 :: 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_quad(nucl_num,num) real (c_double ) , intent(in) :: en_rescaled_quad_gl(3, nucl_num,num) real (c_double ) , intent(out) :: forces_jastrow_quad_en(3,nucl_num,num) integer(qmckl_exit_code) :: info integer*8 :: i, a, k, nw, ii 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 :: dx(3), dx_old(3) 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, num forces_jastrow_quad_en(:,:,nw) = 0.0d0 do a = 1, nucl_num x_old = en_distance_rescaled(indices(nw)+1,a,1) x = en_rescaled_quad(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_quad_gl(1,a,nw) dx(2) = -en_rescaled_quad_gl(2,a,nw) dx(3) = -en_rescaled_quad_gl(3,a,nw) dx_old(1) = -en_distance_rescaled_gl(1,indices(nw)+1,a,1) dx_old(2) = -en_distance_rescaled_gl(2,indices(nw)+1,a,1) dx_old(3) = -en_distance_rescaled_gl(3,indices(nw)+1,a,1) f = a_vector(1, type_nucl_vector(a)+1) * invdenom2 f_old = a_vector(1, type_nucl_vector(a)+1) * invdenom2_old forces_jastrow_quad_en(1,a,nw) = forces_jastrow_quad_en(1,a,nw) + f * dx(1) - f_old * dx_old(1) forces_jastrow_quad_en(2,a,nw) = forces_jastrow_quad_en(2,a,nw) + f * dx(2) - f_old * dx_old(2) forces_jastrow_quad_en(3,a,nw) = forces_jastrow_quad_en(3,a,nw) + f * dx(3) - f_old * dx_old(3) 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 forces_jastrow_quad_en(1,a,nw) = forces_jastrow_quad_en(1,a,nw) + f * x1 * dx(1) - f_old * x1_old * dx_old(1) forces_jastrow_quad_en(2,a,nw) = forces_jastrow_quad_en(2,a,nw) + f * x1 * dx(2) - f_old * x1_old * dx_old(2) forces_jastrow_quad_en(3,a,nw) = forces_jastrow_quad_en(3,a,nw) + f * x1 * dx(3) - f_old * x1_old * dx_old(3) x = x*x1 x_old = x_old*x1_old kf = kf + 1.d0 end do end do end do end function qmckl_compute_forces_jastrow_quad_en_doc #+end_src #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_forces_jastrow_quad_en_doc ( const qmckl_context context, const int64_t num, const int64_t* indices, 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_quad, const double* en_rescaled_quad_gl, double* const forces_jastrow_quad_en ); qmckl_exit_code qmckl_compute_forces_jastrow_quad_en ( const qmckl_context context, const int64_t num, const int64_t* indices, 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_quad, const double* en_rescaled_quad_gl, double* const forces_jastrow_quad_en ); #+end_src #+begin_src c :tangle (eval c) :comments org :exports none qmckl_exit_code qmckl_compute_forces_jastrow_quad_en (const qmckl_context context, const int64_t num, const int64_t* indices, 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_quad, const double* en_rescaled_quad_gl, double* const forces_jastrow_quad_en ) { #ifdef HAVE_HPC return qmckl_compute_forces_jastrow_quad_en_doc #else return qmckl_compute_forces_jastrow_quad_en_doc #endif (context, num, indices, 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_quad, en_rescaled_quad_gl, forces_jastrow_quad_en ); } #+end_src ** Test * Force of delta_p matrix ** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_forces_jastrow_quad_delta_p(qmckl_context context, double* const forces_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_forces_jastrow_quad_delta_p(qmckl_context context, double* const forces_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_forces_jastrow_quad_delta_p(context); if (rc != QMCKL_SUCCESS) return rc; qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int64_t sze = 3 * ctx->quad_point.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_forces_jastrow_delta_p", "Array too small."); } memcpy(forces_delta_p, ctx->quad_point.forces_delta_p, 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_forces_jastrow_quad_delta_p(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_provide_forces_jastrow_quad_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_quad_e(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_een_rescaled_n(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_een_rescaled_e(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_een_rescaled_n_gl(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_een_rescaled_quad_n(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_een_rescaled_quad_e(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_een_rescaled_quad_n_gl(context); if(rc != QMCKL_SUCCESS) return rc; /* Compute if necessary */ if (ctx->quad_point.date > ctx->quad_point.forces_delta_p_date) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = 3 * ctx->quad_point.num * ctx->jastrow_champ.cord_num * (ctx->jastrow_champ.cord_num + 1) * ctx->nucleus.num * ctx->electron.num * sizeof(double); if (mem_info.size > ctx->quad_point.forces_delta_p_maxsize) { if (ctx->quad_point.forces_delta_p != NULL) { rc = qmckl_free(context, ctx->quad_point.forces_delta_p); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_forces_jastrow_delta_p", "Unable to free ctx->quad_point.forces_delta_p"); } ctx->quad_point.forces_delta_p = NULL; } } /* Allocate array */ if (ctx->quad_point.forces_delta_p == NULL) { double* forces_delta_p = (double*) qmckl_malloc(context, mem_info); ctx->quad_point.forces_delta_p_maxsize = mem_info.size; if (forces_delta_p == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_forces_jastrow_delta_p", NULL); } ctx->quad_point.forces_delta_p = forces_delta_p; } rc = qmckl_compute_forces_jastrow_quad_delta_p(context, ctx->quad_point.num, ctx->quad_point.indices, 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->quad_point.een_rescaled_quad_n, ctx->quad_point.een_rescaled_quad_e, ctx->jastrow_champ.een_rescaled_n_gl, ctx->quad_point.een_rescaled_quad_n_gl, ctx->quad_point.forces_delta_p); if (rc != QMCKL_SUCCESS) { return rc; } ctx->quad_point.forces_delta_p_date = ctx->quad_point.date; } return QMCKL_SUCCESS; } #+end_src ** Compute :PROPERTIES: :Name: qmckl_compute_forces_jastrow_quad_delta_p_doc :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_forces_quad_delta_p_args | Variable | Type | In/Out | Description | |----------------------------+---------------------------------------------------------------------+--------+---------------------------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~num~ | ~int64_t~ | in | Index of quad electron | | ~indices~ | ~int64_t[num]~ | in | Indices of quad electrons | | ~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_quad_n~ | ~double[num][0:cord_num][nucl_num]~ | in | Electron-nucleus quad rescaled distances | | ~een_rescaled_quad_e~ | ~double[num][0:cord_num][elec_num]~ | in | Electron-electron quad rescaled distances | | ~een_rescaled_n_gl~ | ~double[walk_num][0:cord_num][nucl_num][4][elec_num]~ | in | Electron-nucleus rescaled distances derivatives | | ~een_rescaled_quad_n_gl~ | ~double[num][0:cord_num][nucl_num][3]~ | in | Electron-nucleus quad rescaled distances derivatives | | ~forces_delta_p~ | ~double[num][0:cord_num-1][0:cord_num][nucl_num][3][elec_num]~ | out | Delta P matrix gradient and Laplacian | |----------------------------+---------------------------------------------------------------------+--------+---------------------------------------------------------| #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer(qmckl_exit_code) function qmckl_compute_forces_jastrow_quad_delta_p_doc( & context, num, indices, walk_num, elec_num, nucl_num, cord_num, & een_rescaled_n, een_rescaled_e, een_rescaled_quad_n, een_rescaled_quad_e, & een_rescaled_n_gl, een_rescaled_quad_n_gl, forces_delta_p) & result(info) bind(C) use, intrinsic :: iso_c_binding use qmckl implicit none integer(qmckl_context), intent(in) :: context integer(c_int64_t) , intent(in), value :: num, walk_num, elec_num, cord_num, nucl_num integer(c_int64_t) , intent(in) :: indices(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_quad_n(nucl_num, 0:cord_num, num) real(c_double) , intent(in) :: een_rescaled_quad_e(elec_num, 0:cord_num, 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_quad_n_gl(3, nucl_num, 0:cord_num, num) real(c_double) , intent(out) :: forces_delta_p(elec_num,3,nucl_num,0:cord_num, 0:cord_num-1, num) double precision :: accu, tmp1, tmp2 double precision :: delta_e(elec_num) integer*8 :: i, a, j, l, k, p, m, n, nw, idx double precision :: tmp integer*8 :: LDA, LDB, LDC 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) then forces_delta_p = 0.d0 return endif do nw=1, num idx = indices(nw)+1 do m=0, cord_num-1 do j = 1, elec_num delta_e(j) = -een_rescaled_quad_e(j,m,nw) + een_rescaled_e(j,idx,m,1) end do do l=1, cord_num do a = 1, nucl_num do k = 1, 3 tmp1 = een_rescaled_n_gl(idx, k, a, l, 1) tmp2 = een_rescaled_quad_n_gl(k,a,l,nw) accu = 0.d0 do j = 1, elec_num forces_delta_p(j,k,a,l,m,nw) = een_rescaled_e(j,idx,m,1) * tmp1 - & een_rescaled_quad_e(j,m,nw) * tmp2 accu = accu + delta_e(j) * een_rescaled_n_gl(j,k,a,l,1) end do forces_delta_p(idx,k,a,l,m,nw) = forces_delta_p(idx,k,a,l,m,nw) + accu end do end do end do end do end do end function qmckl_compute_forces_jastrow_quad_delta_p_doc #+end_src #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_forces_jastrow_quad_delta_p_doc (const qmckl_context context, const int64_t num, const int64_t* indices, 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_quad_n, const double* een_rescaled_quad_e, const double* een_rescaled_n_gl, const double* een_rescaled_quad_n_gl, double* const forces_delta_p ); qmckl_exit_code qmckl_compute_forces_jastrow_quad_delta_p (const qmckl_context context, const int64_t num, const int64_t* indices, 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_quad_n, const double* een_rescaled_quad_e, const double* een_rescaled_n_gl, const double* een_rescaled_quad_n_gl, double* const forces_delta_p ); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_compute_forces_jastrow_quad_delta_p (const qmckl_context context, const int64_t num, const int64_t* indices, 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_quad_n, const double* een_rescaled_quad_e, const double* een_rescaled_n_gl, const double* een_rescaled_quad_n_gl, double* const forces_delta_p ) { #ifdef HAVE_HPC return qmckl_compute_forces_jastrow_quad_delta_p_doc #else return qmckl_compute_forces_jastrow_quad_delta_p_doc #endif (context, num, indices, walk_num, elec_num, nucl_num, cord_num, een_rescaled_n, een_rescaled_e, een_rescaled_quad_n, een_rescaled_quad_e, een_rescaled_n_gl, een_rescaled_quad_n_gl, forces_delta_p); } #+end_src ** Test * Force of quad een Jastrow ** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_forces_jastrow_quad_een(qmckl_context context, double* const forces_jastrow_quad_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_forces_jastrow_quad_een(qmckl_context context, double* const forces_jastrow_quad_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_forces_jastrow_quad_een(context); if (rc != QMCKL_SUCCESS) return rc; qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); const int64_t sze = ctx->quad_point.num * 3 * ctx->nucleus.num; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_forces_jastrow_quad_een", "input array too small"); } memcpy(forces_jastrow_quad_een, ctx->quad_point.forces_jastrow_quad_een, sze * sizeof(double)); return QMCKL_SUCCESS; } #+end_src #+begin_src f90 :tangle (eval fh_func) :comments org :exports none interface integer(qmckl_exit_code) function qmckl_get_forces_jastrow_quad_een (context, & forces_jastrow_quad_een, size_max) bind(C) use, intrinsic :: iso_c_binding import implicit none integer (c_int64_t) , intent(in) , value :: context integer (c_int64_t) , intent(in) , value :: size_max real(c_double), intent(out) :: forces_jastrow_quad_een(size_max) end function qmckl_get_forces_jastrow_quad_een 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_forces_jastrow_quad_een(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :export none qmckl_exit_code qmckl_provide_forces_jastrow_quad_een(qmckl_context context) { qmckl_exit_code rc = QMCKL_SUCCESS; if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return qmckl_failwith( context, QMCKL_INVALID_CONTEXT, "qmckl_provide_forces_jastrow_quad_een", 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_forces_jastrow_quad_een", NULL); } /* Check if en rescaled distance is provided */ rc = qmckl_provide_een_rescaled_n(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if en rescaled distance derivatives is provided */ rc = qmckl_provide_een_rescaled_n_gl(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if en rescaled distance is provided */ rc = qmckl_provide_een_rescaled_quad_n(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if en rescaled distance derivatives is provided */ rc = qmckl_provide_een_rescaled_quad_n_gl(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_tmp_c(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_forces_tmp_c(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_jastrow_champ_quad_delta_p(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_forces_jastrow_quad_delta_p(context); if(rc != QMCKL_SUCCESS) return rc; /* Compute if necessary */ if (ctx->quad_point.date > ctx->quad_point.forces_jastrow_quad_een_date) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = ctx->quad_point.num * 3 * ctx->nucleus.num * sizeof(double); if (mem_info.size > ctx->quad_point.forces_jastrow_quad_een_maxsize) { if (ctx->quad_point.forces_jastrow_quad_een != NULL) { rc = qmckl_free(context, ctx->quad_point.forces_jastrow_quad_een); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_forces_jastrow_quad_een", "Unable to free ctx->quad_point.forces_jastrow_quad_een"); } ctx->quad_point.forces_jastrow_quad_een = NULL; } } /* Allocate array */ if (ctx->quad_point.forces_jastrow_quad_een == NULL) { double* forces_jastrow_quad_een = (double*) qmckl_malloc(context, mem_info); ctx->quad_point.forces_jastrow_quad_een_maxsize = mem_info.size; if (forces_jastrow_quad_een == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_forces_jastrow_quad_een", NULL); } ctx->quad_point.forces_jastrow_quad_een = forces_jastrow_quad_een; } rc = qmckl_compute_forces_jastrow_quad_een(context, ctx->quad_point.num, ctx->quad_point.indices, 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->quad_point.delta_p, ctx->quad_point.forces_delta_p, ctx->jastrow_champ.tmp_c, ctx->forces.forces_tmp_c, ctx->jastrow_champ.een_rescaled_n, ctx->quad_point.een_rescaled_quad_n, ctx->jastrow_champ.een_rescaled_n_gl, ctx->quad_point.een_rescaled_quad_n_gl, ctx->quad_point.forces_jastrow_quad_een); if (rc != QMCKL_SUCCESS) { return rc; } ctx->quad_point.forces_jastrow_quad_een_date = ctx->quad_point.date; } return QMCKL_SUCCESS; } #+end_src ** Compute :PROPERTIES: :Name: qmckl_compute_forces_jastrow_quad_een :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_forces_quad_een_args | Variable | Type | In/Out | Description | |---------------------------+-------------------------------------------+--------+---------------------------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~num~ | ~int64_t~ | in | Index of quad electron | | ~indices~ | ~int64_t[num]~ | in | Indices of quad electrons | | ~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 | | ~delta_p~ | ~double[num][0:cord_num-1][0:cord_num][nucl_num][elec_num]~ | in | quad electron P matrix | | ~forces_delta_p~ | ~double[num][0:cord_num-1][0:cord_num][nucl_num][3][elec_num]~ | in | quad electron P matrix | | ~tmp_c~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]~ | in | quad electron P matrix | | ~forces_tmp_c~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][4][elec_num]~ | in | quad electron P matrix | | ~een_rescaled_n~ | ~double[walk_num][0:cord_num][nucl_num][elec_num]~ | in | Electron-nucleus rescaled distances | | ~een_rescaled_quad_n~ | ~double[num][0:cord_num][nucl_num]~ | in | Electron-nucleus quad rescaled distances | | ~een_rescaled_n_gl~ | ~double[walk_num][0:cord_num][nucl_num][4][elec_num]~ | in | Electron-nucleus rescaled distances derivatives | | ~een_rescaled_quad_n_gl~ | ~double[num][0:cord_num][nucl_num][3]~ | in | Electron-nucleus quad rescaled distances derivatives | | ~forces_jastrow_quad_een~ | ~double[num][nucl_num][3]~ | out | quad electron-nucleus forces | |---------------------------+-------------------------------------------+--------+---------------------------------------------------------| #+begin_src f90 :comments org :tangle (eval f) :noweb yes function qmckl_compute_forces_jastrow_quad_een_doc( & context, num, indices, walk_num, elec_num, nucl_num, cord_num, & dim_c_vector, c_vector_full, lkpm_combined_index, & delta_p, forces_delta_p, tmp_c, forces_tmp_c, & een_rescaled_n, een_rescaled_quad_n, een_rescaled_n_gl, een_rescaled_quad_n_gl, forces_jastrow_quad_een) & bind(C) result(info) use qmckl implicit none integer (qmckl_context), intent(in), value :: context integer (c_int64_t) , intent(in) , value :: num integer (c_int64_t) , intent(in) :: indices(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 :: dim_c_vector integer (c_int64_t) , intent(in) , value :: cord_num integer(c_int64_t) , intent(in) :: lkpm_combined_index(dim_c_vector,4) real(c_double) , intent(in) :: c_vector_full(nucl_num, dim_c_vector) real (c_double ) , intent(in) :: delta_p(elec_num, nucl_num,0:cord_num, 0:cord_num-1, num) real (c_double ) , intent(in) :: forces_delta_p(elec_num, 3, nucl_num,0:cord_num, 0:cord_num-1, num) 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) :: forces_tmp_c(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_quad_n(nucl_num, 0:cord_num, 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_quad_n_gl(3, nucl_num, 0:cord_num, num) real (c_double ) , intent(out) :: forces_jastrow_quad_een(3,nucl_num,num) integer(qmckl_exit_code) :: info double precision :: een_rescaled_delta_n(nucl_num, 0:cord_num), een_rescaled_delta_n_gl(3,nucl_num, 0:cord_num) integer*8 :: i, a, j, l, k, p, m, n, nw, kk, idx double precision :: accu, accu2, cn integer*8 :: LDA, LDB, LDC 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 forces_jastrow_quad_een = 0.0d0 if (cord_num == 0) return do nw =1, num idx = indices(nw)+1 een_rescaled_delta_n(:,:) = een_rescaled_quad_n(:,:,nw) - een_rescaled_n(idx,:,:,1) do kk = 1,3 een_rescaled_delta_n_gl(kk,:,:) = een_rescaled_quad_n_gl(kk,:,:,nw) - een_rescaled_n_gl(idx,kk,:,:,1) end do 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 kk = 1, 3 accu = 0.0d0 do j = 1, elec_num accu = accu - een_rescaled_n_gl(j,kk,a,m,1) * delta_p(j,a,m+l,k,nw) + & een_rescaled_n(j,a,m,1) * forces_delta_p(j,kk,a,m+l,k,nw) end do accu = accu - een_rescaled_delta_n_gl(kk,a,m) * (tmp_c(idx,a,m+l,k,1) + delta_p(idx,a,m+l,k,nw)) + & een_rescaled_delta_n(a,m) * (forces_tmp_c(idx,kk,a,m+l,k,1) + forces_delta_p(idx,kk,a,m+l,k,nw)) forces_jastrow_quad_een(kk,a,nw) = forces_jastrow_quad_een(kk,a,nw) + accu * cn end do end do end do end do end function qmckl_compute_forces_jastrow_quad_een_doc #+end_src #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_forces_jastrow_quad_een_doc ( const qmckl_context context, const int64_t num, const int64_t* indices, 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* delta_p, const double* forces_delta_p, const double* tmp_c, const double* forces_tmp_c, const double* een_rescaled_n, const double* een_rescaled_quad_n, const double* een_rescaled_n_gl, const double* een_rescaled_quad_n_gl, double* const forces_jastrow_quad_een ); qmckl_exit_code qmckl_compute_forces_jastrow_quad_een ( const qmckl_context context, const int64_t num, const int64_t* indices, 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* delta_p, const double* forces_delta_p, const double* tmp_c, const double* forces_tmp_c, const double* een_rescaled_n, const double* een_rescaled_quad_n, const double* een_rescaled_n_gl, const double* een_rescaled_quad_n_gl, double* const forces_jastrow_quad_een ); #+end_src #+begin_src c :tangle (eval c) :comments org :exports none qmckl_exit_code qmckl_compute_forces_jastrow_quad_een (const qmckl_context context, const int64_t num, const int64_t* indices, 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* delta_p, const double* forces_delta_p, const double* tmp_c, const double* forces_tmp_c, const double* een_rescaled_n, const double* een_rescaled_quad_n, const double* een_rescaled_n_gl, const double* een_rescaled_quad_n_gl, double* const forces_jastrow_quad_een ) { #ifdef HAVE_HPC return qmckl_compute_forces_jastrow_quad_een_doc #else return qmckl_compute_forces_jastrow_quad_een_doc #endif (context, num, indices, walk_num, elec_num, nucl_num, cord_num, dim_c_vector, c_vector_full, lkpm_combined_index, delta_p, forces_delta_p, tmp_c, forces_tmp_c, een_rescaled_n, een_rescaled_quad_n, een_rescaled_n_gl, een_rescaled_quad_n_gl, forces_jastrow_quad_een ); } #+end_src ** Test * 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 :noexport: #+begin_src c :tangle (eval c_test) rc = qmckl_context_destroy(context); assert (rc == QMCKL_SUCCESS); return 0; } #+end_src