diff --git a/org/qmckl_jastrow.org b/org/qmckl_jastrow.org index 5830ebb..df09b9c 100644 --- a/org/qmckl_jastrow.org +++ b/org/qmckl_jastrow.org @@ -81,26 +81,26 @@ int main() { | ~double~ | ~aord_vector[aord_num + 1][type_nucl_num]~ | in | Order of a polynomial coefficients | | ~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~ | out | Jastrow factor: electron-electron part | + | ~double~ | ~factor_ee[walk_num]~ | out | Jastrow factor: electron-electron part | | ~double~ | ~factor_ee_date~ | out | Jastrow factor: electron-electron part | - | ~double~ | ~factor_en~ | out | Jastrow factor: electron-nucleus part | + | ~double~ | ~factor_en[walk_num]~ | out | Jastrow factor: electron-nucleus part | | ~double~ | ~factor_en_date~ | out | Jastrow factor: electron-nucleus part | - | ~double~ | ~factor_een~ | out | Jastrow factor: electron-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 | - | ~double~ | ~factor_ee_deriv_e[4][nelec]~ | out | Derivative of the Jastrow factor: electron-electron-nucleus part | - | ~double~ | ~factor_en_deriv_e[4][nelec]~ | out | Derivative of the Jastrow factor: electron-electron-nucleus part | - | ~double~ | ~factor_een_deriv_e[4][nelec]~ | out | Derivative of the 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_en_deriv_e[4][nelec][walk_num]~ | out | Derivative of the Jastrow factor: electron-electron-nucleus part | + | ~double~ | ~factor_een_deriv_e[4][nelec][walk_num]~ | out | Derivative of the Jastrow factor: electron-electron-nucleus part | 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]~ | vector of non-zero coefficients | - | ~double~ | ~dtmp_c[elec_num][4][nucl_num][ncord + 1][ncord]~ | 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 | For H2O we have the following data: @@ -108,6 +108,68 @@ int main() { #+BEGIN_SRC python :results output import numpy as np +elec_num = 10 +nucl_num = 2 +up_num = 5 +down_num = 5 +nucl_coord = [ [0.000000, 0.000000 ], + [0.000000, 0.000000 ], + [2.059801, 2.059801 ] ] + +elec_coord = [[[-0.250655104764153 , 0.503070975550133 , -0.166554344502303], + [-0.587812193472177 , -0.128751981129274 , 0.187773606533075], + [ 1.61335569047166 , -0.615556732874863 , -1.43165470979934 ], + [-4.901239896295210E-003 , -1.120440036458986E-002 , 1.99761909330422 ], + [ 0.766647499681200 , -0.293515395797937 , 3.66454589201239 ], + [-0.127732483187947 , -0.138975497694196 , -8.669850480215846E-002], + [-0.232271834949124 , -1.059321673434182E-002 , -0.504862241464867], + [ 1.09360863531826 , -2.036103063808752E-003 , -2.702796910818986E-002], + [-0.108090166832043 , 0.189161729653261 , 2.15398313919894], + [ 0.397978144318712 , -0.254277292595981 , 2.54553335476344]]]; + +ee_distance_rescaled = [ +[ 0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000, + 0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000, + 0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000, + 0.000000000000000E+000], +[ 0.550227800352402 ,0.000000000000000E+000 ,0.000000000000000E+000, + 0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000, + 0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000, + 0.000000000000000E+000], +[ 0.919155060185168 ,0.937695909123175 ,0.000000000000000E+000, + 0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000, + 0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000, + 0.000000000000000E+000], +[ 0.893325429242815 ,0.851181978173561 ,0.978501685226877 , + 0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000, + 0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000, + 0.000000000000000E+000], +[ 0.982457268305353 ,0.976125002619471 ,0.994349933143149 , + 0.844077311588328 ,0.000000000000000E+000 ,0.000000000000000E+000, + 0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000, + 0.000000000000000E+000], +[ 0.482407528408731 ,0.414816073699124 ,0.894716035479343 , + 0.876540187084407 ,0.978921170036895 ,0.000000000000000E+000, + 0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000, + 0.000000000000000E+000], +[ 0.459541909660400 ,0.545007215761510 ,0.883752955884551 , + 0.918958134888791 ,0.986386936267237 ,0.362209822236419 , + 0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000, + 0.000000000000000E+000], +[ 0.763732576854455 ,0.817282762358449 ,0.801802919535959 , + 0.900089095449775 ,0.975704636491453 ,0.707836537586060 , + 0.755705808346586 ,0.000000000000000E+000 ,0.000000000000000E+000, + 0.000000000000000E+000], +[ 0.904249454052971 ,0.871097965261373 ,0.982717262706270 , + 0.239901207363622 ,0.836519456769083 ,0.896135326270534 , + 0.930694340243023 ,0.917708540815567 ,0.000000000000000E+000, + 0.000000000000000E+000], +[ 0.944400908070716 ,0.922589018494961 ,0.984615718580670 , + 0.514328661540623 ,0.692362267147064 ,0.931894098453677 , + 0.956034127544344 ,0.931221472309472 ,0.540903688625053 , + 0.000000000000000E+000]] + + type_nucl_num = 1 aord_num = 5 bord_num = 5 @@ -167,8 +229,10 @@ lkpm_of_cindex = [[1 , 1 , 2 , 0], [3 , 0 , 5 , 1], [1 , 0 , 5 , 2]] +kappa = 1.0 +kappa_inv = 1.0/kappa #+END_SRC - + #+RESULTS: jastrow_data ** Data structure @@ -976,7 +1040,7 @@ assert(qmckl_nucleus_provided(context)); compute. If it is the case, then the data is recomputed and the current date is stored. -** Asymptotic component for \(J_{ee}\) +** Asymptotic component for \(f_{ee}\) Calculate the asymptotic component ~asymp_jasb~ to be substracted from the final electron-electron jastrow factor \(f_{ee}\). The asymptotic componenet is calculated @@ -1107,13 +1171,11 @@ integer function qmckl_compute_asymp_jasb_f(context, bord_num, bord_vector, resc info = QMCKL_INVALID_CONTEXT return endif - print *,"In Compute 1", asym_one if (bord_num <= 0) then info = QMCKL_INVALID_ARG_2 return endif - print *,"In Compute 2", asym_one asym_one = bord_vector(1) * kappa_inv / (1.0d0 + bord_vector(2) * kappa_inv) asymp_jasb(:) = (/asym_one, 0.5d0 * asym_one/) @@ -1167,14 +1229,12 @@ end function qmckl_compute_asymp_jasb_f #+end_src *** Test + #+name: asymp_jasb #+begin_src python :results output :exports none :noweb yes import numpy as np <> -kappa = 1.0 -kappa_inv = 1.0 - asym_one = bord_vector[0] * kappa_inv / (1.0 + bord_vector[1]*kappa_inv) asymp_jasb = np.array([asym_one, 0.5 * asym_one]) @@ -1189,6 +1249,11 @@ print("asymp_jasb[0] : ", asymp_jasb[0]) print("asymp_jasb[1] : ", asymp_jasb[1]) #+end_src + #+RESULTS: asymp_jasb + : asym_one : 0.6634291325000664 + : asymp_jasb[0] : 1.043287918508297 + : asymp_jasb[1] : 0.7115733522582638 + #+RESULTS: : asym_one : 0.43340325572525706 : asymp_jasb[0] : 0.5323750557252571 @@ -1232,6 +1297,307 @@ assert(fabs(asymp_jasb[1]-0.31567342786262853) < 1.e-12); #+end_src + +** Electron-electron component \(f_{ee}\) + + Calculate the asymptotic component ~asymp_jasb~ to be substracted from the final + electron-electron jastrow factor \(f_{ee}\). The asymptotic componenet is calculated + via the ~bord_vector~ and the electron-electron rescale factor ~rescale_factor_kappa~. + + \[ +f_{ee} = \sum_{i,jjastrow.factor_ee, ctx->electron.walk_num*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_factor_ee(qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_provide_factor_ee(qmckl_context context) +{ + + qmckl_exit_code rc; + + 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 rescaled distance is provided */ + rc = qmckl_provide_ee_distance_rescaled(context); + if(rc != QMCKL_SUCCESS) return rc; + + /* Compute if necessary */ + if (ctx->date > ctx->jastrow.factor_ee_date) { + + /* Allocate array */ + if (ctx->jastrow.factor_ee == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->electron.walk_num * sizeof(double); + double* factor_ee = (double*) qmckl_malloc(context, mem_info); + + if (factor_ee == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_provide_factor_ee", + NULL); + } + ctx->jastrow.factor_ee = factor_ee; + } + + qmckl_exit_code rc = + qmckl_compute_factor_ee(context, + ctx->electron.walk_num, + ctx->electron.num, + ctx->electron.up_num, + ctx->jastrow.bord_num, + ctx->jastrow.bord_vector, + ctx->electron.ee_distance_rescaled, + ctx->jastrow.asymp_jasb, + ctx->jastrow.factor_ee); + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->jastrow.factor_ee_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + #+end_src + +*** Compute + :PROPERTIES: + :Name: qmckl_compute_factor_ee + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_factor_ee_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 | up_num | in | Number of alpha electrons | + | int64_t | bord_num | in | Number of coefficients | + | double | bord_vector[bord_num + 1] | in | List of coefficients | + | double | ee_distance_rescaled[walk_num][elec_num][elec_num] | in | Electron-electron distances | + | double | asymp_jasb[2] | in | Electron-electron distances | + | double | factor_ee[walk_num] | out | Electron-electron distances | + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer function qmckl_compute_factor_ee_f(context, walk_num, elec_num, up_num, bord_num, & + bord_vector, ee_distance_rescaled, asymp_jasb, factor_ee) & + result(info) + use qmckl + implicit none + integer(qmckl_context), intent(in) :: context + integer*8 , intent(in) :: walk_num, elec_num, bord_num, up_num + double precision , intent(in) :: bord_vector(bord_num) + double precision , intent(in) :: ee_distance_rescaled(walk_num, elec_num, elec_num) + double precision , intent(in) :: asymp_jasb(2) + double precision , intent(out) :: factor_ee(walk_num) + + integer*8 :: i, j, p, ipar, nw + double precision :: pow_ser, x, spin_fact, power_ser + + 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 (bord_num <= 0) then + info = QMCKL_INVALID_ARG_4 + return + endif + + factor_ee = 0.0d0 + + do nw =1, walk_num + do j = 1, elec_num + do i = 1, j - 1 + x = ee_distance_rescaled(nw,i,j) + power_ser = 0.0d0 + spin_fact = 1.0d0 + ipar = 1 + + do p = 2, bord_num + x = x * ee_distance_rescaled(nw,i,j) + power_ser = power_ser + bord_vector(p + 1) * x + end do + + if(j .LE. up_num .OR. i .GT. up_num) then + spin_fact = 0.5d0 + ipar = 2 + endif + + factor_ee(nw) = factor_ee(nw) + spin_fact * bord_vector(1) * & + ee_distance_rescaled(nw,i,j) / & + (1.0d0 + bord_vector(2) * & + ee_distance_rescaled(nw,i,j)) & + -asymp_jasb(ipar) + power_ser + + end do + end do + end do + +end function qmckl_compute_factor_ee_f + #+end_src + + #+CALL: generate_c_header(table=qmckl_factor_ee_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + + #+RESULTS: + #+begin_src c :tangle (eval h_func) :comments org + qmckl_exit_code qmckl_compute_factor_ee ( + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t up_num, + const int64_t bord_num, + const double* bord_vector, + const double* ee_distance_rescaled, + const double* asymp_jasb, + double* const factor_ee ); + #+end_src + + + #+CALL: generate_c_interface(table=qmckl_factor_ee_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_factor_ee & + (context, & + walk_num, & + elec_num, & + up_num, & + bord_num, & + bord_vector, & + ee_distance_rescaled, & + asymp_jasb, & + factor_ee) & + 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 :: up_num + integer (c_int64_t) , intent(in) , value :: bord_num + real (c_double ) , intent(in) :: bord_vector(bord_num + 1) + real (c_double ) , intent(in) :: ee_distance_rescaled(elec_num,elec_num,walk_num) + real (c_double ) , intent(in) :: asymp_jasb(2) + real (c_double ) , intent(out) :: factor_ee(walk_num) + + integer(c_int32_t), external :: qmckl_compute_factor_ee_f + info = qmckl_compute_factor_ee_f & + (context, & + walk_num, & + elec_num, & + up_num, & + bord_num, & + bord_vector, & + ee_distance_rescaled, & + asymp_jasb, & + factor_ee) + + end function qmckl_compute_factor_ee + #+end_src + +*** Test + #+begin_src python :results output :exports none :noweb yes +import numpy as np + +<> + +<> + +factor_ee = 0.0 +for i in range(0,elec_num): + for j in range(0,i): + x = ee_distance_rescaled[i][j] + pow_ser = 0.0 + spin_fact = 1.0 + ipar = 0 + + for p in range(1,bord_num): + x = x * ee_distance_rescaled[i][j] + pow_ser = pow_ser + bord_vector[p + 1] * x + + if(i < up_num or j >= up_num): + spin_fact = 0.5 + ipar = 1 + + factor_ee = factor_ee + spin_fact * bord_vector[0] * ee_distance_rescaled[i][j] \ + / (1.0 + bord_vector[1] * ee_distance_rescaled[i][j]) \ + - asymp_jasb[ipar] + pow_ser +print("factor_ee :",factor_ee) + + #+end_src + + #+RESULTS: + : asym_one : 0.43340325572525706 + : asymp_jasb[0] : 0.5323750557252571 + : asymp_jasb[1] : 0.31567342786262853 + : factor_ee : -4.282760865958113 + + + #+begin_src c :tangle (eval c_test) +/* Check if Jastrow is properly initialized */ +assert(qmckl_jastrow_provided(context)); + +double factor_ee[walk_num]; +rc = qmckl_get_jastrow_factor_ee(context, factor_ee); + +// calculate factor_ee +assert(fabs(factor_ee[0]+4.282760865958113) < 1.e-12); + + #+end_src + + * End of files :noexport: #+begin_src c :tangle (eval h_private_type)