From 348845511077cd5e4b9044b0a4f41ea55ce49063 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 31 Jan 2022 16:47:28 +0100 Subject: [PATCH] Work on Jastrow --- org/qmckl_jastrow.org | 329 ++++++++++----------- share/qmckl/test_data/chbrclf/metadata.txt | 1 + 2 files changed, 150 insertions(+), 180 deletions(-) diff --git a/org/qmckl_jastrow.org b/org/qmckl_jastrow.org index ad6826f..0d92e00 100644 --- a/org/qmckl_jastrow.org +++ b/org/qmckl_jastrow.org @@ -67,7 +67,7 @@ int main() { #include "qmckl_jastrow_private_func.h" #include "qmckl_jastrow_private_type.h" #+end_src - + * Context :PROPERTIES: :Name: qmckl_jastrow @@ -116,7 +116,7 @@ int main() { | ~lkpm_combined_index_date~ | ~uint64_t~ | Transform l,k,p, and m into consecutive indices | | | ~tmp_c~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]~ | out | vector of non-zero coefficients | | ~dtmp_c~ | ~double[walk_num][elec_num][4][nucl_num][0:cord_num][0:cord_num-1]~ | vector of non-zero coefficients | | - + | ~een_rescaled_n~ | ~double[walk_num][elec_num][nucl_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord | | | ~een_rescaled_n_date~ | ~uint64_t~ | Keep track of the date of creation | | | ~een_rescaled_e_deriv_e~ | ~double[walk_num][elec_num][4][elec_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons | | @@ -139,7 +139,7 @@ nucl_coord = np.array([ [0.000000, 0.000000 ], [0.000000, 2.059801 ] ]) elec_coord = [[[-0.250655104764153 , 0.503070975550133 , -0.166554344502303], - [-0.587812193472177 , -0.128751981129274 , 0.187773606533075], + [-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 ], @@ -192,13 +192,13 @@ ee_distance_rescaled = [ 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.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.899360150637444 , 0.860035135365386 , 0.979659405613798 , + 6.140678415933776E-002, 0.835118398056681 , 0.884071658981068 , + 0.923860000907362 , 0.905203414522289 , 0.211286300932359 , 0.492104840907350 ]])) # symmetrize it @@ -212,7 +212,7 @@ bord_num = 5 cord_num = 5 dim_cord_vect= 23 type_nucl_vector = [ 1, 1] -aord_vector = [ +aord_vector = [ [0.000000000000000E+000], [0.000000000000000E+000], [-0.380512000000000E+000], @@ -222,7 +222,7 @@ aord_vector = [ 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, +cord_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, -5.270400000000000E-003, -7.499999999999999E-005, -8.301649999999999E-002, @@ -231,7 +231,7 @@ cord_vector = [ 0.571702000000000E-000, -0.514253000000000E-000, -0.513043000000 2.614000000000000E-003, -1.477000000000000E-003, -1.137000000000000E-003, -4.010475000000000E-002, 6.106710000000000E-003 ] cord_vector_full = [ -[ 0.571702000000000E-000, -0.514253000000000E-000, -0.513043000000000E-000, +[ 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, -5.270400000000000E-003, -7.499999999999999E-005, -8.301649999999999E-002, @@ -239,7 +239,7 @@ cord_vector_full = [ -4.099100000000000E-003, 4.327600000000000E-003, -1.654470000000000E-003, 2.614000000000000E-003, -1.477000000000000E-003, -1.137000000000000E-003, -4.010475000000000E-002, 6.106710000000000E-003 ], -[ 0.571702000000000E-000, -0.514253000000000E-000, -0.513043000000000E-000, +[ 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, -5.270400000000000E-003, -7.499999999999999E-005, -8.301649999999999E-002, @@ -248,7 +248,7 @@ cord_vector_full = [ 2.614000000000000E-003, -1.477000000000000E-003, -1.137000000000000E-003, -4.010475000000000E-002, 6.106710000000000E-003 ], ] -lkpm_combined_index = [[1 , 1 , 2 , 0], +lkpm_combined_index = [[1 , 1 , 2 , 0], [0 , 0 , 2 , 1], [1 , 2 , 3 , 0], [2 , 1 , 3 , 0], @@ -275,7 +275,7 @@ lkpm_combined_index = [[1 , 1 , 2 , 0], kappa = 1.0 kappa_inv = 1.0/kappa #+END_SRC - + #+RESULTS: jastrow_data ** Data structure @@ -336,10 +336,10 @@ typedef struct qmckl_jastrow_struct{ Some values are initialized by default, and are not concerned by this mechanism. - #+begin_src c :comments org :tangle (eval h_func) + #+begin_src c :comments org :tangle (eval h_func) qmckl_exit_code qmckl_init_jastrow(qmckl_context context); #+end_src - + #+begin_src c :comments org :tangle (eval c) qmckl_exit_code qmckl_init_jastrow(qmckl_context context) { @@ -350,14 +350,14 @@ 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 << 6) - 1; + ctx->jastrow.uninitialized = (1 << 5) - 1; /* Default values */ return QMCKL_SUCCESS; } #+end_src - + ** Access functions #+begin_src c :comments org :tangle (eval h_func) :exports none @@ -409,7 +409,7 @@ qmckl_exit_code qmckl_get_jastrow_aord_num (const qmckl_context context, int64_t if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return (char) 0; } - + if (aord_num == NULL) { return qmckl_failwith( context, QMCKL_INVALID_ARG_2, @@ -436,7 +436,7 @@ qmckl_exit_code qmckl_get_jastrow_bord_num (const qmckl_context context, int64_t if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return (char) 0; } - + if (bord_num == NULL) { return qmckl_failwith( context, QMCKL_INVALID_ARG_2, @@ -463,7 +463,7 @@ qmckl_exit_code qmckl_get_jastrow_cord_num (const qmckl_context context, int64_t if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return (char) 0; } - + if (cord_num == NULL) { return qmckl_failwith( context, QMCKL_INVALID_ARG_2, @@ -490,7 +490,7 @@ qmckl_exit_code qmckl_get_jastrow_type_nucl_num (const qmckl_context context, in if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return (char) 0; } - + if (type_nucl_num == NULL) { return qmckl_failwith( context, QMCKL_INVALID_ARG_2, @@ -524,7 +524,7 @@ qmckl_exit_code qmckl_get_jastrow_type_nucl_vector (const qmckl_context context, "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); @@ -551,7 +551,7 @@ qmckl_exit_code qmckl_get_jastrow_aord_vector (const qmckl_context context, doub "qmckl_get_jastrow_aord_vector", "aord_vector is a null pointer"); } - + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); @@ -578,7 +578,7 @@ qmckl_exit_code qmckl_get_jastrow_bord_vector (const qmckl_context context, doub "qmckl_get_jastrow_bord_vector", "bord_vector is a null pointer"); } - + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); @@ -605,7 +605,7 @@ qmckl_exit_code qmckl_get_jastrow_cord_vector (const qmckl_context context, doub "qmckl_get_jastrow_cord_vector", "cord_vector is a null pointer"); } - + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); @@ -634,7 +634,6 @@ qmckl_exit_code qmckl_set_jastrow_type_nucl_vector (qmckl_context context, con 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 @@ -651,8 +650,8 @@ qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; ctx->jastrow.uninitialized &= ~mask; ctx->jastrow.provided = (ctx->jastrow.uninitialized == 0); if (ctx->jastrow.provided) { - //qmckl_exit_code rc_ = qmckl_set_jastrow_dependencies(context); - //if (rc_ != QMCKL_SUCCESS) return rc_; + qmckl_exit_code rc_ = qmckl_finalize_jastrow(context); + if (rc_ != QMCKL_SUCCESS) return rc_; } return QMCKL_SUCCESS; @@ -716,7 +715,7 @@ qmckl_exit_code qmckl_set_jastrow_type_nucl_vector(qmckl_context context, int64_ 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, @@ -730,7 +729,7 @@ qmckl_exit_code qmckl_set_jastrow_type_nucl_vector(qmckl_context context, int64_ "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) { @@ -766,11 +765,11 @@ qmckl_exit_code qmckl_set_jastrow_aord_vector(qmckl_context context, double cons int64_t aord_num; qmckl_exit_code rc = qmckl_get_jastrow_aord_num(context, &aord_num); if (rc != QMCKL_SUCCESS) return rc; - + int64_t type_nucl_num; rc = qmckl_get_jastrow_type_nucl_num(context, &type_nucl_num); if (rc != QMCKL_SUCCESS) return rc; - + if (aord_num == 0) { return qmckl_failwith( context, QMCKL_FAILURE, @@ -784,7 +783,7 @@ qmckl_exit_code qmckl_set_jastrow_aord_vector(qmckl_context context, double cons "qmckl_set_jastrow_aord_vector", "aord_vector = NULL"); } - + if (ctx->jastrow.aord_vector != NULL) { qmckl_exit_code rc = qmckl_free(context, ctx->jastrow.aord_vector); if (rc != QMCKL_SUCCESS) { @@ -834,7 +833,7 @@ qmckl_exit_code qmckl_set_jastrow_bord_vector(qmckl_context context, double cons "qmckl_set_jastrow_bord_vector", "bord_vector = NULL"); } - + if (ctx->jastrow.bord_vector != NULL) { qmckl_exit_code rc = qmckl_free(context, ctx->jastrow.bord_vector); if (rc != QMCKL_SUCCESS) { @@ -873,7 +872,7 @@ qmckl_exit_code qmckl_set_jastrow_cord_vector(qmckl_context context, double cons int64_t dim_cord_vect; rc = qmckl_get_jastrow_dim_cord_vect(context, &dim_cord_vect); if (rc != QMCKL_SUCCESS) return rc; - + int64_t type_nucl_num; rc = qmckl_get_jastrow_type_nucl_num(context, &type_nucl_num); if (rc != QMCKL_SUCCESS) return rc; @@ -891,7 +890,7 @@ qmckl_exit_code qmckl_set_jastrow_cord_vector(qmckl_context context, double cons "qmckl_set_jastrow_cord_vector", "cord_vector = NULL"); } - + if (ctx->jastrow.cord_vector != NULL) { qmckl_exit_code rc = qmckl_free(context, ctx->jastrow.cord_vector); if (rc != QMCKL_SUCCESS) { @@ -919,30 +918,6 @@ qmckl_exit_code qmckl_set_jastrow_cord_vector(qmckl_context context, double cons <> } -qmckl_exit_code qmckl_set_jastrow_dependencies(qmckl_context context) { -<> - - /* Check for electron data */ - if (!(ctx->electron.provided)) { - return qmckl_failwith( context, - QMCKL_NOT_PROVIDED, - "qmckl_provide_ee_distance", - NULL); - } - - /* Check for nucleus data */ - if (!(ctx->nucleus.provided)) { - return qmckl_failwith( context, - QMCKL_NOT_PROVIDED, - "qmckl_provide_en_distance", - NULL); - } - - int32_t mask = 1 << 6; - - <> -} - #+end_src When the required information is completely entered, other data structures are @@ -989,17 +964,13 @@ qmckl_exit_code qmckl_finalize_jastrow(qmckl_context context) { NULL); } - qmckl_exit_code rc = QMCKL_FAILURE; + qmckl_exit_code rc = QMCKL_SUCCESS; return rc; - /* ----------------------------------- */ - /* Start calculation of data */ - /* ----------------------------------- */ - } #+end_src - + ** Test #+begin_src c :tangle (eval c_test) /* Reference input data */ @@ -1169,7 +1140,7 @@ for (int64_t i=0 ; ijastrow.factor_ee = factor_ee; } - + qmckl_exit_code rc = qmckl_compute_factor_ee(context, ctx->electron.walk_num, @@ -1599,7 +1568,7 @@ integer function qmckl_compute_factor_ee_f(context, walk_num, elec_num, up_num, endif factor_ee = 0.0d0 - + do nw =1, walk_num do j = 1, elec_num do i = 1, j - 1 @@ -1617,7 +1586,7 @@ integer function qmckl_compute_factor_ee_f(context, walk_num, elec_num, up_num, 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) * & @@ -1644,7 +1613,7 @@ end function qmckl_compute_factor_ee_f const double* bord_vector, const double* ee_distance_rescaled, const double* asymp_jasb, - double* const factor_ee ); + double* const factor_ee ); #+end_src @@ -1707,7 +1676,7 @@ for i in range(0,elec_num): 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 @@ -1715,7 +1684,7 @@ for i in range(0,elec_num): 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 @@ -1728,7 +1697,7 @@ print("factor_ee :",factor_ee) : 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 */ @@ -1741,17 +1710,17 @@ rc = qmckl_get_jastrow_factor_ee(context, factor_ee); assert(fabs(factor_ee[0]+4.282760865958113) < 1.e-12); #+end_src - + ** Electron-electron component derivative \(f'_{ee}\) Calculate the derivative of the ~factor_ee~ using the ~ee_distance_rescaled~ and the electron-electron rescaled distances derivatives ~ee_distance_rescaled_deriv_e~. - There are four components, the gradient which has 3 components in the \(x, y, z\) + There are four components, the gradient which has 3 components in the \(x, y, z\) directions and the laplacian as the last component. TODO: Add equation - + *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_jastrow_factor_ee_deriv_e(qmckl_context context, double* const factor_ee_deriv_e); @@ -1783,7 +1752,7 @@ qmckl_exit_code qmckl_get_jastrow_factor_ee_deriv_e(qmckl_context context, doubl #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_provide_factor_ee_deriv_e(qmckl_context context); #+end_src - + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_provide_factor_ee_deriv_e(qmckl_context context) { @@ -1823,7 +1792,7 @@ qmckl_exit_code qmckl_provide_factor_ee_deriv_e(qmckl_context context) } ctx->jastrow.factor_ee_deriv_e = factor_ee_deriv_e; } - + qmckl_exit_code rc = qmckl_compute_factor_ee_deriv_e(context, ctx->electron.walk_num, @@ -1913,7 +1882,7 @@ integer function qmckl_compute_factor_ee_deriv_e_f(context, walk_num, elec_num, factor_ee_deriv_e = 0.0d0 third = 1.0d0 / 3.0d0 - + do nw =1, walk_num do j = 1, elec_num do i = 1, elec_num @@ -1936,7 +1905,7 @@ integer function qmckl_compute_factor_ee_deriv_e_f(context, walk_num, elec_num, (i .GT. up_num .AND. j .GT. up_num)) then spin_fact = 0.5d0 endif - + lap1 = 0.0d0 lap2 = 0.0d0 lap3 = 0.0d0 @@ -1950,7 +1919,7 @@ integer function qmckl_compute_factor_ee_deriv_e_f(context, walk_num, elec_num, lap2 = lap2 + y x = x * ee_distance_rescaled(nw, i, j) end do - + lap3 = lap3 - 2.0d0 * bord_vector(2) * dx(ii) * dx(ii) factor_ee_deriv_e( j, ii, nw) = factor_ee_deriv_e( j, ii, nw) + spin_fact * bord_vector(1) * & @@ -1984,7 +1953,7 @@ end function qmckl_compute_factor_ee_deriv_e_f const double* ee_distance_rescaled, const double* ee_distance_rescaled_deriv_e, const double* asymp_jasb, - double* const factor_ee_deriv_e ); + double* const factor_ee_deriv_e ); #+end_src @@ -2097,7 +2066,7 @@ for j in range(elec_num): if((i <= (up_num-1) and j <= (up_num-1) ) or \ (i > (up_num-1) and j > (up_num-1))): spin_fact = 0.5 - + lap1 = 0.0 lap2 = 0.0 lap3 = 0.0 @@ -2111,7 +2080,7 @@ for j in range(elec_num): lap1 = lap1 + (p - 1) * y * xinv * dx[ii] * dx[ii] lap2 = lap2 + y x = x * ee_distance_rescaled[j][i] - + lap3 = lap3 - 2.0 * bord_vector[1] * dx[ii] * dx[ii] factor_ee_deriv_e[ii][j] = factor_ee_deriv_e[ii][j] + spin_fact * bord_vector[0] * \ @@ -2137,7 +2106,7 @@ print("factor_ee_deriv_e[3][0]:",factor_ee_deriv_e[3][0]) : factor_ee_deriv_e[1][0]: -0.6927548119830084 : factor_ee_deriv_e[2][0]: 0.073267755223968 : factor_ee_deriv_e[3][0]: 1.5111672803213185 - + #+begin_src c :tangle (eval c_test) /* Check if Jastrow is properly initialized */ @@ -2158,13 +2127,13 @@ assert(fabs(factor_ee_deriv_e[0][3][0]-1.5111672803213185 ) < 1.e-12); ** 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~. + coeffecients and the electron-nucleus rescaled distances ~en_distance_rescaled~. \[ f_{en} = \sum_{i,jjastrow.factor_en = factor_en; } - + qmckl_exit_code rc = qmckl_compute_factor_en(context, ctx->electron.walk_num, @@ -2320,7 +2289,7 @@ integer function qmckl_compute_factor_en_f(context, walk_num, elec_num, nucl_num endif factor_en = 0.0d0 - + do nw =1, walk_num do a = 1, nucl_num do i = 1, elec_num @@ -2359,7 +2328,7 @@ end function qmckl_compute_factor_en_f const int64_t aord_num, const double* aord_vector, const double* en_distance_rescaled, - double* const factor_en ); + double* const factor_en ); #+end_src @@ -2421,7 +2390,7 @@ 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 @@ -2432,10 +2401,10 @@ for a in range(0,nucl_num): 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 */ @@ -2455,7 +2424,7 @@ assert(fabs(factor_en[0]+5.865822569188727) < 1.e-12); ~en_distance_rescaled_deriv_e~ which are already calculated previously. TODO: write equations. - + *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_jastrow_factor_en_deriv_e(qmckl_context context, double* const factor_en_deriv_e); @@ -2487,7 +2456,7 @@ qmckl_exit_code qmckl_get_jastrow_factor_en_deriv_e(qmckl_context context, doubl #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_provide_factor_en_deriv_e(qmckl_context context); #+end_src - + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_provide_factor_en_deriv_e(qmckl_context context) { @@ -2527,7 +2496,7 @@ qmckl_exit_code qmckl_provide_factor_en_deriv_e(qmckl_context context) } ctx->jastrow.factor_en_deriv_e = factor_en_deriv_e; } - + qmckl_exit_code rc = qmckl_compute_factor_en_deriv_e(context, ctx->electron.walk_num, @@ -2550,7 +2519,7 @@ qmckl_exit_code qmckl_provide_factor_en_deriv_e(qmckl_context context) return QMCKL_SUCCESS; } #+end_src - + *** Compute :PROPERTIES: :Name: qmckl_compute_factor_en_deriv_e @@ -2623,7 +2592,7 @@ integer function qmckl_compute_factor_en_deriv_e_f(context, walk_num, elec_num, factor_en_deriv_e = 0.0d0 third = 1.0d0 / 3.0d0 - + do nw =1, walk_num do a = 1, nucl_num do i = 1, elec_num @@ -2689,7 +2658,7 @@ end function qmckl_compute_factor_en_deriv_e_f const double* aord_vector, const double* en_distance_rescaled, const double* en_distance_rescaled_deriv_e, - double* const factor_en_deriv_e ); + double* const factor_en_deriv_e ); #+end_src @@ -2811,7 +2780,7 @@ for a in range(nucl_num): lap1 = lap1 + (p - 1) * y * xinv * dx[ii] * dx[ii] lap2 = lap2 + y x = x * en_distance_rescaled[i][a] - + lap3 = lap3 - 2.0 * aord_vector[1][type_nucl_vector[a]-1] * dx[ii] * dx[ii] factor_en_deriv_e[ii][i] = factor_en_deriv_e[ii][i] + aord_vector[0][type_nucl_vector[a]-1] * \ @@ -2836,7 +2805,7 @@ print("factor_en_deriv_e[3][0]:",factor_en_deriv_e[3][0]) : factor_en_deriv_e[1][0]: -0.23301394780804574 : factor_en_deriv_e[2][0]: 0.17548337641865783 : factor_en_deriv_e[3][0]: -0.9667363412285741 - + #+begin_src c :tangle (eval c_test) /* Check if Jastrow is properly initialized */ @@ -2853,7 +2822,7 @@ assert(fabs(factor_en_deriv_e[0][2][0]-0.17548337641865783) < 1.e-12); 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 @@ -2864,7 +2833,7 @@ assert(fabs(factor_en_deriv_e[0][3][0]+0.9667363412285741 ) < 1.e-12); \] where \(C_{ij}\) is the matrix of electron-electron distances. - + *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes @@ -2892,7 +2861,7 @@ qmckl_exit_code qmckl_get_jastrow_een_rescaled_e(qmckl_context context, double* return QMCKL_SUCCESS; } #+end_src - + *** Provide :noexport: #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none @@ -3050,7 +3019,7 @@ integer function qmckl_compute_een_rescaled_e_f(context, walk_num, elec_num, cor do l = 0, cord_num do j = 1, elec_num een_rescaled_e(l, j, j, nw) = 0.0d0 - end do + end do end do end do @@ -3069,7 +3038,7 @@ end function qmckl_compute_een_rescaled_e_f const int64_t cord_num, const double rescale_factor_kappa_ee, const double* ee_distance, - double* const een_rescaled_e ); + 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")) @@ -3176,15 +3145,15 @@ assert(fabs(een_rescaled_e[0][1][5][2]-0.3424402276009091) < 1.e-12); #+end_src ** Electron-electron rescaled distances for each order and derivatives - + ~een_rescaled_e_deriv_e~ stores the table of the derivatives of the rescaled distances between all pairs of electrons and raised to the power \(p\) defined by ~cord_num~. Here we take its derivatives required for the een jastrow. - + TODO: write formulae - + *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes @@ -3212,7 +3181,7 @@ qmckl_exit_code qmckl_get_jastrow_een_rescaled_e_deriv_e(qmckl_context context, return QMCKL_SUCCESS; } #+end_src - + *** Provide :noexport: #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none @@ -3391,7 +3360,7 @@ end function qmckl_compute_factor_een_rescaled_e_deriv_e_f const double* coord_new, const double* ee_distance, const double* een_rescaled_e, - double* const een_rescaled_e_deriv_e ); + double* const een_rescaled_e_deriv_e ); #+end_src @@ -3495,11 +3464,11 @@ for l in range(0,cord_num+1): for j in range(0,elec_num): for i in range(0,elec_num): for ii in range(0,4): - een_rescaled_e_deriv_e[i,ii,j,l] = kappa_l * elec_dist_deriv_e[ii,i,j] + een_rescaled_e_deriv_e[i,ii,j,l] = kappa_l * elec_dist_deriv_e[ii,i,j] een_rescaled_e_deriv_e[i,3,j,l] = een_rescaled_e_deriv_e[i,3,j,l] + \ een_rescaled_e_deriv_e[i,0,j,l] * een_rescaled_e_deriv_e[i,0,j,l] + \ een_rescaled_e_deriv_e[i,1,j,l] * een_rescaled_e_deriv_e[i,1,j,l] + \ - een_rescaled_e_deriv_e[i,2,j,l] * een_rescaled_e_deriv_e[i,2,j,l] + een_rescaled_e_deriv_e[i,2,j,l] * een_rescaled_e_deriv_e[i,2,j,l] for ii in range(0,4): een_rescaled_e_deriv_e[i,ii,j,l] = een_rescaled_e_deriv_e[i,ii,j,l] * een_rescaled_e[i,j,l] @@ -3535,8 +3504,8 @@ assert(fabs(een_rescaled_e_deriv_e[0][1][0][5][2] + 0.5880599146214673 ) < 1. #+end_src ** Electron-nucleus rescaled distances for each order - - ~een_rescaled_n~ stores the table of the rescaled distances between + + ~een_rescaled_n~ stores the table of the rescaled distances between electrons and nucleii raised to the power \(p\) defined by ~cord_num~: \[ @@ -3544,7 +3513,7 @@ assert(fabs(een_rescaled_e_deriv_e[0][1][0][5][2] + 0.5880599146214673 ) < 1. \] where \(C_{ia}\) is the matrix of electron-nucleus distances. - + *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes @@ -3572,7 +3541,7 @@ qmckl_exit_code qmckl_get_jastrow_een_rescaled_n(qmckl_context context, double* return QMCKL_SUCCESS; } #+end_src - + *** Provide :noexport: #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none @@ -3734,7 +3703,7 @@ end function qmckl_compute_een_rescaled_n_f const int64_t cord_num, const double rescale_factor_kappa_en, const double* en_distance, - double* const een_rescaled_n ); + double* const een_rescaled_n ); #+end_src #+CALL: generate_c_interface(table=qmckl_factor_een_rescaled_n_args,rettyp=get_value("CRetType"),fname=get_value("Name")) @@ -3840,10 +3809,10 @@ assert(fabs(een_rescaled_n[0][5][1][2]-0.01343938025140174) < 1.e-12); ** Electron-nucleus rescaled distances for each order and derivatives - ~een_rescaled_n_deriv_e~ stores the table of the rescaled distances between + ~een_rescaled_n_deriv_e~ stores the table of the rescaled distances between electrons and nucleii raised to the power \(p\) defined by ~cord_num~: - + *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes @@ -3871,7 +3840,7 @@ qmckl_exit_code qmckl_get_jastrow_een_rescaled_n_deriv_e(qmckl_context context, return QMCKL_SUCCESS; } #+end_src - + *** Provide :noexport: #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none @@ -4068,11 +4037,11 @@ end function qmckl_compute_factor_een_rescaled_n_deriv_e_f const double* coord, const double* en_distance, const double* een_rescaled_n, - double* const een_rescaled_n_deriv_e ); + double* const een_rescaled_n_deriv_e ); #+end_src #+CALL: generate_c_interface(table=qmckl_compute_factor_een_rescaled_n_deriv_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_factor_een_rescaled_n_deriv_e & @@ -4163,11 +4132,11 @@ for l in range(0,cord_num+1): for j in range(0,elec_num): for a in range(0,nucl_num): for ii in range(0,4): - een_rescaled_n_deriv_e[j,ii,a,l] = kappa_l * elnuc_dist_deriv_e[ii,j,a] + een_rescaled_n_deriv_e[j,ii,a,l] = kappa_l * elnuc_dist_deriv_e[ii,j,a] een_rescaled_n_deriv_e[j,3,a,l] = een_rescaled_n_deriv_e[j,3,a,l] + \ een_rescaled_n_deriv_e[j,0,a,l] * een_rescaled_n_deriv_e[j,0,a,l] + \ een_rescaled_n_deriv_e[j,1,a,l] * een_rescaled_n_deriv_e[j,1,a,l] + \ - een_rescaled_n_deriv_e[j,2,a,l] * een_rescaled_n_deriv_e[j,2,a,l] + een_rescaled_n_deriv_e[j,2,a,l] * een_rescaled_n_deriv_e[j,2,a,l] for ii in range(0,4): een_rescaled_n_deriv_e[j,ii,a,l] = een_rescaled_n_deriv_e[j,ii,a,l] * een_rescaled_n[a,j,l] @@ -4208,7 +4177,7 @@ assert(fabs(een_rescaled_n_deriv_e[0][5][0][1][2]-0.001593334817691633 ) < 1.e Prepare ~cord_vect_full~ and ~lkpm_combined_index~ tables required for the calculation of the three-body jastrow ~factor_een~ and its derivative - ~factor_een_deriv_e~. + ~factor_een_deriv_e~. *** Get @@ -4306,7 +4275,7 @@ qmckl_exit_code qmckl_get_jastrow_tmp_c(qmckl_context context, double* const tmp qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - size_t sze = (ctx->jastrow.cord_num) * (ctx->jastrow.cord_num + 1) + size_t sze = (ctx->jastrow.cord_num) * (ctx->jastrow.cord_num + 1) * ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num; memcpy(tmp_c, ctx->jastrow.tmp_c, sze * sizeof(double)); @@ -4333,14 +4302,14 @@ qmckl_exit_code qmckl_get_jastrow_dtmp_c(qmckl_context context, double* const dt qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - size_t sze = (ctx->jastrow.cord_num) * (ctx->jastrow.cord_num + 1) + size_t sze = (ctx->jastrow.cord_num) * (ctx->jastrow.cord_num + 1) *4* ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num; memcpy(dtmp_c, ctx->jastrow.dtmp_c, sze * sizeof(double)); return QMCKL_SUCCESS; } #+end_src - + *** Provide :noexport: #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none @@ -4499,7 +4468,7 @@ qmckl_exit_code qmckl_provide_tmp_c(qmckl_context context) if (ctx->jastrow.tmp_c == NULL) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; - mem_info.size = (ctx->jastrow.cord_num) * (ctx->jastrow.cord_num + 1) + mem_info.size = (ctx->jastrow.cord_num) * (ctx->jastrow.cord_num + 1) ,* ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num * sizeof(double); double* tmp_c = (double*) qmckl_malloc(context, mem_info); @@ -4552,7 +4521,7 @@ qmckl_exit_code qmckl_provide_dtmp_c(qmckl_context context) if (ctx->jastrow.dtmp_c == NULL) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; - mem_info.size = (ctx->jastrow.cord_num) * (ctx->jastrow.cord_num + 1) + mem_info.size = (ctx->jastrow.cord_num) * (ctx->jastrow.cord_num + 1) * 4 * ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num * sizeof(double); double* dtmp_c = (double*) qmckl_malloc(context, mem_info); @@ -4648,7 +4617,7 @@ end function qmckl_compute_dim_cord_vect_f qmckl_exit_code qmckl_compute_dim_cord_vect ( const qmckl_context context, const int64_t cord_num, - int64_t* const dim_cord_vect ); + int64_t* const dim_cord_vect ); #+end_src @@ -4691,7 +4660,7 @@ end function qmckl_compute_dim_cord_vect_f | ~type_nucl_vector~ | ~int64_t[nucl_num]~ | in | dimension of cord full table | | ~cord_vector~ | ~double[dim_cord_vect][type_nucl_num]~ | in | dimension of cord full table | | ~cord_vect_full~ | ~double[dim_cord_vect][nucl_num]~ | out | Full list of coefficients | - + #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_cord_vect_full_f(context, nucl_num, dim_cord_vect, type_nucl_num, & type_nucl_vector, cord_vector, cord_vect_full) & @@ -4749,7 +4718,7 @@ end function qmckl_compute_cord_vect_full_f const int64_t type_nucl_num, const int64_t* type_nucl_vector, const double* cord_vector, - double* const cord_vect_full ); + double* const cord_vect_full ); #+end_src @@ -4856,7 +4825,7 @@ end function qmckl_compute_lkpm_combined_index_f const qmckl_context context, const int64_t cord_num, const int64_t dim_cord_vect, - int64_t* const lpkm_combined_index ); + int64_t* const lpkm_combined_index ); #+end_src @@ -5013,7 +4982,7 @@ end function qmckl_compute_tmp_c_f const int64_t walk_num, const double* een_rescaled_e, const double* een_rescaled_n, - double* const tmp_c ); + double* const tmp_c ); #+end_src @@ -5167,7 +5136,7 @@ end function qmckl_compute_dtmp_c_f const int64_t walk_num, const double* een_rescaled_e_deriv_e, const double* een_rescaled_n, - double* const dtmp_c ); + double* const dtmp_c ); #+end_src @@ -5283,12 +5252,12 @@ assert(fabs(dtmp_c[0][1][0][0][0][0] + 0.237440520852232) < 1e-12); #+end_src ** Electron-electron-nucleus Jastrow \(f_{een}\) - - Calculate the electron-electron-nuclear three-body jastrow component ~factor_een~ + + Calculate the electron-electron-nuclear three-body jastrow component ~factor_een~ using the above prepared tables. TODO: write equations. - + *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_jastrow_factor_een(qmckl_context context, double* const factor_een); @@ -5320,7 +5289,7 @@ qmckl_exit_code qmckl_get_jastrow_factor_een(qmckl_context context, double* cons #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_provide_factor_een(qmckl_context context); #+end_src - + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_provide_factor_een(qmckl_context context) { @@ -5372,7 +5341,7 @@ qmckl_exit_code qmckl_provide_factor_een(qmckl_context context) } ctx->jastrow.factor_een = factor_een; } - + qmckl_exit_code rc = qmckl_compute_factor_een(context, ctx->electron.walk_num, @@ -5510,7 +5479,7 @@ end function qmckl_compute_factor_een_naive_f const int64_t* lkpm_combined_index, const double* een_rescaled_e, const double* een_rescaled_n, - double* const factor_een ); + double* const factor_een ); #+end_src #+CALL: generate_c_interface(table=qmckl_factor_een_naive_args,rettyp=get_value("CRetType"),fname=get_value("Name")) @@ -5584,7 +5553,7 @@ end function qmckl_compute_factor_een_naive_f | ~tmp_c~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]~ | vector of non-zero coefficients | | | ~een_rescaled_n~ | ~double[walk_num][elec_num][nucl_num][0:cord_num]~ | in | Electron-nucleus rescaled factor | | ~factor_een~ | ~double[walk_num]~ | out | Electron-nucleus jastrow | - + #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_factor_een_f(context, walk_num, elec_num, nucl_num, cord_num, & dim_cord_vect, cord_vect_full, lkpm_combined_index, & @@ -5670,7 +5639,7 @@ end function qmckl_compute_factor_een_f const int64_t* lkpm_combined_index, const double* een_rescaled_e, const double* een_rescaled_n, - double* const factor_een ); + double* const factor_een ); #+end_src #+CALL: generate_c_interface(table=qmckl_factor_een_args,rettyp=get_value("CRetType"),fname=get_value("Name")) @@ -5757,7 +5726,7 @@ print("factor_een:",factor_een) #+RESULTS: : factor_een: -0.37407972141304213 - + #+begin_src c :tangle (eval c_test) /* Check if Jastrow is properly initialized */ @@ -5768,14 +5737,14 @@ rc = qmckl_get_jastrow_factor_een(context, &(factor_een[0])); assert(fabs(factor_een[0] + 0.37407972141304213) < 1e-12); #+end_src - + ** Electron-electron-nucleus Jastrow \(f_{een}\) derivative - Calculate the electron-electron-nuclear three-body jastrow component ~factor_een_deriv_e~ + Calculate the electron-electron-nuclear three-body jastrow component ~factor_een_deriv_e~ using the above prepared tables. TODO: write equations. - + *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_jastrow_factor_een_deriv_e(qmckl_context context, double* const factor_een_deriv_e); @@ -5807,7 +5776,7 @@ qmckl_exit_code qmckl_get_jastrow_factor_een_deriv_e(qmckl_context context, doub #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_provide_factor_een_deriv_e(qmckl_context context); #+end_src - + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_provide_factor_een_deriv_e(qmckl_context context) { @@ -5871,7 +5840,7 @@ qmckl_exit_code qmckl_provide_factor_een_deriv_e(qmckl_context context) } ctx->jastrow.factor_een_deriv_e = factor_een_deriv_e; } - + qmckl_exit_code rc = qmckl_compute_factor_een_deriv_e(context, ctx->electron.walk_num, @@ -5971,13 +5940,13 @@ integer function qmckl_compute_factor_een_deriv_e_naive_f(context, walk_num, ele endif factor_een_deriv_e = 0.0d0 - + do nw =1, walk_num do n = 1, dim_cord_vect - l = lkpm_combined_index(n, 1) - k = lkpm_combined_index(n, 2) - p = lkpm_combined_index(n, 3) - m = lkpm_combined_index(n, 4) + l = lkpm_combined_index(n, 1) + k = lkpm_combined_index(n, 2) + p = lkpm_combined_index(n, 3) + m = lkpm_combined_index(n, 4) do a = 1, nucl_num cn = cord_vect_full(a, n) @@ -6009,7 +5978,7 @@ integer function qmckl_compute_factor_een_deriv_e_naive_f(context, walk_num, ele daccu2(1) * een_rescaled_n_deriv_e(m, a, 1, j, nw ) + & daccu2(2) * een_rescaled_n_deriv_e(m, a, 2, j, nw ) + & daccu2(3) * een_rescaled_n_deriv_e(m, a, 3, j, nw ) ) * cn - + end do end do end do @@ -6035,7 +6004,7 @@ end function qmckl_compute_factor_een_deriv_e_naive_f const double* een_rescaled_n, const double* een_rescaled_e_deriv_e, const double* een_rescaled_n_deriv_e, - double* const factor_een_deriv_e ); + double* const factor_een_deriv_e ); #+end_src @@ -6118,8 +6087,8 @@ end function qmckl_compute_factor_een_deriv_e_naive_f | ~een_rescaled_n~ | ~double[walk_num][elec_num][nucl_num][0:cord_num]~ | in | Electron-nucleus rescaled factor | | ~een_rescaled_n_deriv_e~ | ~double[walk_num][elec_num][4][nucl_num][0:cord_num]~ | in | Derivative of Electron-nucleus rescaled factor | | ~factor_een_deriv_e~ | ~double[walk_num][4][elec_num]~ | out | Derivative of Electron-nucleus jastrow | - - + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_factor_een_deriv_e_f(context, walk_num, elec_num, nucl_num, cord_num, dim_cord_vect, & cord_vect_full, lkpm_combined_index, & @@ -6226,7 +6195,7 @@ end function qmckl_compute_factor_een_deriv_e_f const double* dtmp_c, const double* een_rescaled_n, const double* een_rescaled_n_deriv_e, - double* const factor_een_deriv_e ); + double* const factor_een_deriv_e ); #+end_src @@ -6336,7 +6305,7 @@ print("factor_een:",factor_een) #+RESULTS: : (6, 10, 4, 10) : factor_een: 0.0 - + #+begin_src c :tangle (eval c_test) /* Check if Jastrow is properly initialized */ @@ -6347,7 +6316,7 @@ rc = qmckl_get_jastrow_factor_een_deriv_e(context, &(factor_een_deriv_e[0][0])); assert(fabs(factor_een_deriv_e[0][0] + 0.0005481671107226865) < 1e-12); #+end_src - + * End of files :noexport: #+begin_src c :tangle (eval h_private_type) diff --git a/share/qmckl/test_data/chbrclf/metadata.txt b/share/qmckl/test_data/chbrclf/metadata.txt index d0180b3..834b9f7 100644 --- a/share/qmckl/test_data/chbrclf/metadata.txt +++ b/share/qmckl/test_data/chbrclf/metadata.txt @@ -9,3 +9,4 @@ len_metadata_description 0 metadata_description metadata_code metadata_author +metadata_unsafe_isSet 0