diff --git a/org/qmckl_electron.org b/org/qmckl_electron.org index 9f2f89f..832ac97 100644 --- a/org/qmckl_electron.org +++ b/org/qmckl_electron.org @@ -63,29 +63,36 @@ int main() { The following data stored in the context: - | ~uninitialized~ | ~int32_t~ | Keeps bit set for uninitialized data | - | ~num~ | ~int64_t~ | Total number of electrons | - | ~up_num~ | ~int64_t~ | Number of up-spin electrons | - | ~down_num~ | ~int64_t~ | Number of down-spin electrons | - | ~walk_num~ | ~int64_t~ | Number of walkers | - | ~rescale_factor_kappa_ee~ | ~double~ | The distance scaling factor | - | ~rescale_factor_kappa_en~ | ~double~ | The distance scaling factor | - | ~provided~ | ~bool~ | If true, ~electron~ is valid | - | ~coord_new~ | ~double[walk_num][3][num]~ | New set of electron coordinates | - | ~coord_old~ | ~double[walk_num][3][num]~ | Old set of electron coordinates | - | ~coord_new_date~ | ~uint64_t~ | Last modification date of the coordinates | - | ~ee_distance~ | ~double[walk_num][num][num]~ | Electron-electron distances | - | ~ee_distance_date~ | ~uint64_t~ | Last modification date of the electron-electron distances | - | ~en_distance~ | ~double[walk_num][nucl_num][num]~ | Electron-nucleus distances | - | ~en_distance_date~ | ~uint64_t~ | Last modification date of the electron-electron distances | - | ~ee_distance_rescaled~ | ~double[walk_num][num][num]~ | Electron-electron rescaled distances | - | ~ee_distance_rescaled_date~ | ~uint64_t~ | Last modification date of the electron-electron distances | - | ~ee_distance_rescaled_deriv_e~ | ~double[walk_num][4][num][num]~ | Electron-electron rescaled distances derivatives | - | ~ee_distance_rescaled_deriv_e_date~ | ~uint64_t~ | Last modification date of the electron-electron distance derivatives | - | ~en_distance_rescaled~ | ~double[walk_num][nucl_num][num]~ | Electron-nucleus distances | - | ~en_distance_rescaled_date~ | ~uint64_t~ | Last modification date of the electron-electron distances | - | ~en_distance_rescaled_deriv_e~ | ~double[walk_num][4][nucl_num][num]~ | Electron-electron rescaled distances derivatives | - | ~en_distance_rescaled_deriv_e_date~ | ~uint64_t~ | Last modification date of the electron-electron distance derivatives | + |-------------------------------------+--------------------------------------+----------------------------------------------------------------------| + | ~uninitialized~ | ~int32_t~ | Keeps bit set for uninitialized data | + | ~num~ | ~int64_t~ | Total number of electrons | + | ~up_num~ | ~int64_t~ | Number of up-spin electrons | + | ~down_num~ | ~int64_t~ | Number of down-spin electrons | + | ~walk_num~ | ~int64_t~ | Number of walkers | + | ~rescale_factor_kappa_ee~ | ~double~ | The distance scaling factor | + | ~rescale_factor_kappa_en~ | ~double~ | The distance scaling factor | + | ~provided~ | ~bool~ | If true, ~electron~ is valid | + | ~coord_new~ | ~double[walk_num][3][num]~ | New set of electron coordinates | + | ~coord_old~ | ~double[walk_num][3][num]~ | Old set of electron coordinates | + | ~coord_new_date~ | ~uint64_t~ | Last modification date of the coordinates | + + Computed data: + + |-------------------------------------+--------------------------------------+----------------------------------------------------------------------| + | ~ee_distance~ | ~double[walk_num][num][num]~ | Electron-electron distances | + | ~ee_distance_date~ | ~uint64_t~ | Last modification date of the electron-electron distances | + | ~en_distance~ | ~double[walk_num][nucl_num][num]~ | Electron-nucleus distances | + | ~en_distance_date~ | ~uint64_t~ | Last modification date of the electron-electron distances | + | ~ee_distance_rescaled~ | ~double[walk_num][num][num]~ | Electron-electron rescaled distances | + | ~ee_distance_rescaled_date~ | ~uint64_t~ | Last modification date of the electron-electron distances | + | ~ee_distance_rescaled_deriv_e~ | ~double[walk_num][4][num][num]~ | Electron-electron rescaled distances derivatives | + | ~ee_distance_rescaled_deriv_e_date~ | ~uint64_t~ | Last modification date of the electron-electron distance derivatives | + | ~ee_pot~ | ~double[walk_num][4][num][num]~ | Electron-electron rescaled distances derivatives | + | ~ee_pot_date~ | ~uint64_t~ | Last modification date of the electron-electron distance derivatives | + | ~en_distance_rescaled~ | ~double[walk_num][nucl_num][num]~ | Electron-nucleus distances | + | ~en_distance_rescaled_date~ | ~uint64_t~ | Last modification date of the electron-electron distances | + | ~en_distance_rescaled_deriv_e~ | ~double[walk_num][4][nucl_num][num]~ | Electron-electron rescaled distances derivatives | + | ~en_distance_rescaled_deriv_e_date~ | ~uint64_t~ | Last modification date of the electron-electron distance derivatives | ** Data structure @@ -100,6 +107,7 @@ typedef struct qmckl_electron_struct { int64_t coord_new_date; int64_t ee_distance_date; int64_t en_distance_date; + int64_t ee_pot_date; int64_t ee_distance_rescaled_date; int64_t ee_distance_rescaled_deriv_e_date; int64_t en_distance_rescaled_date; @@ -108,6 +116,7 @@ typedef struct qmckl_electron_struct { double* coord_old; double* ee_distance; double* en_distance; + double* ee_pot; double* ee_distance_rescaled; double* ee_distance_rescaled_deriv_e; double* en_distance_rescaled; @@ -1378,7 +1387,7 @@ qmckl_exit_code qmckl_get_electron_ee_distance_rescaled_deriv_e(qmckl_context co #+end_src *** Provide :noexport: - + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_provide_ee_distance_rescaled_deriv_e(qmckl_context context); #+end_src @@ -1565,6 +1574,194 @@ rc = qmckl_get_electron_ee_distance_rescaled_deriv_e(context, ee_distance_rescal #+end_src +** Electron-electron potential + + ~ee_pot~ calculates the ~ee~ potential energy. + + \[ + \mathcal{V}_{ee} = \sum_{i=1}^{N_e}\sum_{j>i}^{N_e}\frac{1}{r_{ij}} + \] + + where \(\mathcal{V}_{ee}\) is the ~ee~ potential and \[r_{ij}\] the ~ee~ + distance. + +*** Get + + #+begin_src c :comments org :tangle (eval h_func) :noweb yes +qmckl_exit_code qmckl_get_electron_ee_potential(qmckl_context context, double* const ee_pot); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_get_electron_ee_potential(qmckl_context context, double* const ee_pot) +{ + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + + qmckl_exit_code rc; + + rc = qmckl_provide_ee_potential(context); + if (rc != QMCKL_SUCCESS) return rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + size_t sze = ctx->electron.walk_num * sizeof(double); + memcpy(ee_pot, ctx->electron.ee_pot, sze); + + 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_potential(qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_provide_ee_potential(qmckl_context context) +{ + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + if (!ctx->electron.provided) return QMCKL_NOT_PROVIDED; + + /* Compute if necessary */ + if (ctx->electron.coord_new_date > ctx->electron.ee_pot_date) { + + /* Allocate array */ + if (ctx->electron.ee_pot == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->electron.walk_num * sizeof(double); + double* ee_pot = (double*) qmckl_malloc(context, mem_info); + + if (ee_pot == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_ee_potential", + NULL); + } + ctx->electron.ee_pot = ee_pot; + } + + qmckl_exit_code rc = + qmckl_compute_ee_potential(context, + ctx->electron.num, + ctx->electron.walk_num, + ctx->electron.ee_distance, + ctx->electron.ee_pot); + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->electron.ee_pot_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + #+end_src + +*** Compute + :PROPERTIES: + :Name: qmckl_compute_ee_potential + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_ee_potential_args + | qmckl_context | context | in | Global state | + | int64_t | elec_num | in | Number of electrons | + | int64_t | walk_num | in | Number of walkers | + | double | ee_distance[walk_num][elec_num][elec_num] | in | Electron-electron rescaled distances | + | double | ee_pot[walk_num] | out | Electron-electron potential | + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer function qmckl_compute_ee_potential_f(context, elec_num, walk_num, & + ee_distance, ee_pot) & + result(info) + use qmckl + implicit none + integer(qmckl_context), intent(in) :: context + integer*8 , intent(in) :: elec_num + integer*8 , intent(in) :: walk_num + double precision , intent(in) :: ee_distance(elec_num,elec_num,walk_num) + double precision , intent(out) :: ee_pot(walk_num) + + integer*8 :: nw, i, j + + 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_pot = 0.0d0 + do nw=1,walk_num + do j=1,elec_num + do i=1,j-1 + ee_pot(nw) = ee_pot(nw) + 1.0d0/(ee_distance(i,j,nw)) + end do + end do + end do + +end function qmckl_compute_ee_potential_f + #+end_src + + #+CALL: generate_c_header(table=qmckl_ee_potential_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + + #+RESULTS: + #+begin_src c :tangle (eval h_func) :comments org + qmckl_exit_code qmckl_compute_ee_potential ( + const qmckl_context context, + const int64_t elec_num, + const int64_t walk_num, + const double* ee_distance, + double* const ee_pot ); + #+end_src + + #+CALL: generate_c_interface(table=qmckl_ee_potential_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + + #+RESULTS: + #+begin_src f90 :tangle (eval f) :comments org :exports none + integer(c_int32_t) function qmckl_compute_ee_potential & + (context, elec_num, walk_num, ee_distance, ee_pot) & + bind(C) result(info) + + use, intrinsic :: iso_c_binding + implicit none + + integer (c_int64_t) , intent(in) , value :: context + integer (c_int64_t) , intent(in) , value :: elec_num + integer (c_int64_t) , intent(in) , value :: walk_num + real (c_double ) , intent(in) :: ee_distance(elec_num,elec_num,walk_num) + real (c_double ) , intent(out) :: ee_pot(walk_num) + + integer(c_int32_t), external :: qmckl_compute_ee_potential_f + info = qmckl_compute_ee_potential_f & + (context, elec_num, walk_num, ee_distance, ee_pot) + + end function qmckl_compute_ee_potential + #+end_src + +*** Test ** Electron-nucleus distances *** Get diff --git a/org/qmckl_nucleus.org b/org/qmckl_nucleus.org index bdd8cb8..82afe85 100644 --- a/org/qmckl_nucleus.org +++ b/org/qmckl_nucleus.org @@ -62,18 +62,26 @@ int main() { The following data stored in the context: - | ~uninitialized~ | int32_t | Keeps bit set for uninitialized data | - | ~num~ | int64_t | Total number of nuclei | - | ~provided~ | bool | If true, ~nucleus~ is valid | - | ~charge~ | double[num] | Nuclear charges | - | ~coord~ | double[3][num] | Nuclear coordinates, in transposed format | - | ~nn_distance~ | double[num][num] | Nucleus-nucleus distances | - | ~nn_distance_date~ | int64_t | Date when Nucleus-nucleus distances were computed | - | ~nn_distance_rescaled~ | double[num][num] | Nucleus-nucleus rescaled distances | - | ~nn_distance_rescaled_date~ | int64_t | Date when Nucleus-nucleus rescaled distances were computed | - | ~repulsion~ | double | Nuclear repulsion energy | - | ~repulsion_date~ | int64_t | Date when the nuclear repulsion energy was computed | - | ~rescale_factor_kappa~ | double | The distance scaling factor | + |------------------------+----------------+-------------------------------------------| + | ~uninitialized~ | int32_t | Keeps bit set for uninitialized data | + | ~num~ | int64_t | Total number of nuclei | + | ~provided~ | bool | If true, ~nucleus~ is valid | + | ~charge~ | double[num] | Nuclear charges | + | ~coord~ | double[3][num] | Nuclear coordinates, in transposed format | + | ~coord_date~ | int64_t | Nuclear coordinates, date if modified | + | ~rescale_factor_kappa~ | double | The distance scaling factor | + + Computed data: + + |-----------------------------+------------------+--------------------------------------------------------------| + | ~nn_distance~ | double[num][num] | Nucleus-nucleus distances | + | ~nn_distance_date~ | int64_t | Date when Nucleus-nucleus distances were computed | + | ~nn_distance_rescaled~ | double[num][num] | Nucleus-nucleus rescaled distances | + | ~nn_distance_rescaled_date~ | int64_t | Date when Nucleus-nucleus rescaled distances were computed | + | ~repulsion~ | double | Nuclear repulsion energy | + | ~repulsion_date~ | int64_t | Date when the nuclear repulsion energy was computed | + | ~en_pot~ | double | Electron-nucleus potential energy | + | ~en_pot_date~ | int64_t | Date when the electron-nucleus potential energy was computed | ** Data structure @@ -83,6 +91,7 @@ typedef struct qmckl_nucleus_struct { int64_t repulsion_date; int64_t nn_distance_date; int64_t nn_distance_rescaled_date; + int64_t coord_date; double* coord; double* charge; double* nn_distance; @@ -91,6 +100,8 @@ typedef struct qmckl_nucleus_struct { double rescale_factor_kappa; int32_t uninitialized; bool provided; + int64_t en_pot_date; + double* en_pot; } qmckl_nucleus_struct; #+end_src @@ -125,8 +136,6 @@ qmckl_exit_code qmckl_init_nucleus(qmckl_context context) { } #+end_src - - ** Access functions #+begin_src c :comments org :tangle (eval h_func) :exports none @@ -785,7 +794,6 @@ assert(fabs(distance[1]-2.070304721365169) < 1.e-12); #+end_src - ** Nucleus-nucleus rescaled distances *** Get @@ -1132,6 +1140,205 @@ assert(rep - 318.2309879436158 < 1.e-10); #+end_src +** Electron-nucleus potential + ~en_potential~ stores the ~en~ potential energy + + \[ + \mathcal{V}_{en} = \sum_{i=1}^{N_e}\sum_{A=1}^{N_n}\frac{Z_A}{r_{iA}} + \] + + where \(\mathcal{V}_{en}\) is the ~en~ potential, \[r_{iA}\] the ~en~ + distance and \[Z_A\] is the nuclear charge. + +*** Get + + #+begin_src c :comments org :tangle (eval h_func) :noweb yes +qmckl_exit_code qmckl_get_electron_en_potential(qmckl_context context, double* const en_pot); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_get_electron_en_potential(qmckl_context context, double* const en_pot) +{ + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + + qmckl_exit_code rc; + + rc = qmckl_provide_en_potential(context); + if (rc != QMCKL_SUCCESS) return rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + size_t sze = ctx->electron.walk_num * sizeof(double); + memcpy(en_pot, ctx->nucleus.en_pot, sze); + + 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_potential(qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_provide_en_potential(qmckl_context context) +{ + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + if (!ctx->electron.provided) return QMCKL_NOT_PROVIDED; + if (!ctx->nucleus.provided) return QMCKL_NOT_PROVIDED; + + /* Compute if necessary */ + if (ctx->nucleus.coord_date > ctx->nucleus.en_pot_date) { + + /* Allocate array */ + if (ctx->nucleus.en_pot == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->electron.walk_num * sizeof(double); + double* en_pot = (double*) qmckl_malloc(context, mem_info); + + if (en_pot == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_en_potential", + NULL); + } + ctx->nucleus.en_pot = en_pot; + } + + qmckl_exit_code rc = + qmckl_compute_en_potential(context, + ctx->electron.num, + ctx->nucleus.num, + ctx->electron.walk_num, + ctx->nucleus.charge, + ctx->electron.en_distance, + ctx->nucleus.en_pot); + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->nucleus.en_pot_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + #+end_src + +*** Compute + :PROPERTIES: + :Name: qmckl_compute_en_potential + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_en_potential_args + | qmckl_context | context | in | Global state | + | int64_t | elec_num | in | Number of electrons | + | int64_t | nucl_num | in | Number of nucleii | + | int64_t | walk_num | in | Number of walkers | + | double | charge[nucl_num] | in | charge of nucleus | + | double | en_distance[walk_num][nucl_num][elec_num] | in | Electron-electron rescaled distances | + | double | en_pot[walk_num] | out | Electron-electron potential | + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer function qmckl_compute_en_potential_f(context, elec_num, nucl_num, walk_num, & + charge, en_distance, en_pot) & + result(info) + use qmckl + implicit none + integer(qmckl_context), intent(in) :: context + integer*8 , intent(in) :: elec_num + integer*8 , intent(in) :: nucl_num + integer*8 , intent(in) :: walk_num + double precision , intent(in) :: charge(nucl_num) + double precision , intent(in) :: en_distance(elec_num,nucl_num,walk_num) + double precision , intent(out) :: en_pot(walk_num) + + integer*8 :: nw, i, j + + 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 + + en_pot = 0.0d0 + do nw=1,walk_num + do j=1,nucl_num + do i=1,elec_num + en_pot(nw) = en_pot(nw) + charge(j)/(en_distance(i,j,nw)) + end do + end do + end do + +end function qmckl_compute_en_potential_f + #+end_src + + #+CALL: generate_c_header(table=qmckl_en_potential_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + + #+RESULTS: + #+begin_src c :tangle (eval h_func) :comments org + qmckl_exit_code qmckl_compute_en_potential ( + const qmckl_context context, + const int64_t elec_num, + const int64_t nucl_num, + const int64_t walk_num, + const double* charge, + const double* en_distance, + double* const en_pot ); + #+end_src + + #+CALL: generate_c_interface(table=qmckl_en_potential_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + + #+RESULTS: + #+begin_src f90 :tangle (eval f) :comments org :exports none + integer(c_int32_t) function qmckl_compute_en_potential & + (context, elec_num, nucl_num, walk_num, charge, en_distance, en_pot) & + bind(C) result(info) + + use, intrinsic :: iso_c_binding + implicit none + + integer (c_int64_t) , intent(in) , value :: context + integer (c_int64_t) , intent(in) , value :: elec_num + integer (c_int64_t) , intent(in) , value :: nucl_num + integer (c_int64_t) , intent(in) , value :: walk_num + real (c_double ) , intent(in) :: charge(nucl_num) + real (c_double ) , intent(in) :: en_distance(elec_num,nucl_num,walk_num) + real (c_double ) , intent(out) :: en_pot(walk_num) + + integer(c_int32_t), external :: qmckl_compute_en_potential_f + info = qmckl_compute_en_potential_f & + (context, elec_num, nucl_num, walk_num, charge, en_distance, en_pot) + + end function qmckl_compute_en_potential + #+end_src + +*** Test + * End of files :noexport: #+begin_src c :tangle (eval h_private_type)