From 7c226d0a99fc33f425dab73e4e1f7872f0bfd814 Mon Sep 17 00:00:00 2001 From: vijay gopal chilkuri Date: Tue, 6 Jul 2021 15:51:51 +0530 Subject: [PATCH] Finished factor_en. #22 --- org/qmckl_distance.org | 2 +- org/qmckl_jastrow.org | 482 +++++++++++++++++++++++++++++++++++++---- tests/n2.h | 6 +- 3 files changed, 444 insertions(+), 46 deletions(-) diff --git a/org/qmckl_distance.org b/org/qmckl_distance.org index 082209b..1465012 100644 --- a/org/qmckl_distance.org +++ b/org/qmckl_distance.org @@ -923,7 +923,7 @@ integer function qmckl_distance_rescaled_f(context, transa, transb, m, n, & if (transb == 'N' .or. transb == 'n') then continue - else if (transa == 'T' .or. transa == 't') then + else if (transb == 'T' .or. transb == 't') then transab = transab + 2 else transab = -100 diff --git a/org/qmckl_jastrow.org b/org/qmckl_jastrow.org index c731d8e..461a445 100644 --- a/org/qmckl_jastrow.org +++ b/org/qmckl_jastrow.org @@ -72,27 +72,28 @@ int main() { The following data stored in the context: #+NAME: qmckl_jastrow_args - |------------+--------------------------------------------+-----+-------------------------------------------------------------------| - | ~int32_t~ | ~uninitialized~ | in | Keeps bit set for uninitialized data | - | ~int64_t~ | ~aord_num~ | in | The number of a coeffecients | - | ~int64_t~ | ~bord_num~ | in | The number of b coeffecients | - | ~int64_t~ | ~cord_num~ | in | The number of c coeffecients | - | ~uint64_t~ | ~type_nucl_num~ | in | Number of Nucleii types | - | ~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[walk_num]~ | out | Jastrow factor: electron-electron part | - | ~double~ | ~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 | - | ~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][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 | - | ~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 | - | ~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 | + |-----------+--------------------------------------------+-----+-------------------------------------------------------------------| + | ~int32_t~ | ~uninitialized~ | in | Keeps bit set for uninitialized data | + | ~int64_t~ | ~aord_num~ | in | The number of a coeffecients | + | ~int64_t~ | ~bord_num~ | in | The number of b coeffecients | + | ~int64_t~ | ~cord_num~ | in | The number of c coeffecients | + | ~int64_t~ | ~type_nucl_num~ | in | Number of Nucleii types | + | ~int64_t~ | ~type_nucl_vector[nucl_num]~ | in | IDs of types of Nucleii | + | ~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[walk_num]~ | out | Jastrow factor: electron-electron part | + | ~double~ | ~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 | + | ~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][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 | + | ~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 | + | ~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 | computed data: @@ -172,6 +173,16 @@ ee_distance_rescaled = [ 0.956034127544344 ,0.931221472309472 ,0.540903688625053 , 0.000000000000000E+000]] +en_distance_rescaled = np.transpose(np.array([ +[ 0.443570948411811 , 0.467602196999105 , 0.893870160799932 , + 0.864347190364447 , 0.976608182392358 , 0.187563183468210 , + 0.426404699872689 , 0.665107090128166 , 0.885246991424583 , + 0.924902909715270 ], +[ 0.899360150637444 , 0.860035135365386 , 0.979659405613798 , + 6.140678415933776E-002, 0.835118398056681 , 0.884071658981068 , + 0.923860000907362 , 0.905203414522289 , 0.211286300932359 , + 0.492104840907350 ]])) + # symmetrize it for i in range(elec_num): for j in range(elec_num): @@ -182,8 +193,15 @@ aord_num = 5 bord_num = 5 cord_num = 23 dim_cord_vec = 23 -aord_vector = [ 0.000000000000000E+000, 0.000000000000000E+000, -0.380512000000000E+000, - -0.157996000000000E+000, -3.155800000000000E-002, 2.151200000000000E-002] +type_nucl_vector = [ 1, 1] +aord_vector = [ +[0.000000000000000E+000], +[0.000000000000000E+000], +[-0.380512000000000E+000], +[-0.157996000000000E+000], +[-3.155800000000000E-002], +[2.151200000000000E-002]] + bord_vector = [ 0.500000000000000E-000, 0.153660000000000E-000, 6.722620000000000E-002, 2.157000000000000E-002, 7.309600000000000E-003, 2.866000000000000E-003] cord_vector = [ 0.571702000000000E-000, -0.514253000000000E-000, -0.513043000000000E-000, @@ -260,6 +278,7 @@ typedef struct qmckl_jastrow_struct{ int64_t factor_ee_deriv_e_date; int64_t factor_en_deriv_e_date; int64_t factor_een_deriv_e_date; + int64_t* type_nucl_vector; double * aord_vector; double * bord_vector; double * cord_vector; @@ -301,7 +320,7 @@ qmckl_exit_code qmckl_init_jastrow(qmckl_context context) { qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - ctx->jastrow.uninitialized = (1 << 5) - 1; + ctx->jastrow.uninitialized = (1 << 6) - 1; /* Default values */ @@ -312,13 +331,14 @@ qmckl_exit_code qmckl_init_jastrow(qmckl_context context) { ** Access functions #+begin_src c :comments org :tangle (eval h_func) :exports none -qmckl_exit_code qmckl_get_jastrow_aord_num (qmckl_context context, int64_t* const aord_num); -qmckl_exit_code qmckl_get_jastrow_bord_num (qmckl_context context, int64_t* const bord_num); -qmckl_exit_code qmckl_get_jastrow_cord_num (qmckl_context context, int64_t* const bord_num); -qmckl_exit_code qmckl_get_jastrow_type_nucl_num (qmckl_context context, int64_t* const type_nucl_num); -qmckl_exit_code qmckl_get_jastrow_aord_vector (qmckl_context context, double * const aord_vector); -qmckl_exit_code qmckl_get_jastrow_bord_vector (qmckl_context context, double * const bord_vector); -qmckl_exit_code qmckl_get_jastrow_cord_vector (qmckl_context context, double * const cord_vector); +qmckl_exit_code qmckl_get_jastrow_aord_num (qmckl_context context, int64_t* const aord_num); +qmckl_exit_code qmckl_get_jastrow_bord_num (qmckl_context context, int64_t* const bord_num); +qmckl_exit_code qmckl_get_jastrow_cord_num (qmckl_context context, int64_t* const bord_num); +qmckl_exit_code qmckl_get_jastrow_type_nucl_num (qmckl_context context, int64_t* const type_nucl_num); +qmckl_exit_code qmckl_get_jastrow_type_nucl_vector (qmckl_context context, int64_t* const type_nucl_num); +qmckl_exit_code qmckl_get_jastrow_aord_vector (qmckl_context context, double * const aord_vector); +qmckl_exit_code qmckl_get_jastrow_bord_vector (qmckl_context context, double * const bord_vector); +qmckl_exit_code qmckl_get_jastrow_cord_vector (qmckl_context context, double * const cord_vector); #+end_src Along with these core functions, calculation of the jastrow factor @@ -462,6 +482,33 @@ qmckl_exit_code qmckl_get_jastrow_type_nucl_num (const qmckl_context context, in return QMCKL_SUCCESS; } +qmckl_exit_code qmckl_get_jastrow_type_nucl_vector (const qmckl_context context, int64_t * const type_nucl_vector) { + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return (char) 0; + } + + if (type_nucl_vector == NULL) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_get_jastrow_type_nucl_vector", + "type_nucl_vector is a null pointer"); + } + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + int32_t mask = 1 << 2; + + if ( (ctx->jastrow.uninitialized & mask) != 0) { + return QMCKL_NOT_PROVIDED; + } + + assert (ctx->jastrow.type_nucl_vector != NULL); + memcpy(type_nucl_vector, ctx->jastrow.type_nucl_vector, ctx->jastrow.type_nucl_num*sizeof(int64_t)); + return QMCKL_SUCCESS; +} + qmckl_exit_code qmckl_get_jastrow_aord_vector (const qmckl_context context, double * const aord_vector) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { @@ -478,7 +525,7 @@ qmckl_exit_code qmckl_get_jastrow_aord_vector (const qmckl_context context, doub qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - int32_t mask = 1 << 2; + int32_t mask = 1 << 3; if ( (ctx->jastrow.uninitialized & mask) != 0) { return QMCKL_NOT_PROVIDED; @@ -505,7 +552,7 @@ qmckl_exit_code qmckl_get_jastrow_bord_vector (const qmckl_context context, doub qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - int32_t mask = 1 << 3; + int32_t mask = 1 << 4; if ( (ctx->jastrow.uninitialized & mask) != 0) { return QMCKL_NOT_PROVIDED; @@ -532,7 +579,7 @@ qmckl_exit_code qmckl_get_jastrow_cord_vector (const qmckl_context context, doub qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - int32_t mask = 1 << 4; + int32_t mask = 1 << 5; if ( (ctx->jastrow.uninitialized & mask) != 0) { return QMCKL_NOT_PROVIDED; @@ -551,12 +598,13 @@ qmckl_exit_code qmckl_get_jastrow_cord_vector (const qmckl_context context, doub called. #+begin_src c :comments org :tangle (eval h_func) -qmckl_exit_code qmckl_set_jastrow_ord_num (qmckl_context context, const int64_t aord_num, const int64_t bord_num, const int64_t cord_num); -qmckl_exit_code qmckl_set_jastrow_type_nucl_num (qmckl_context context, const int64_t type_nucl_num); -qmckl_exit_code qmckl_set_jastrow_aord_vector (qmckl_context context, const double * aord_vector); -qmckl_exit_code qmckl_set_jastrow_bord_vector (qmckl_context context, const double * bord_vector); -qmckl_exit_code qmckl_set_jastrow_cord_vector (qmckl_context context, const double * cord_vector); -qmckl_exit_code qmckl_set_jastrow_dependencies (qmckl_context context); +qmckl_exit_code qmckl_set_jastrow_ord_num (qmckl_context context, const int64_t aord_num, const int64_t bord_num, const int64_t cord_num); +qmckl_exit_code qmckl_set_jastrow_type_nucl_num (qmckl_context context, const int64_t type_nucl_num); +qmckl_exit_code qmckl_set_jastrow_type_nucl_vector (qmckl_context context, const int64_t* type_nucl_vector, const int64_t nucl_num); +qmckl_exit_code qmckl_set_jastrow_aord_vector (qmckl_context context, const double * aord_vector); +qmckl_exit_code qmckl_set_jastrow_bord_vector (qmckl_context context, const double * bord_vector); +qmckl_exit_code qmckl_set_jastrow_cord_vector (qmckl_context context, const double * cord_vector); +qmckl_exit_code qmckl_set_jastrow_dependencies (qmckl_context context); #+end_src #+NAME:pre2 @@ -630,11 +678,61 @@ qmckl_exit_code qmckl_set_jastrow_type_nucl_num(qmckl_context context, const int <> } -qmckl_exit_code qmckl_set_jastrow_aord_vector(qmckl_context context, double const * aord_vector) { +qmckl_exit_code qmckl_set_jastrow_type_nucl_vector(qmckl_context context, int64_t const * type_nucl_vector, const int64_t nucl_num) { <> int32_t mask = 1 << 2; + int64_t type_nucl_num; + qmckl_exit_code rc = qmckl_get_jastrow_type_nucl_num(context, &type_nucl_num); + if (rc != QMCKL_SUCCESS) return rc; + + if (type_nucl_num == 0) { + return qmckl_failwith( context, + QMCKL_FAILURE, + "qmckl_set_jastrow_type_nucl_vector", + "type_nucl_num is not set"); + } + + if (type_nucl_vector == NULL) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_set_jastrow_type_nucl_vector", + "type_nucl_vector = NULL"); + } + + if (ctx->jastrow.type_nucl_vector != NULL) { + qmckl_exit_code rc = qmckl_free(context, ctx->jastrow.type_nucl_vector); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_set_type_nucl_vector", + NULL); + } + } + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = nucl_num * sizeof(int64_t); + int64_t* new_array = (int64_t*) qmckl_malloc(context, mem_info); + + if(new_array == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_set_jastrow_type_nucl_vector", + NULL); + } + + memcpy(new_array, type_nucl_vector, mem_info.size); + + ctx->jastrow.type_nucl_vector = new_array; + + <> +} + +qmckl_exit_code qmckl_set_jastrow_aord_vector(qmckl_context context, double const * aord_vector) { +<> + + int32_t mask = 1 << 3; + int64_t aord_num; qmckl_exit_code rc = qmckl_get_jastrow_aord_num(context, &aord_num); if (rc != QMCKL_SUCCESS) return rc; @@ -687,7 +785,7 @@ qmckl_exit_code qmckl_set_jastrow_aord_vector(qmckl_context context, double cons qmckl_exit_code qmckl_set_jastrow_bord_vector(qmckl_context context, double const * bord_vector) { <> - int32_t mask = 1 << 3; + int32_t mask = 1 << 4; int64_t bord_num; qmckl_exit_code rc = qmckl_get_jastrow_bord_num(context, &bord_num); @@ -737,7 +835,7 @@ qmckl_exit_code qmckl_set_jastrow_bord_vector(qmckl_context context, double cons qmckl_exit_code qmckl_set_jastrow_cord_vector(qmckl_context context, double const * cord_vector) { <> - int32_t mask = 1 << 4; + int32_t mask = 1 << 5; int64_t cord_num; qmckl_exit_code rc = qmckl_get_jastrow_cord_num(context, &cord_num); @@ -807,7 +905,7 @@ qmckl_exit_code qmckl_set_jastrow_dependencies(qmckl_context context) { NULL); } - int32_t mask = 1 << 5; + int32_t mask = 1 << 6; <> } @@ -1272,6 +1370,7 @@ print("asymp_jasb[1] : ", asymp_jasb[1]) #+begin_src c :tangle (eval c_test) assert(qmckl_electron_provided(context)); +int64_t* type_nucl_vector = &(n2_type_nucl_vector[0]); double* aord_vector = &(n2_aord_vector[0][0]); double* bord_vector = &(n2_bord_vector[0]); double* cord_vector = &(n2_cord_vector[0][0]); @@ -1285,6 +1384,8 @@ rc = qmckl_set_jastrow_ord_num(context, n2_aord_num, n2_bord_num, n2_cord_num); assert(rc == QMCKL_SUCCESS); rc = qmckl_set_jastrow_type_nucl_num(context, n2_type_nucl_num); assert(rc == QMCKL_SUCCESS); +rc = qmckl_set_jastrow_type_nucl_vector(context, n2_type_nucl_vector, nucl_num); +assert(rc == QMCKL_SUCCESS); rc = qmckl_set_jastrow_aord_vector(context, aord_vector); assert(rc == QMCKL_SUCCESS); rc = qmckl_set_jastrow_bord_vector(context, bord_vector); @@ -2029,6 +2130,299 @@ assert(fabs(factor_ee_deriv_e[0][3][0]-1.5111672803213185 ) < 1.e-12); #+end_src +** Electron-nucleus component \(f_{en}\) + + Calculate the electron-electron jastrow component ~factor_en~ using the ~aord_vector~ + coeffecients and the electron-nucleus rescaled distances ~en_distance_rescaled~. + + \[ +f_{en} = \sum_{i,jjastrow.factor_en, 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_en(qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_provide_factor_en(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 en rescaled distance is provided */ + rc = qmckl_provide_en_distance_rescaled(context); + if(rc != QMCKL_SUCCESS) return rc; + + /* Compute if necessary */ + if (ctx->date > ctx->jastrow.factor_en_date) { + + /* Allocate array */ + if (ctx->jastrow.factor_en == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->electron.walk_num * sizeof(double); + double* factor_en = (double*) qmckl_malloc(context, mem_info); + + if (factor_en == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_provide_factor_en", + NULL); + } + ctx->jastrow.factor_en = factor_en; + } + + qmckl_exit_code rc = + qmckl_compute_factor_en(context, + ctx->electron.walk_num, + ctx->electron.num, + ctx->nucleus.num, + ctx->jastrow.type_nucl_num, + ctx->jastrow.type_nucl_vector, + ctx->jastrow.aord_num, + ctx->jastrow.aord_vector, + ctx->electron.en_distance_rescaled, + ctx->jastrow.factor_en); + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->jastrow.factor_en_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + #+end_src + +*** Compute + :PROPERTIES: + :Name: qmckl_compute_factor_en + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_factor_en_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 | nucl_num | in | Number of nucleii | + | int64_t | type_nucl_num | in | Number of unique nuclei | + | int64_t | type_nucl_vector[type_nucl_num] | in | IDs of unique nucleii | + | int64_t | aord_num | in | Number of coefficients | + | double | aord_vector[aord_num + 1][type_nucl_num] | in | List of coefficients | + | double | en_distance_rescaled[walk_num][nucl_num][elec_num] | in | Electron-nucleus distances | + | double | factor_en[walk_num] | out | Electron-nucleus jastrow | + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer function qmckl_compute_factor_en_f(context, walk_num, elec_num, nucl_num, type_nucl_num, & + type_nucl_vector, aord_num, aord_vector, & + en_distance_rescaled, factor_en) & + result(info) + use qmckl + implicit none + integer(qmckl_context), intent(in) :: context + integer*8 , intent(in) :: walk_num, elec_num, aord_num, nucl_num, type_nucl_num + integer*8 , intent(in) :: type_nucl_vector(type_nucl_num) + double precision , intent(in) :: aord_vector(aord_num, nucl_num) + double precision , intent(in) :: en_distance_rescaled(walk_num, elec_num, nucl_num) + double precision , intent(out) :: factor_en(walk_num) + + integer*8 :: i, a, p, ipar, nw + double precision :: 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 (nucl_num <= 0) then + info = QMCKL_INVALID_ARG_4 + return + endif + + if (aord_num <= 0) then + info = QMCKL_INVALID_ARG_7 + return + endif + + factor_en = 0.0d0 + + do nw =1, walk_num + do a = 1, nucl_num + do i = 1, elec_num + x = en_distance_rescaled(nw, i, a) + power_ser = 0.0d0 + + do p = 2, aord_num + x = x * en_distance_rescaled(nw, i, a) + power_ser = power_ser + aord_vector(p + 1, type_nucl_vector(a)) * x + end do + + factor_en(nw) = factor_en(nw) + aord_vector(1, type_nucl_vector(a)) * & + en_distance_rescaled(nw, i, a) / & + (1.0d0 + aord_vector(2, type_nucl_vector(a)) * & + en_distance_rescaled(nw, i, a)) & + + power_ser + + end do + end do + end do + +end function qmckl_compute_factor_en_f + #+end_src + + #+CALL: generate_c_header(table=qmckl_factor_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_factor_en ( + const qmckl_context context, + 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* aord_vector, + const double* en_distance_rescaled, + double* const factor_en ); + #+end_src + + + #+CALL: generate_c_interface(table=qmckl_factor_en_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_en & + (context, & + walk_num, & + elec_num, & + nucl_num, & + type_nucl_num, & + type_nucl_vector, & + aord_num, & + aord_vector, & + en_distance_rescaled, & + factor_en) & + 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 :: nucl_num + integer (c_int64_t) , intent(in) , value :: type_nucl_num + integer (c_int64_t) , intent(in) :: type_nucl_vector(type_nucl_num) + integer (c_int64_t) , intent(in) , value :: aord_num + real (c_double ) , intent(in) :: aord_vector(type_nucl_num,aord_num + 1) + real (c_double ) , intent(in) :: en_distance_rescaled(elec_num,nucl_num,walk_num) + real (c_double ) , intent(out) :: factor_en(walk_num) + + integer(c_int32_t), external :: qmckl_compute_factor_en_f + info = qmckl_compute_factor_en_f & + (context, & + walk_num, & + elec_num, & + nucl_num, & + type_nucl_num, & + type_nucl_vector, & + aord_num, & + aord_vector, & + en_distance_rescaled, & + factor_en) + + end function qmckl_compute_factor_en + #+end_src + +*** Test + #+begin_src python :results output :exports none :noweb yes +import numpy as np + +<> + +factor_en = 0.0 +for a in range(0,nucl_num): + for i in range(0,elec_num): + x = en_distance_rescaled[i][a] + pow_ser = 0.0 + + for p in range(2,aord_num+1): + x = x * en_distance_rescaled[i][a] + pow_ser = pow_ser + aord_vector[(p-1) + 1][type_nucl_vector[a]-1] * x + + factor_en = factor_en + aord_vector[0][type_nucl_vector[a]-1] * en_distance_rescaled[i][a] \ + / (1.0 + aord_vector[1][type_nucl_vector[a]-1] * en_distance_rescaled[i][a]) \ + + pow_ser +print("factor_en :",factor_en) + + #+end_src + + #+RESULTS: + : factor_en : -5.865822569188727 + + + #+begin_src c :tangle (eval c_test) +/* Check if Jastrow is properly initialized */ +assert(qmckl_jastrow_provided(context)); + +double factor_en[walk_num]; +rc = qmckl_get_jastrow_factor_en(context, factor_en); + +// calculate factor_en +assert(fabs(factor_en[0]+5.865822569188727) < 1.e-12); + + #+end_src + + * End of files :noexport: #+begin_src c :tangle (eval h_private_type) diff --git a/tests/n2.h b/tests/n2.h index 80f887a..ad675c0 100644 --- a/tests/n2.h +++ b/tests/n2.h @@ -5,7 +5,7 @@ double n2_charge[n2_nucl_num] = { 5., 5.}; double n2_nucl_coord[3][n2_nucl_num] = { {0.000000, 0.000000 }, {0.000000, 0.000000 }, - {2.059801, 2.059801 } }; + {0.000000, 2.059801 } }; #define n2_elec_up_num ((int64_t) 5) #define n2_elec_dn_num ((int64_t) 5) @@ -32,6 +32,10 @@ double n2_elec_coord[n2_walk_num][n2_elec_num][3] = { { #define n2_cord_num ((int64_t) 23) #define n2_dim_cord_vec ((int64_t) 23) +int64_t n2_type_nucl_vector[n2_nucl_num] = { + 1, + 1}; + double n2_aord_vector[n2_aord_num + 1][n2_type_nucl_num] = { { 0. }, { 0. },