diff --git a/org/qmckl_electron.org b/org/qmckl_electron.org index dab2ae6..d125803 100644 --- a/org/qmckl_electron.org +++ b/org/qmckl_electron.org @@ -63,33 +63,29 @@ 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 | - | ~een_distance_rescaled~ | ~double[walk_num][nucl_num][num]~ | Electron-nucleus distances | - | ~een_distance_rescaled_date~ | ~uint64_t~ | Last modification date of the electron-electron distances | - | ~een_distance_rescaled_deriv_e~ | ~double[walk_num][4][nucl_num][num]~ | Electron-electron rescaled distances derivatives | - | ~een_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 | + | ~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 | ** Data structure @@ -104,12 +100,10 @@ typedef struct qmckl_electron_struct { int64_t coord_new_date; int64_t ee_distance_date; int64_t en_distance_date; - int64_t een_distance_date; int64_t ee_distance_rescaled_date; int64_t ee_distance_rescaled_deriv_e_date; int64_t en_distance_rescaled_date; int64_t en_distance_rescaled_deriv_e_date; - int64_t een_distance_rescaled_deriv_e_date; double* coord_new; double* coord_old; double* ee_distance; @@ -118,8 +112,6 @@ typedef struct qmckl_electron_struct { double* ee_distance_rescaled_deriv_e; double* en_distance_rescaled; double* en_distance_rescaled_deriv_e; - double* een_distance_rescaled; - double* een_distance_rescaled_deriv_e; int32_t uninitialized; bool provided; } qmckl_electron_struct; @@ -1559,7 +1551,6 @@ rc = qmckl_get_electron_ee_distance_rescaled_deriv_e(context, ee_distance_rescal #+end_src - ** Electron-nucleus distances *** Get diff --git a/org/qmckl_jastrow.org b/org/qmckl_jastrow.org index e05c37a..340fd36 100644 --- a/org/qmckl_jastrow.org +++ b/org/qmckl_jastrow.org @@ -83,28 +83,32 @@ int main() { | ~double~ | ~bord_vector[bord_num + 1]~ | in | Order of b polynomial coefficients | | ~double~ | ~cord_vector[cord_num][type_nucl_num]~ | in | Order of c polynomial coefficients | | ~double~ | ~factor_ee[walk_num]~ | out | Jastrow factor: electron-electron part | - | ~double~ | ~factor_ee_date~ | out | Jastrow factor: electron-electron part | + | ~int64_t~ | ~factor_ee_date~ | out | Jastrow factor: electron-electron part | | ~double~ | ~factor_en[walk_num]~ | out | Jastrow factor: electron-nucleus part | - | ~double~ | ~factor_en_date~ | out | Jastrow factor: electron-nucleus part | + | ~int64_t~ | ~factor_en_date~ | out | Jastrow factor: electron-nucleus part | | ~double~ | ~factor_een[walk_num]~ | out | Jastrow factor: electron-electron-nucleus part | - | ~double~ | ~factor_een_date~ | out | Jastrow factor: electron-electron-nucleus part | + | ~int64_t~ | ~factor_een_date~ | out | Jastrow factor: electron-electron-nucleus part | | ~double~ | ~factor_ee_deriv_e[4][nelec][walk_num]~ | out | Derivative of the Jastrow factor: electron-electron-nucleus part | - | ~double~ | ~factor_ee_deriv_e_date~ | out | Keep track of the date for the derivative | + | ~int64_t~ | ~factor_ee_deriv_e_date~ | out | Keep track of the date for the derivative | | ~double~ | ~factor_en_deriv_e[4][nelec][walk_num]~ | out | Derivative of the Jastrow factor: electron-electron-nucleus part | - | ~double~ | ~factor_en_deriv_e_date~ | out | Keep track of the date for the en derivative | + | ~int64_t~ | ~factor_en_deriv_e_date~ | out | Keep track of the date for the en derivative | | ~double~ | ~factor_een_deriv_e[4][nelec][walk_num]~ | out | Derivative of the Jastrow factor: electron-electron-nucleus part | - | ~double~ | ~factor_een_deriv_e_date~ | out | Keep track of the date for the een derivative | + | ~int64_t~ | ~factor_een_deriv_e_date~ | out | Keep track of the date for the een derivative | computed data: - |-----------+-------------------------------------------------------------+-------------------------------------------------| - | ~int64_t~ | ~dim_cord_vec~ | Number of unique C coefficients | - | ~double~ | ~asymp_jasb[2]~ | Asymptotic component | - | ~int64_t~ | ~asymp_jasb_date~ | Asymptotic component | - | ~double~ | ~coord_vect_full[dim_cord_vec][nucl_num]~ | vector of non-zero coefficients | - | ~int64_t~ | ~lkpm_of_cindex[4][dim_cord_vec]~ | Transform l,k,p, and m into consecutive indices | - | ~double~ | ~tmp_c[elec_num][nucl_num][ncord + 1][ncord][walk_num]~ | vector of non-zero coefficients | - | ~double~ | ~dtmp_c[elec_num][4][nucl_num][ncord + 1][ncord][walk_num]~ | vector of non-zero coefficients | + |-----------+-------------------------------------------------------------+-------------------------------------------------------------------------------| + | ~int64_t~ | ~dim_cord_vec~ | Number of unique C coefficients | + | ~double~ | ~asymp_jasb[2]~ | Asymptotic component | + | ~int64_t~ | ~asymp_jasb_date~ | Asymptotic component | + | ~double~ | ~coord_vect_full[dim_cord_vec][nucl_num]~ | vector of non-zero coefficients | + | ~int64_t~ | ~lkpm_of_cindex[4][dim_cord_vec]~ | Transform l,k,p, and m into consecutive indices | + | ~double~ | ~tmp_c[elec_num][nucl_num][ncord + 1][ncord][walk_num]~ | vector of non-zero coefficients | + | ~double~ | ~dtmp_c[elec_num][4][nucl_num][ncord + 1][ncord][walk_num]~ | vector of non-zero coefficients | + | ~double~ | ~een_rescaled_e[walk_num][elec_num][elec_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord | + | ~int64_t~ | ~een_rescaled_e_date~ | Keep track of the date of creation | + | ~double~ | ~een_rescaled_n[walk_num][elec_num][nucl_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord | + | ~int64_t~ | ~een_rescaled_n_date~ | Keep track of the date of creation | For H2O we have the following data: @@ -293,6 +297,10 @@ typedef struct qmckl_jastrow_struct{ double * coord_vect_full; double * tmp_c; double * dtmp_c; + double * een_rescaled_e; + double * een_rescaled_n; + int64_t een_rescaled_e_date; + int64_t een_rescaled_n_date; bool provided; char * type; } qmckl_jastrow_struct; @@ -1147,7 +1155,7 @@ assert(qmckl_nucleus_provided(context)); the dependencies are more recent than the date of the data to compute. If it is the case, then the data is recomputed and the current date is stored. - + ** Asymptotic component for \(f_{ee}\) Calculate the asymptotic component ~asymp_jasb~ to be substracted from the final @@ -1411,9 +1419,8 @@ assert(fabs(asymp_jasb[1]-0.31567342786262853) < 1.e-12); #+end_src - ** Electron-electron component \(f_{ee}\) - + Calculate the electron-electron jastrow component ~factor_ee~ using the ~asymp_jasb~ componenet and the electron-electron rescaled distances ~ee_distance_rescaled~. @@ -2829,6 +2836,320 @@ assert(fabs(factor_en_deriv_e[0][3][0]+0.9667363412285741 ) < 1.e-12); #+end_src +** Electron-electron rescaled distances for each order + + ~een_rescaled_e~ stores the table of the rescaled distances between all + pairs of electrons and raised to the power \(p\) defined by ~cord_num~: + + \[ + C_{ij,p} = \left( 1 - \exp{-\kappa C_{ij}} \right)^p + \] + + where \(C_{ij}\) is the matrix of electron-electron distances. + +*** Get + + #+begin_src c :comments org :tangle (eval h_func) :noweb yes +qmckl_exit_code qmckl_get_jastrow_een_rescaled_e(qmckl_context context, double* const distance_rescaled); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_get_jastrow_een_rescaled_e(qmckl_context context, double* const distance_rescaled) +{ + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + + qmckl_exit_code rc; + + rc = qmckl_provide_een_rescaled_e(context); + if (rc != QMCKL_SUCCESS) return rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + size_t sze = ctx->electron.num * ctx->electron.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1); + memcpy(distance_rescaled, ctx->jastrow.een_rescaled_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_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_e(qmckl_context context) +{ + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + /* Check if ee distance is provided */ + qmckl_exit_code rc = qmckl_provide_ee_distance(context); + if(rc != QMCKL_SUCCESS) return rc; + + /* Compute if necessary */ + if (ctx->date > ctx->jastrow.een_rescaled_e_date) { + + /* Allocate array */ + if (ctx->jastrow.een_rescaled_e == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->electron.num * ctx->electron.num * + ctx->electron.walk_num * (ctx->jastrow.cord_num + 1) * sizeof(double); + double* een_rescaled_e = (double*) qmckl_malloc(context, mem_info); + + if (een_rescaled_e == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_een_rescaled_e", + NULL); + } + ctx->jastrow.een_rescaled_e = een_rescaled_e; + } + + qmckl_exit_code rc = + qmckl_compute_een_rescaled_e(context, + ctx->electron.walk_num, + ctx->electron.num, + ctx->jastrow.cord_num, + ctx->electron.rescale_factor_kappa_ee, + ctx->electron.ee_distance, + ctx->jastrow.een_rescaled_e); + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->jastrow.een_rescaled_e_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + #+end_src + +*** Compute + :PROPERTIES: + :Name: qmckl_compute_een_rescaled_e + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_factor_een_rescaled_e_args + | qmckl_context | context | in | Global state | + | int64_t | walk_num | in | Number of walkers | + | int64_t | elec_num | in | Number of electrons | + | int64_t | cord_num | in | Order of polynomials | + | double | rescale_factor_kappa_ee | in | Factor to rescale ee distances | + | double | ee_distance[walk_num][elec_num][elec_num] | in | Electron-electron distances | + | double | een_rescaled_e[walk_num][elec_num][elec_num][0:cord_num] | out | Electron-electron rescaled distances | + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer function qmckl_compute_een_rescaled_e_f(context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, & + ee_distance, een_rescaled_e) & + result(info) + use qmckl + implicit none + integer(qmckl_context), intent(in) :: context + integer*8 , intent(in) :: walk_num + integer*8 , intent(in) :: elec_num + integer*8 , intent(in) :: cord_num + double precision , intent(in) :: rescale_factor_kappa_ee + double precision , intent(in) :: ee_distance(elec_num,elec_num,walk_num) + double precision , intent(out) :: een_rescaled_e(0:cord_num,elec_num,elec_num,walk_num) + double precision,dimension(:,:),allocatable :: een_rescaled_e_ij + double precision :: x + integer*8 :: i, j, k, l, nw + + allocate(een_rescaled_e_ij(elec_num * (elec_num - 1) / 2, cord_num + 1)) + + + info = QMCKL_SUCCESS + + if (context == QMCKL_NULL_CONTEXT) then + info = QMCKL_INVALID_CONTEXT + return + endif + + if (walk_num <= 0) then + info = QMCKL_INVALID_ARG_2 + return + endif + + if (elec_num <= 0) then + info = QMCKL_INVALID_ARG_3 + return + endif + + if (cord_num <= 0) then + info = QMCKL_INVALID_ARG_4 + return + endif + + een_rescaled_e_ij(:, 1) = 1.0d0 + + ! Prepare table of exponentiated distances raised to appropriate power + do nw = 1, walk_num + k = 0 + do j = 1, elec_num + do i = 1, j - 1 + k = k + 1 + een_rescaled_e_ij(k, 2) = dexp(-rescale_factor_kappa_ee * ee_distance(i, j, nw)) + end do + end do + + do l = 2, cord_num + do k = 1, elec_num * (elec_num - 1)/2 + een_rescaled_e_ij(k, l + 1) = een_rescaled_e_ij(k, l + 1 - 1) * een_rescaled_e_ij(k, 1) + end do + end do + + ! prepare the actual een table + een_rescaled_e = 0.0d0 + een_rescaled_e(0, :, :, :) = 1.0d0 + do l = 1, cord_num + k = 0 + do j = 1, elec_num + do i = 1, j - 1 + k = k + 1 + x = een_rescaled_e_ij(k, l + 1) + een_rescaled_e(l, i, j, nw) = x + een_rescaled_e(l, j, i, nw) = x + end do + end do + end do + end do + +end function qmckl_compute_een_rescaled_e_f + #+end_src + + #+CALL: generate_c_header(table=qmckl_factor_een_rescaled_e_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + + #+RESULTS: + #+begin_src c :tangle (eval h_func) :comments org + qmckl_exit_code qmckl_compute_een_rescaled_e ( + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t cord_num, + const double rescale_factor_kappa_ee, + const double* ee_distance, + double* const een_rescaled_e ); + #+end_src + + #+CALL: generate_c_interface(table=qmckl_factor_een_rescaled_e_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + + #+RESULTS: + #+begin_src f90 :tangle (eval f) :comments org :exports none + integer(c_int32_t) function qmckl_compute_een_rescaled_e & + (context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, ee_distance, een_rescaled_e) & + bind(C) result(info) + + use, intrinsic :: iso_c_binding + implicit none + + integer (c_int64_t) , intent(in) , value :: context + integer (c_int64_t) , intent(in) , value :: 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_kappa_ee + real (c_double ) , intent(in) :: ee_distance(elec_num,elec_num,walk_num) + real (c_double ) , intent(out) :: een_rescaled_e(0:cord_num,elec_num,elec_num,walk_num) + + integer(c_int32_t), external :: qmckl_compute_een_rescaled_e_f + info = qmckl_compute_een_rescaled_e_f & + (context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, ee_distance, een_rescaled_e) + + end function qmckl_compute_een_rescaled_e + #+end_src + +*** Test + + #+begin_src python :results output :exports none :noweb yes +import numpy as np + +<> + +elec_coord = np.array(elec_coord)[0] +elec_dist = np.zeros(shape=(elec_num, elec_num),dtype=float) +for i in range(elec_num): + for j in range(elec_num): + elec_dist[i, j] = np.linalg.norm(elec_coord[i] - elec_coord[j]) + +kappa = 1.0 + +een_rescaled_e_ij = np.zeros(shape=(elec_num * (elec_num - 1)//2, cord_num+1), dtype=float) +een_rescaled_e_ij[:,0] = 1.0 + +k = 0 +for j in range(elec_num): + for i in range(j - 1): + een_rescaled_e_ij[k, 1] = np.exp(-kappa * elec_dist[i, j]) + k = k + 1 + +for l in range(2, cord_num + 1): + for k in range(elec_num * (elec_num - 1)//2): + een_rescaled_e_ij[k, l] = een_rescaled_e_ij[k, l - 1] * een_rescaled_e_ij[k, 1] + +een_rescaled_e = np.zeros(shape=(elec_num, elec_num, cord_num + 1), dtype=float) +een_rescaled_e[:,:,0] = 1.0 + +for l in range(1,cord_num+1): + k = 0 + for j in range(elec_num): + for i in range(j - 1): + x = een_rescaled_e_ij[k, l] + een_rescaled_e[i, j, l] = x + een_rescaled_e[j, i, l] = x + k = k + 1 + +print(" een_rescaled_e[0, 2, 1] = ",een_rescaled_e[0, 2, 1]) +print(" een_rescaled_e[0, 3, 1] = ",een_rescaled_e[0, 3, 1]) +print(" een_rescaled_e[0, 4, 1] = ",een_rescaled_e[0, 4, 1]) +print(" een_rescaled_e[1, 3, 2] = ",een_rescaled_e[1, 3, 2]) +print(" een_rescaled_e[1, 4, 2] = ",een_rescaled_e[1, 4, 2]) +print(" een_rescaled_e[1, 5, 2] = ",een_rescaled_e[1, 5, 2]) + #+end_src + + #+RESULTS: + : een_rescaled_e[0, 2, 1] = 0.08084493981483197 + : een_rescaled_e[0, 3, 1] = 0.1066745707571846 + : een_rescaled_e[0, 4, 1] = 0.01754273169464735 + : een_rescaled_e[1, 3, 2] = 0.02214680362033448 + : een_rescaled_e[1, 4, 2] = 0.0005700154999202759 + : een_rescaled_e[1, 5, 2] = 0.3424402276009091 + + #+begin_src c :tangle (eval c_test) +assert(qmckl_electron_provided(context)); + + +double een_rescaled_e[walk_num][elec_num][elec_num][(cord_num + 1)]; +rc = qmckl_get_jastrow_een_rescaled_e(context, &(een_rescaled_e[0][0][0][0])); + +// value of (0,2,1) +//printf("%10.15f = \n", een_rescaled_e[0][0][2][1]); +//printf("%10.15f = \n", een_rescaled_e[0][0][3][1]); +//printf("%10.15f = \n", een_rescaled_e[0][0][4][1]); +//printf("%10.15f = \n", een_rescaled_e[0][1][3][2]); +//printf("%10.15f = \n", een_rescaled_e[0][1][4][2]); +//printf("%10.15f = \n", een_rescaled_e[0][1][5][2]); +//assert(fabs(een_rescaled_e[0][0][2][1]-) < 1.e-12); +//assert(fabs(een_rescaled_e[0][0][3][1]-) < 1.e-12); +//assert(fabs(een_rescaled_e[0][0][4][1]-) < 1.e-12); +//assert(fabs(een_rescaled_e[0][1][3][2]-) < 1.e-12); +//assert(fabs(een_rescaled_e[0][1][4][2]-) < 1.e-12); +//assert(fabs(een_rescaled_e[0][1][5][2]-) < 1.e-12); + + #+end_src + + * End of files :noexport: #+begin_src c :tangle (eval h_private_type)