diff --git a/org/qmckl_jastrow.org b/org/qmckl_jastrow.org index 3b1fcd6..4c37e1a 100644 --- a/org/qmckl_jastrow.org +++ b/org/qmckl_jastrow.org @@ -18,9 +18,10 @@ $J_{\text{eN}}$ contains electron-nucleus terms: \[ - J_{\text{eN}}(\mathbf{r},\mathbf{R}) = \sum_{i=1}^{N_\text{elec}} \sum_{\alpha=1}^{N_\text{nucl}} + J_{\text{eN}}(\mathbf{r},\mathbf{R}) = + \sum_{\alpha=1}^{N_\text{nucl}} \sum_{i=1}^{N_\text{elec}} \frac{a_{1,\alpha}\, g_\alpha(R_{i\alpha})}{1+a_{2,\alpha}\, g_\alpha(R_{i\alpha})} + - \sum_{p=2}^{N_\text{ord}^a} a_{p+1,\alpha}\, [g_\alpha(R_{i\alpha})]^p - J_{eN}^\infty + \sum_{p=2}^{N_\text{ord}^a} a_{p+1,\alpha}\, [g_\alpha(R_{i\alpha})]^p - J_{eN}^{\infty \alpha} \] $J_{\text{ee}}$ contains electron-electron terms: @@ -166,6 +167,8 @@ int main() { |-------------------------------------+-------------------------------------------------------------------+---------------------------------------------------------------------------------------------------------| | ~dim_c_vector~ | ~int64_t~ | Number of unique C coefficients | | ~dim_c_vector_date~ | ~uint64_t~ | Number of unique C coefficients | + | ~asymp_jasa~ | ~double[nucl_num]~ | Asymptotic component | + | ~asymp_jasa_date~ | ~uint64_t~ | Ladt modification of the asymptotic component | | ~asymp_jasb~ | ~double[2]~ | Asymptotic component (up- or down-spin) | | ~asymp_jasb_date~ | ~uint64_t~ | Ladt modification of the asymptotic component | | ~c_vector_full~ | ~double[dim_c_vector][nucl_num]~ | vector of non-zero coefficients | @@ -202,7 +205,7 @@ nucl_coord = np.array([ [0.000000, 0.000000 ], [0.000000, 0.000000 ], [0.000000, 2.059801 ] ]) -elec_coord = [[[-0.250655104764153 , 0.503070975550133 , -0.166554344502303], +elec_coord = np.array( [[[-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 ], @@ -211,9 +214,9 @@ elec_coord = [[[-0.250655104764153 , 0.503070975550133 , -0.16655 [-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]]]; + [ 0.397978144318712 , -0.254277292595981 , 2.54553335476344]]]) -ee_distance_rescaled = [ +ee_distance_rescaled = np.array([ [ 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, @@ -253,7 +256,7 @@ ee_distance_rescaled = [ [ 0.944400908070716, 0.922589018494961, 0.984615718580670, 0.514328661540623, 0.692362267147064, 0.931894098453677, 0.956034127544344, 0.931221472309472, 0.540903688625053, - 0.000000000000000 ]] + 0.000000000000000 ]]) en_distance_rescaled = np.transpose(np.array([ [ 0.443570948411811 , 0.467602196999105 , 0.893870160799932 , @@ -276,16 +279,16 @@ bord_num = 5 cord_num = 5 dim_c_vector= 23 type_nucl_vector = [ 1, 1] -a_vector = [ +a_vector = np.array([ [0.000000000000000E+000], [0.000000000000000E+000], [-0.380512000000000E+000], [-0.157996000000000E+000], [-3.155800000000000E-002], -[2.151200000000000E-002]] +[2.151200000000000E-002]]) -b_vector = [ 0.500000000000000E-000, 0.153660000000000E-000, 6.722620000000000E-002, - 2.157000000000000E-002, 7.309600000000000E-003, 2.866000000000000E-003] +b_vector =np.array( [ 0.500000000000000E-000, 0.153660000000000E-000, 6.722620000000000E-002, + 2.157000000000000E-002, 7.309600000000000E-003, 2.866000000000000E-003]) c_vector = [ 0.571702000000000E-000, -0.514253000000000E-000, -0.513043000000000E-000, 9.486000000000000E-003, -4.205000000000000E-003, 0.426325800000000E-000, 8.288150000000000E-002, 5.118600000000000E-003, -2.997800000000000E-003, @@ -349,6 +352,7 @@ typedef struct qmckl_jastrow_struct{ int64_t bord_num; int64_t cord_num; int64_t type_nucl_num; + uint64_t asymp_jasa_date; uint64_t asymp_jasb_date; uint64_t tmp_c_date; uint64_t dtmp_c_date; @@ -364,6 +368,7 @@ typedef struct qmckl_jastrow_struct{ double * a_vector; double * b_vector; double * c_vector; + double * asymp_jasa; double * asymp_jasb; double * factor_ee; double * factor_en; @@ -1696,6 +1701,22 @@ qmckl_get_jastrow_asymp_jasb(qmckl_context context, } #+end_src +**** Fortran interface + +#+begin_src f90 :tangle (eval fh_func) :comments org +interface + integer(qmckl_exit_code) function qmckl_get_jastrow_asymp_jasb(context, & + asymp_jasb, 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 + double precision, intent(out) :: asymp_jasb(size_max) + end function qmckl_get_jastrow_asymp_jasb +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_asymp_jasb(qmckl_context context); @@ -2957,6 +2978,277 @@ assert(fabs(factor_ee_deriv_e[0][2][0]-0.073267755223968 ) < 1.e-12); assert(fabs(factor_ee_deriv_e[0][3][0]-1.5111672803213185 ) < 1.e-12); #+end_src +** Asymptotic component for \(J_{eN}\) + + Calculate the asymptotic component ~asymp_jasa~ to be substracted from the final + electron-nucleus jastrow factor \(J_{\text{eN}}\). The asymptotic component is calculated + via the ~a_vector~ and the electron-nucleus rescale factors ~rescale_factor_en~. + + \[ + J_{\text{en}}^{\infty \alpha} = \frac{a_1 \kappa_\alpha^{-1}}{1 + a_2 \kappa_\alpha^{-1}} + \] + +*** Get + #+begin_src c :comments org :tangle (eval h_func) :noweb yes +qmckl_exit_code +qmckl_get_jastrow_asymp_jasa(qmckl_context context, + double* const asymp_jasa, + 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_asymp_jasa(qmckl_context context, + double* const asymp_jasa, + const int64_t size_max) +{ + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return qmckl_failwith( context, + QMCKL_INVALID_CONTEXT, + "qmckl_get_jastrow_asymp_jasa", + NULL); + } + + qmckl_exit_code rc; + + rc = qmckl_provide_jastrow_asymp_jasa(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; + if (size_max < sze) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_get_jastrow_asymp_jasa", + "Array too small. Expected nucleus.num"); + } + memcpy(asymp_jasa, ctx->jastrow.asymp_jasa, sze * sizeof(double)); + + return QMCKL_SUCCESS; +} + #+end_src + +**** Fortran interface + +#+begin_src f90 :tangle (eval fh_func) :comments org +interface + integer(qmckl_exit_code) function qmckl_get_jastrow_asymp_jasa(context, & + asymp_jasa, 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 + double precision, intent(out) :: asymp_jasa(size_max) + end function qmckl_get_jastrow_asymp_jasa +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_asymp_jasa(qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_provide_jastrow_asymp_jasa(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_asymp_jasa", + NULL); + } + + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; + assert (ctx != NULL); + + if (!ctx->jastrow.provided) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_provide_jastrow_asymp_jasa", + NULL); + } + + /* Compute if necessary */ + if (ctx->date > ctx->jastrow.asymp_jasa_date) { + + /* Allocate array */ + if (ctx->jastrow.asymp_jasa == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->nucleus.num * sizeof(double); + double* asymp_jasa = (double*) qmckl_malloc(context, mem_info); + + if (asymp_jasa == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_asymp_jasa", + NULL); + } + ctx->jastrow.asymp_jasa = asymp_jasa; + } + + rc = qmckl_compute_jastrow_asymp_jasa(context, + ctx->jastrow.aord_num, + ctx->jastrow.a_vector, + ctx->jastrow.rescale_factor_en, + ctx->nucleus.num, + ctx->jastrow.asymp_jasa); + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->jastrow.asymp_jasa_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + #+end_src + +*** Compute + :PROPERTIES: + :Name: qmckl_compute_jastrow_asymp_jasa + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_asymp_jasa_args + | Variable | Type | In/Out | Description | + |---------------------+----------------------+--------+----------------------------| + | ~context~ | ~qmckl_context~ | in | Global state | + | ~aord_num~ | ~int64_t~ | in | Order of the polynomial | + | ~a_vector~ | ~double[aord_num+1]~ | in | Values of a | + | ~rescale_factor_en~ | ~double~ | in | Electron nucleus distances | + | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | + | ~asymp_jasa~ | ~double[nucl_num]~ | out | Asymptotic value | + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer function qmckl_compute_jastrow_asymp_jasa_f(context, aord_num, a_vector, & + rescale_factor_en, nucl_num, asymp_jasa) & + result(info) + use qmckl + implicit none + integer(qmckl_context), intent(in) :: context + integer*8 , intent(in) :: aord_num + double precision , intent(in) :: a_vector(aord_num + 1) + integer*8 , intent(in) :: nucl_num + double precision , intent(in) :: rescale_factor_en(nucl_num) + double precision , intent(out) :: asymp_jasa(nucl_num) + + integer*8 :: i, j, p + double precision :: kappa_inv, x, asym_one + + + info = QMCKL_SUCCESS + + if (context == QMCKL_NULL_CONTEXT) then + info = QMCKL_INVALID_CONTEXT + return + endif + + if (aord_num < 0) then + info = QMCKL_INVALID_ARG_2 + return + endif + + do i=1,nucl_num + + kappa_inv = 1.0d0 / rescale_factor_en(i) + + asymp_jasa(i) = a_vector(1) * kappa_inv / (1.0d0 + a_vector(2) * kappa_inv) + + x = kappa_inv + do p = 1, aord_num + x = x * kappa_inv + asymp_jasa(i) = asymp_jasa(i) + a_vector(p + 1) * x + end do + + end do + +end function qmckl_compute_jastrow_asymp_jasa_f + #+end_src + +#+begin_src c :comments org :tangle (eval c) :noweb yes +qmckl_exit_code qmckl_compute_jastrow_asymp_jasa ( + const qmckl_context context, + const int64_t aord_num, + const double* a_vector, + double* const rescale_factor_en, + const int64_t nucl_num, + double* const asymp_jasa ) { + + if (context == QMCKL_NULL_CONTEXT){ + return QMCKL_INVALID_CONTEXT; + } + + if (aord_num < 0) { + return QMCKL_INVALID_ARG_2; + } + + for (int i = 0 ; i <= nucl_num; ++i) { + const double kappa_inv = 1.0 / rescale_factor_en[i]; + asymp_jasa[i] = a_vector[0] * kappa_inv / (1.0 + a_vector[1] * kappa_inv); + + double x = kappa_inv; + for (int p = 1; p < aord_num; ++p){ + x *= kappa_inv; + asymp_jasa[i] = asymp_jasa[i] + a_vector[p + 1] * x; + } + } + + return QMCKL_SUCCESS; +} +#+end_src + +# #+CALL: generate_c_header(table=qmckl_asymp_jasa_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none + qmckl_exit_code qmckl_compute_jastrow_asymp_jasa ( + const qmckl_context context, + const int64_t aord_num, + const double* a_vector, + double* const rescale_factor_en, + const int64_t nucl_num, + double* const asymp_jasa ); + #+end_src + +*** Test + #+name: asymp_jasa + #+begin_src python :results output :exports none :noweb yes +import numpy as np + +<> + + +asymp_jasa = a_vector[0] * kappa_inv / (1.0 + a_vector[1]*kappa_inv) +x = kappa_inv +for p in range(1,aord_num): + x = x * kappa_inv + asymp_jasa += a_vector[p + 1] * x +print("asymp_jasa[i] : ", asymp_jasa) + + #+end_src + + #+RESULTS: asymp_jasa + : asymp_jasa[i] : [-0.548554] + + #+begin_src c :tangle (eval c_test) +double asymp_jasa[2]; +rc = qmckl_get_jastrow_asymp_jasa(context, asymp_jasa, nucl_num); + +// calculate asymp_jasb +printf("%e %e\n", asymp_jasa[0], -0.548554); +printf("%e %e\n", asymp_jasa[1], -0.548554); +assert(fabs(-0.548554 - asymp_jasa[0]) < 1.e-12); +assert(fabs(-0.548554 - asymp_jasa[1]) < 1.e-12); + + #+end_src + ** Electron-nucleus component \(f_{en}\) Calculate the electron-electron jastrow component ~factor_en~ using the ~a_vector~