diff --git a/Makefile.am b/Makefile.am index 57ad20b..ca045b9 100644 --- a/Makefile.am +++ b/Makefile.am @@ -52,7 +52,7 @@ test_qmckl_f = tests/qmckl_f.f90 test_qmckl_fo = tests/qmckl_f.o src_qmckl_f = src/qmckl_f.f90 src_qmckl_fo = src/qmckl_f.o -header_tests = tests/chbrclf.h +header_tests = tests/chbrclf.h tests/n2.h fortrandir = $(datadir)/$(PACKAGE_NAME)/fortran/ dist_fortran_DATA = $(qmckl_f) @@ -157,6 +157,7 @@ $(htmlize_el): tests/chbrclf.h: $(qmckl_h) +tests/n2.h: $(qmckl_h) generated.mk: $(ORG_FILES) $(PYTHON) $(srcdir)/tools/build_makefile.py diff --git a/configure.ac b/configure.ac index c313f56..a2a2438 100644 --- a/configure.ac +++ b/configure.ac @@ -98,7 +98,6 @@ AC_CHECK_LIB([pthread], [pthread_create]) # CFLAGS="${CFLAGS} ${OPENMP_CFLAGS}" #fi - ## BLAS #AX_BLAS([], [AC_MSG_ERROR([BLAS was not found.])]) @@ -221,13 +220,13 @@ ${PACKAGE_NAME} Version ${PACKAGE_VERSION} ${QMCKL_DEVEL} Prefix: '${prefix}'. -CC..........: ${CC} -CPPFLAGS....: ${CPPFLAGS} -CFLAGS......: ${CFLAGS} -FC..........: ${FC} -FCLAGS......: ${FCFLAGS} -LDFLAGS:....: ${LDFLAGS} -LIBS........: ${LIBS} +CC..............: ${CC} +CPPFLAGS........: ${CPPFLAGS} +CFLAGS..........: ${CFLAGS} +FC..............: ${FC} +FCLAGS..........: ${FCFLAGS} +LDFLAGS:........: ${LDFLAGS} +LIBS............: ${LIBS} Package features: ${ARGS} diff --git a/org/qmckl_ao.org b/org/qmckl_ao.org index 21d9b72..0b10745 100644 --- a/org/qmckl_ao.org +++ b/org/qmckl_ao.org @@ -1921,6 +1921,13 @@ qmckl_exit_code qmckl_provide_ao_basis_shell_vgl(qmckl_context context) "qmckl_ao_basis_shell_vgl", NULL); } + + if(!(ctx->electron.provided)) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_electron", + NULL); + } /* Compute if necessary */ if (ctx->electron.coord_new_date > ctx->ao_basis.shell_vgl_date) { diff --git a/org/qmckl_context.org b/org/qmckl_context.org index 11eb235..d2ca467 100644 --- a/org/qmckl_context.org +++ b/org/qmckl_context.org @@ -31,6 +31,7 @@ int main() { #include "qmckl_nucleus_private_type.h" #include "qmckl_electron_private_type.h" #include "qmckl_ao_private_type.h" +#include "qmckl_jastrow_private_type.h" #include "qmckl_nucleus_private_func.h" #include "qmckl_electron_private_func.h" #include "qmckl_ao_private_func.h" @@ -118,6 +119,7 @@ typedef struct qmckl_context_struct { qmckl_nucleus_struct nucleus; qmckl_electron_struct electron; qmckl_ao_basis_struct ao_basis; + qmckl_jastrow_struct jastrow; /* To be implemented: qmckl_mo_struct mo; 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_electron.org b/org/qmckl_electron.org index 6587f4a..945a214 100644 --- a/org/qmckl_electron.org +++ b/org/qmckl_electron.org @@ -63,30 +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][num][num]~ | Electron-electron rescaled distances derivatives | - | ~en_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 @@ -1566,7 +1565,6 @@ rc = qmckl_get_electron_ee_distance_rescaled_deriv_e(context, ee_distance_rescal #+end_src - ** Electron-nucleus distances *** Get @@ -2140,12 +2138,12 @@ assert(fabs(en_distance_rescaled[1][0][1] - 0.9584331688679852) < 1.e-12); The rescaled distances which is given as $R = (1 - \exp{-\kappa r})/\kappa$ needs to be perturbed with respect to the nuclear coordinates. This data is stored in the ~en_distance_rescaled_deriv_e~ tensor. The - The first three elements of this three index tensor ~[4][num][num]~ gives the + The first three elements of this three index tensor ~[4][nucl_num][elec_num]~ gives the derivatives in the x, y, and z directions $dx, dy, dz$ and the last index gives the Laplacian $\partial x^2 + \partial y^2 + \partial z^2$. *** Get - + #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_electron_en_distance_rescaled_deriv_e(qmckl_context context, double* distance_rescaled_deriv_e); #+end_src diff --git a/org/qmckl_jastrow.org b/org/qmckl_jastrow.org new file mode 100644 index 0000000..118b6f2 --- /dev/null +++ b/org/qmckl_jastrow.org @@ -0,0 +1,5509 @@ +#+TITLE: Jastrow Factor +#+SETUPFILE: ../tools/theme.setup +#+INCLUDE: ../tools/lib.org + +Functions for the calculation of the Jastrow factor \(f_{ee}, f_{en}, f_{een}\). +These are stored in the ~factor_ee~, ~factor_en~, and ~factor_een~ variables. +The ~jastrow~ structure contains all the information required to build +these factors along with their derivatives. + +* Headers :noexport: + #+begin_src elisp :noexport :results none +(org-babel-lob-ingest "../tools/lib.org") +#+end_src + + + #+begin_src c :tangle (eval h_private_type) +#ifndef QMCKL_JASTROW_HPT +#define QMCKL_JASTROW_HPT +#include + #+end_src + + #+begin_src c :tangle (eval c_test) :noweb yes +#include "qmckl.h" +#include +#include +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include +#include "n2.h" + +int main() { + qmckl_context context; + context = qmckl_context_create(); + #+end_src + + #+begin_src c :tangle (eval c) +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#ifdef HAVE_STDINT_H +#include +#elif HAVE_INTTYPES_H +#include +#endif + +#include +#include +#include +#include +#include + +#include + +#include "qmckl.h" +#include "qmckl_context_private_type.h" +#include "qmckl_memory_private_type.h" +#include "qmckl_memory_private_func.h" +#include "qmckl_jastrow_private_func.h" +#include "qmckl_jastrow_private_type.h" + #+end_src + +* Context + :PROPERTIES: + :Name: qmckl_jastrow + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + 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 | + | ~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 | + | ~uint64_t~ | ~factor_ee_date~ | out | Jastrow factor: electron-electron part | + | ~double~ | ~factor_en[walk_num]~ | out | Jastrow factor: electron-nucleus part | + | ~uint64_t~ | ~factor_en_date~ | out | Jastrow factor: electron-nucleus part | + | ~double~ | ~factor_een[walk_num]~ | out | Jastrow factor: electron-electron-nucleus part | + | ~uint64_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 | + | ~uint64_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 | + | ~uint64_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 | + | ~uint64_t~ | ~factor_een_deriv_e_date~ | out | Keep track of the date for the een derivative | + + computed data: + + |------------+-----------------------------------------------------------------------+---------------------------------------------------------------------------------------------------------| + | ~int64_t~ | ~dim_cord_vect~ | Number of unique C coefficients | + | ~uint64_t~ | ~dim_cord_vect_date~ | Number of unique C coefficients | + | ~double~ | ~asymp_jasb[2]~ | Asymptotic component | + | ~uint64_t~ | ~asymp_jasb_date~ | Asymptotic component | + | ~double~ | ~cord_vect_full[dim_cord_vect][nucl_num]~ | vector of non-zero coefficients | + | ~uint64_t~ | ~cord_vect_full_date~ | Keep track of changes here | + | ~int64_t~ | ~lkpm_combined_index[4][dim_cord_vect]~ | Transform l,k,p, and m into consecutive indices | + | ~uint64_t~ | ~lkpm_combined_index_date~ | 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 | + | ~uint64_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 | + | ~uint64_t~ | ~een_rescaled_n_date~ | Keep track of the date of creation | + | ~double~ | ~een_rescaled_e_deriv_e[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 | + | ~uint64_t~ | ~een_rescaled_e_deriv_e_date~ | Keep track of the date of creation | + | ~double~ | ~een_rescaled_n_deriv_e[walk_num][elec_num][4][nucl_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons | + | ~uint64_t~ | ~een_rescaled_n_deriv_e_date~ | Keep track of the date of creation | + + For H2O we have the following data: + + #+NAME: jastrow_data + #+BEGIN_SRC python :results output +import numpy as np + +elec_num = 10 +nucl_num = 2 +up_num = 5 +down_num = 5 +nucl_coord = np.array([ [0.000000, 0.000000 ], + [0.000000, 0.000000 ], + [0.000000, 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]] + +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): + ee_distance_rescaled[i][j] = ee_distance_rescaled[j][i] + +type_nucl_num = 1 +aord_num = 5 +bord_num = 5 +cord_num = 23 +dim_cord_vect= 23 +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, + 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, + 1.454340000000000E-002, 5.143510000000000E-002, 9.250000000000000E-004, + -4.099100000000000E-003, 4.327600000000000E-003, -1.654470000000000E-003, + 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, + 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, + 1.454340000000000E-002, 5.143510000000000E-002, 9.250000000000000E-004, + -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, + 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, + 1.454340000000000E-002, 5.143510000000000E-002, 9.250000000000000E-004, + -4.099100000000000E-003, 4.327600000000000E-003, -1.654470000000000E-003, + 2.614000000000000E-003, -1.477000000000000E-003, -1.137000000000000E-003, + -4.010475000000000E-002, 6.106710000000000E-003 ], +] +lkpm_combined_index = [[1 , 1 , 2 , 0], + [0 , 0 , 2 , 1], + [1 , 2 , 3 , 0], + [2 , 1 , 3 , 0], + [0 , 1 , 3 , 1], + [1 , 0 , 3 , 1], + [1 , 3 , 4 , 0], + [2 , 2 , 4 , 0], + [0 , 2 , 4 , 1], + [3 , 1 , 4 , 0], + [1 , 1 , 4 , 1], + [2 , 0 , 4 , 1], + [0 , 0 , 4 , 2], + [1 , 4 , 5 , 0], + [2 , 3 , 5 , 0], + [0 , 3 , 5 , 1], + [3 , 2 , 5 , 0], + [1 , 2 , 5 , 1], + [4 , 1 , 5 , 0], + [2 , 1 , 5 , 1], + [0 , 1 , 5 , 2], + [3 , 0 , 5 , 1], + [1 , 0 , 5 , 2]] + +kappa = 1.0 +kappa_inv = 1.0/kappa + #+END_SRC + + #+RESULTS: jastrow_data + +** Data structure + + #+begin_src c :comments org :tangle (eval h_private_type) +typedef struct qmckl_jastrow_struct{ + int32_t uninitialized; + int64_t aord_num; + int64_t bord_num; + int64_t cord_num; + int64_t type_nucl_num; + uint64_t asymp_jasb_date; + uint64_t tmp_c_date; + uint64_t dtmp_c_date; + uint64_t factor_ee_date; + uint64_t factor_en_date; + uint64_t factor_een_date; + uint64_t factor_ee_deriv_e_date; + uint64_t factor_en_deriv_e_date; + uint64_t factor_een_deriv_e_date; + int64_t* type_nucl_vector; + double * aord_vector; + double * bord_vector; + double * cord_vector; + double * asymp_jasb; + double * factor_ee; + double * factor_en; + double * factor_een; + double * factor_ee_deriv_e; + double * factor_en_deriv_e; + double * factor_een_deriv_e; + int64_t dim_cord_vect; + uint64_t dim_cord_vect_date; + double * cord_vect_full; + uint64_t cord_vect_full_date; + int64_t* lkpm_combined_index; + uint64_t lkpm_combined_index_date; + double * tmp_c; + double * dtmp_c; + double * een_rescaled_e; + double * een_rescaled_n; + uint64_t een_rescaled_e_date; + uint64_t een_rescaled_n_date; + double * een_rescaled_e_deriv_e; + double * een_rescaled_n_deriv_e; + uint64_t een_rescaled_e_deriv_e_date; + uint64_t een_rescaled_n_deriv_e_date; + bool provided; + char * type; +} qmckl_jastrow_struct; + #+end_src + + + The ~uninitialized~ integer contains one bit set to one for each + initialization function which has not been called. It becomes equal + to zero after all initialization functions have been called. The + struct is then initialized and ~provided == true~. + Some values are initialized by default, and are not concerned by + this mechanism. + + #+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) { + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return false; + } + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + ctx->jastrow.uninitialized = (1 << 6) - 1; + + /* Default values */ + + return QMCKL_SUCCESS; +} + #+end_src + +** 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_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 + requires the following additional information to be set: + + + When all the data for the AOs have been provided, the following + function returns ~true~. + + #+begin_src c :comments org :tangle (eval h_func) +bool qmckl_jastrow_provided (const qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +bool qmckl_jastrow_provided(const qmckl_context context) { + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return false; + } + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + return ctx->jastrow.provided; +} + #+end_src + + #+NAME:post + #+begin_src c :exports none +if ( (ctx->jastrow.uninitialized & mask) != 0) { + return NULL; +} + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_get_jastrow_aord_num (const qmckl_context context, int64_t* const aord_num) { + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return (char) 0; + } + + if (aord_num == NULL) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_get_jastrow_aord_num", + "aord_num is a null pointer"); + } + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + int32_t mask = 1 << 0; + + if ( (ctx->jastrow.uninitialized & mask) != 0) { + return QMCKL_NOT_PROVIDED; + } + + assert (ctx->jastrow.aord_num > 0); + *aord_num = ctx->jastrow.aord_num; + return QMCKL_SUCCESS; +} + +qmckl_exit_code qmckl_get_jastrow_bord_num (const qmckl_context context, int64_t* const bord_num) { + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return (char) 0; + } + + if (bord_num == NULL) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_get_jastrow_bord_num", + "aord_num is a null pointer"); + } + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + int32_t mask = 1 << 0; + + if ( (ctx->jastrow.uninitialized & mask) != 0) { + return QMCKL_NOT_PROVIDED; + } + + assert (ctx->jastrow.bord_num > 0); + *bord_num = ctx->jastrow.bord_num; + return QMCKL_SUCCESS; +} + +qmckl_exit_code qmckl_get_jastrow_cord_num (const qmckl_context context, int64_t* const cord_num) { + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return (char) 0; + } + + if (cord_num == NULL) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_get_jastrow_cord_num", + "aord_num is a null pointer"); + } + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + int32_t mask = 1 << 0; + + if ( (ctx->jastrow.uninitialized & mask) != 0) { + return QMCKL_NOT_PROVIDED; + } + + assert (ctx->jastrow.cord_num > 0); + *cord_num = ctx->jastrow.cord_num; + return QMCKL_SUCCESS; +} + +qmckl_exit_code qmckl_get_jastrow_type_nucl_num (const qmckl_context context, int64_t* const type_nucl_num) { + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return (char) 0; + } + + if (type_nucl_num == NULL) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_get_jastrow_type_nucl_num", + "type_nucl_num is a null pointer"); + } + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + int32_t mask = 1 << 1; + + if ( (ctx->jastrow.uninitialized & mask) != 0) { + return QMCKL_NOT_PROVIDED; + } + + assert (ctx->jastrow.type_nucl_num > 0); + *type_nucl_num = ctx->jastrow.type_nucl_num; + 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) { + return (char) 0; + } + + if (aord_vector == NULL) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_get_jastrow_aord_vector", + "aord_vector is a null pointer"); + } + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + int32_t mask = 1 << 3; + + if ( (ctx->jastrow.uninitialized & mask) != 0) { + return QMCKL_NOT_PROVIDED; + } + + assert (ctx->jastrow.aord_vector != NULL); + memcpy(aord_vector, ctx->jastrow.aord_vector, ctx->jastrow.aord_num*sizeof(double)); + return QMCKL_SUCCESS; +} + +qmckl_exit_code qmckl_get_jastrow_bord_vector (const qmckl_context context, double * const bord_vector) { + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return (char) 0; + } + + if (bord_vector == NULL) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_get_jastrow_bord_vector", + "bord_vector is a null pointer"); + } + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + int32_t mask = 1 << 4; + + if ( (ctx->jastrow.uninitialized & mask) != 0) { + return QMCKL_NOT_PROVIDED; + } + + assert (ctx->jastrow.bord_vector != NULL); + memcpy(bord_vector, ctx->jastrow.bord_vector, ctx->jastrow.bord_num*sizeof(double)); + return QMCKL_SUCCESS; +} + +qmckl_exit_code qmckl_get_jastrow_cord_vector (const qmckl_context context, double * const cord_vector) { + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return (char) 0; + } + + if (cord_vector == NULL) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_get_jastrow_cord_vector", + "cord_vector is a null pointer"); + } + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + int32_t mask = 1 << 5; + + if ( (ctx->jastrow.uninitialized & mask) != 0) { + return QMCKL_NOT_PROVIDED; + } + + assert (ctx->jastrow.cord_vector != NULL); + memcpy(cord_vector, ctx->jastrow.cord_vector, ctx->jastrow.cord_num*sizeof(double)); + return QMCKL_SUCCESS; +} + + #+end_src + +** Initialization functions + + To prepare for the Jastrow and its derivative, all the following functions need to be + 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_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 + #+begin_src c :exports none +if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + +qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + #+end_src + + #+NAME:post2 + #+begin_src c :exports none +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_; + } + +return QMCKL_SUCCESS; + #+end_src + + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +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) { +<> + + if (aord_num <= 0) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_set_jastrow_ord_num", + "aord_num <= 0"); + } + + if (bord_num <= 0) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_set_jastrow_ord_num", + "bord_num <= 0"); + } + + if (cord_num <= 0) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_set_jastrow_ord_num", + "cord_num <= 0"); + } + + int32_t mask = 1 << 0; + ctx->jastrow.aord_num = aord_num; + ctx->jastrow.bord_num = bord_num; + ctx->jastrow.cord_num = cord_num; + + <> +} + +qmckl_exit_code qmckl_set_jastrow_type_nucl_num(qmckl_context context, const int64_t type_nucl_num) { +<> + + if (type_nucl_num <= 0) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_set_jastrow_type_nucl_num", + "type_nucl_num < 0"); + } + + int32_t mask = 1 << 1; + ctx->jastrow.type_nucl_num = type_nucl_num; + + <> +} + +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; + + 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, + "qmckl_set_jastrow_coefficient", + "aord_num is not set"); + } + + if (aord_vector == NULL) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "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) { + return qmckl_failwith( context, rc, + "qmckl_set_ord_vector", + NULL); + } + } + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = (aord_num + 1) * type_nucl_num * sizeof(double); + double* new_array = (double*) qmckl_malloc(context, mem_info); + + if(new_array == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_set_jastrow_coefficient", + NULL); + } + + memcpy(new_array, aord_vector, mem_info.size); + + ctx->jastrow.aord_vector = new_array; + + <> +} + +qmckl_exit_code qmckl_set_jastrow_bord_vector(qmckl_context context, double const * bord_vector) { +<> + + int32_t mask = 1 << 4; + + int64_t bord_num; + qmckl_exit_code rc = qmckl_get_jastrow_bord_num(context, &bord_num); + if (rc != QMCKL_SUCCESS) return rc; + + if (bord_num == 0) { + return qmckl_failwith( context, + QMCKL_FAILURE, + "qmckl_set_jastrow_coefficient", + "bord_num is not set"); + } + + if (bord_vector == NULL) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "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) { + return qmckl_failwith( context, rc, + "qmckl_set_ord_vector", + NULL); + } + } + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = (bord_num + 1) * sizeof(double); + double* new_array = (double*) qmckl_malloc(context, mem_info); + + if(new_array == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_set_jastrow_coefficient", + NULL); + } + + memcpy(new_array, bord_vector, mem_info.size); + + ctx->jastrow.bord_vector = new_array; + + <> +} + +qmckl_exit_code qmckl_set_jastrow_cord_vector(qmckl_context context, double const * cord_vector) { +<> + + int32_t mask = 1 << 5; + + int64_t cord_num; + qmckl_exit_code rc = qmckl_get_jastrow_cord_num(context, &cord_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 (cord_num == 0) { + return qmckl_failwith( context, + QMCKL_FAILURE, + "qmckl_set_jastrow_coefficient", + "cord_num is not set"); + } + + if (cord_vector == NULL) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "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) { + return qmckl_failwith( context, rc, + "qmckl_set_ord_vector", + NULL); + } + } + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = cord_num * type_nucl_num * sizeof(double); + double* new_array = (double*) qmckl_malloc(context, mem_info); + + if(new_array == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_set_jastrow_coefficient", + NULL); + } + + memcpy(new_array, cord_vector, mem_info.size); + + ctx->jastrow.cord_vector = new_array; + + <> +} + +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 + computed to accelerate the calculations. The intermediates factors + are precontracted using BLAS LEVEL 3 operations for an optimal FLOP count. + + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none +qmckl_exit_code qmckl_finalize_jastrow(qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_finalize_jastrow(qmckl_context context) { + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_INVALID_CONTEXT; + } + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + /* ----------------------------------- */ + /* Check for the necessary information */ + /* ----------------------------------- */ + + /* Check for the electron data + 1. elec_num + 2. ee_distances_rescaled + */ + if (!(ctx->electron.provided)) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_electron", + NULL); + } + + /* Check for the nucleus data + 1. nucl_num + 2. en_distances_rescaled + */ + if (!(ctx->nucleus.provided)) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_nucleus", + NULL); + } + + qmckl_exit_code rc = QMCKL_FAILURE; + + + /* ----------------------------------- */ + /* Start calculation of data */ + /* ----------------------------------- */ + + +} + #+end_src + +** Test + #+begin_src c :tangle (eval c_test) +/* Reference input data */ +int64_t walk_num = n2_walk_num; +int64_t elec_num = n2_elec_num; +int64_t elec_up_num = n2_elec_up_num; +int64_t elec_dn_num = n2_elec_dn_num; +double rescale_factor_kappa_ee = 1.0; +double rescale_factor_kappa_en = 1.0; +double nucl_rescale_factor_kappa = 1.0; +double* elec_coord = &(n2_elec_coord[0][0][0]); + +const double* nucl_charge = n2_charge; +int64_t nucl_num = n2_nucl_num; +double* charge = n2_charge; +double* nucl_coord = &(n2_nucl_coord[0][0]); + +/* Provide Electron data */ + +qmckl_exit_code rc; + +assert(!qmckl_electron_provided(context)); + +int64_t n; +rc = qmckl_get_electron_num (context, &n); +assert(rc == QMCKL_NOT_PROVIDED); + +rc = qmckl_get_electron_up_num (context, &n); +assert(rc == QMCKL_NOT_PROVIDED); + +rc = qmckl_get_electron_down_num (context, &n); +assert(rc == QMCKL_NOT_PROVIDED); + + +rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num); +assert(rc == QMCKL_SUCCESS); +assert(!qmckl_electron_provided(context)); + +rc = qmckl_get_electron_up_num (context, &n); +assert(rc == QMCKL_SUCCESS); +assert(n == elec_up_num); + +rc = qmckl_get_electron_down_num (context, &n); +assert(rc == QMCKL_SUCCESS); +assert(n == elec_dn_num); + +rc = qmckl_get_electron_num (context, &n); +assert(rc == QMCKL_SUCCESS); +assert(n == elec_num); + +double k_ee = 0.; +double k_en = 0.; +rc = qmckl_get_electron_rescale_factor_ee (context, &k_ee); +assert(rc == QMCKL_SUCCESS); +assert(k_ee == 1.0); + +rc = qmckl_get_electron_rescale_factor_en (context, &k_en); +assert(rc == QMCKL_SUCCESS); +assert(k_en == 1.0); + +rc = qmckl_set_electron_rescale_factor_en(context, rescale_factor_kappa_en); +assert(rc == QMCKL_SUCCESS); + +rc = qmckl_set_electron_rescale_factor_ee(context, rescale_factor_kappa_ee); +assert(rc == QMCKL_SUCCESS); + +rc = qmckl_get_electron_rescale_factor_ee (context, &k_ee); +assert(rc == QMCKL_SUCCESS); +assert(k_ee == rescale_factor_kappa_ee); + +rc = qmckl_get_electron_rescale_factor_en (context, &k_en); +assert(rc == QMCKL_SUCCESS); +assert(k_en == rescale_factor_kappa_en); + + +int64_t w; +rc = qmckl_get_electron_walk_num (context, &w); +assert(rc == QMCKL_NOT_PROVIDED); + + +rc = qmckl_set_electron_walk_num (context, walk_num); +assert(rc == QMCKL_SUCCESS); + +rc = qmckl_get_electron_walk_num (context, &w); +assert(rc == QMCKL_SUCCESS); +assert(w == walk_num); + +assert(qmckl_electron_provided(context)); + +rc = qmckl_set_electron_coord (context, 'N', elec_coord); +assert(rc == QMCKL_SUCCESS); + +double elec_coord2[walk_num*3*elec_num]; + +rc = qmckl_get_electron_coord (context, 'N', elec_coord2); +assert(rc == QMCKL_SUCCESS); +for (int64_t i=0 ; i<3*elec_num ; ++i) { + assert( elec_coord[i] == elec_coord2[i] ); + } + + +/* Provide Nucleus data */ + +assert(!qmckl_nucleus_provided(context)); + +rc = qmckl_get_nucleus_num (context, &n); +assert(rc == QMCKL_NOT_PROVIDED); + + +rc = qmckl_set_nucleus_num (context, nucl_num); +assert(rc == QMCKL_SUCCESS); +assert(!qmckl_nucleus_provided(context)); + +rc = qmckl_get_nucleus_num (context, &n); +assert(rc == QMCKL_SUCCESS); +assert(n == nucl_num); + +double k; +rc = qmckl_get_nucleus_rescale_factor (context, &k); +assert(rc == QMCKL_SUCCESS); +assert(k == 1.0); + + +rc = qmckl_set_nucleus_rescale_factor (context, nucl_rescale_factor_kappa); +assert(rc == QMCKL_SUCCESS); + +rc = qmckl_get_nucleus_rescale_factor (context, &k); +assert(rc == QMCKL_SUCCESS); +assert(k == nucl_rescale_factor_kappa); + +double nucl_coord2[3*nucl_num]; + +rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2); +assert(rc == QMCKL_NOT_PROVIDED); + +rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0])); +assert(rc == QMCKL_SUCCESS); + +assert(!qmckl_nucleus_provided(context)); + +rc = qmckl_get_nucleus_coord (context, 'N', nucl_coord2); +assert(rc == QMCKL_SUCCESS); +for (int64_t k=0 ; k<3 ; ++k) { + for (int64_t i=0 ; ijastrow.asymp_jasb, 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_asymp_jasb(qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_provide_asymp_jasb(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 kappa is provided */ + double rescale_factor_kappa_ee; + rc = qmckl_get_electron_rescale_factor_ee(context, &rescale_factor_kappa_ee); + if(rc != QMCKL_SUCCESS) return rc; + + /* Compute if necessary */ + if (ctx->date > ctx->jastrow.asymp_jasb_date) { + + /* Allocate array */ + if (ctx->jastrow.asymp_jasb == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = 2 * sizeof(double); + double* asymp_jasb = (double*) qmckl_malloc(context, mem_info); + + if (asymp_jasb == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_asymp_jasb", + NULL); + } + ctx->jastrow.asymp_jasb = asymp_jasb; + } + + qmckl_exit_code rc = + qmckl_compute_asymp_jasb(context, + ctx->jastrow.bord_num, + ctx->jastrow.bord_vector, + rescale_factor_kappa_ee, + ctx->jastrow.asymp_jasb); + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->jastrow.asymp_jasb_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + #+end_src + +*** Compute + :PROPERTIES: + :Name: qmckl_compute_asymp_jasb + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_asymp_jasb_args + | qmckl_context | context | in | Global state | + | int64_t | bord_num | in | Number of electrons | + | double | bord_vector[bord_num + 1] | in | Number of walkers | + | double | rescale_factor_kappa_ee | in | Electron coordinates | + | double | asymp_jasb[2] | out | Electron-electron distances | + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer function qmckl_compute_asymp_jasb_f(context, bord_num, bord_vector, rescale_factor_kappa_ee, asymp_jasb) & + result(info) + use qmckl + implicit none + integer(qmckl_context), intent(in) :: context + integer*8 , intent(in) :: bord_num + double precision , intent(in) :: bord_vector(bord_num + 1) + double precision , intent(in) :: rescale_factor_kappa_ee + double precision , intent(out) :: asymp_jasb(2) + + integer*8 :: i, p + double precision :: kappa_inv, x, asym_one + kappa_inv = 1.0d0 / rescale_factor_kappa_ee + + info = QMCKL_SUCCESS + + if (context == QMCKL_NULL_CONTEXT) then + info = QMCKL_INVALID_CONTEXT + return + endif + + if (bord_num <= 0) then + info = QMCKL_INVALID_ARG_2 + return + endif + + asym_one = bord_vector(1) * kappa_inv / (1.0d0 + bord_vector(2) * kappa_inv) + asymp_jasb(:) = (/asym_one, 0.5d0 * asym_one/) + + do i = 1, 2 + x = kappa_inv + do p = 2, bord_num + x = x * kappa_inv + asymp_jasb(i) = asymp_jasb(i) + bord_vector(p + 1) * x + end do + end do + +end function qmckl_compute_asymp_jasb_f + #+end_src + + #+CALL: generate_c_header(table=qmckl_asymp_jasb_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + + #+RESULTS: + #+begin_src c :tangle (eval h_private_func) :comments org + qmckl_exit_code qmckl_compute_asymp_jasb ( + const qmckl_context context, + const int64_t bord_num, + const double* bord_vector, + const double rescale_factor_kappa_ee, + double* const asymp_jasb ); + #+end_src + + + #+CALL: generate_c_interface(table=qmckl_asymp_jasb_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_asymp_jasb & + (context, bord_num, bord_vector, rescale_factor_kappa_ee, asymp_jasb) & + 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 :: bord_num + real (c_double ) , intent(in) :: bord_vector(bord_num + 1) + real (c_double ) , intent(in) , value :: rescale_factor_kappa_ee + real (c_double ) , intent(out) :: asymp_jasb(2) + + integer(c_int32_t), external :: qmckl_compute_asymp_jasb_f + info = qmckl_compute_asymp_jasb_f & + (context, bord_num, bord_vector, rescale_factor_kappa_ee, asymp_jasb) + + end function qmckl_compute_asymp_jasb + #+end_src + +*** Test + #+name: asymp_jasb + #+begin_src python :results output :exports none :noweb yes +import numpy as np + +<> + +asym_one = bord_vector[0] * kappa_inv / (1.0 + bord_vector[1]*kappa_inv) +asymp_jasb = np.array([asym_one, 0.5 * asym_one]) + +for i in range(2): + x = kappa_inv + for p in range(1,bord_num): + x = x * kappa_inv + asymp_jasb[i] += bord_vector[p + 1] * x + +print("asym_one : ", asym_one) +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 + : asymp_jasb[1] : 0.31567342786262853 + + #+begin_src c :tangle (eval c_test) +assert(qmckl_electron_provided(context)); + +int64_t type_nucl_num = n2_type_nucl_num; +int64_t* type_nucl_vector = &(n2_type_nucl_vector[0]); +int64_t aord_num = n2_aord_num; +int64_t bord_num = n2_bord_num; +int64_t cord_num = n2_cord_num; +double* aord_vector = &(n2_aord_vector[0][0]); +double* bord_vector = &(n2_bord_vector[0]); +double* cord_vector = &(n2_cord_vector[0][0]); + +/* Initialize the Jastrow data */ +rc = qmckl_init_jastrow(context); +assert(!qmckl_jastrow_provided(context)); + +/* Set the data */ +rc = qmckl_set_jastrow_ord_num(context, aord_num, bord_num, cord_num); +assert(rc == QMCKL_SUCCESS); +rc = qmckl_set_jastrow_type_nucl_num(context, type_nucl_num); +assert(rc == QMCKL_SUCCESS); +rc = qmckl_set_jastrow_type_nucl_vector(context, 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); +assert(rc == QMCKL_SUCCESS); +rc = qmckl_set_jastrow_cord_vector(context, cord_vector); +assert(rc == QMCKL_SUCCESS); +rc = qmckl_set_jastrow_dependencies(context); +assert(rc == QMCKL_SUCCESS); + +/* Check if Jastrow is properly initialized */ +assert(qmckl_jastrow_provided(context)); + +double asymp_jasb[2]; +rc = qmckl_get_jastrow_asymp_jasb(context, asymp_jasb); + +// calculate asymp_jasb +assert(fabs(asymp_jasb[0]-0.5323750557252571) < 1.e-12); +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~. + + \[ +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 + 1) + 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 + +** 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\) + 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); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_get_jastrow_factor_ee_deriv_e(qmckl_context context, double* const factor_ee_deriv_e) +{ + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + + qmckl_exit_code rc; + + rc = qmckl_provide_factor_ee_deriv_e(context); + if (rc != QMCKL_SUCCESS) return rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + int64_t sze = ctx->electron.walk_num * 4 * ctx->electron.num; + memcpy(factor_ee_deriv_e, ctx->jastrow.factor_ee_deriv_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_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) +{ + + 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; + + /* Check if ee rescaled distance deriv e is provided */ + rc = qmckl_provide_ee_distance_rescaled_deriv_e(context); + if(rc != QMCKL_SUCCESS) return rc; + + /* Compute if necessary */ + if (ctx->date > ctx->jastrow.factor_ee_deriv_e_date) { + + /* Allocate array */ + if (ctx->jastrow.factor_ee_deriv_e == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->electron.walk_num * 4 * ctx->electron.num * sizeof(double); + double* factor_ee_deriv_e = (double*) qmckl_malloc(context, mem_info); + + if (factor_ee_deriv_e == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_provide_factor_ee_deriv_e", + NULL); + } + 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, + ctx->electron.num, + ctx->electron.up_num, + ctx->jastrow.bord_num, + ctx->jastrow.bord_vector, + ctx->electron.ee_distance_rescaled, + ctx->electron.ee_distance_rescaled_deriv_e, + ctx->jastrow.asymp_jasb, + ctx->jastrow.factor_ee_deriv_e); + 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_deriv_e + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_factor_ee_deriv_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 | 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 | ee_distance_rescaled_deriv_e[walk_num][4][elec_num][elec_num] | in | Electron-electron distances | + | double | asymp_jasb[2] | in | Electron-electron distances | + | double | factor_ee_deriv_e[walk_num][4][elec_num] | out | Electron-electron distances | + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer function qmckl_compute_factor_ee_deriv_e_f(context, walk_num, elec_num, up_num, bord_num, & + bord_vector, ee_distance_rescaled, ee_distance_rescaled_deriv_e, & + asymp_jasb, factor_ee_deriv_e) & + 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 + 1) + double precision , intent(in) :: ee_distance_rescaled(walk_num, elec_num, elec_num) + double precision , intent(in) :: ee_distance_rescaled_deriv_e(walk_num, 4, elec_num, elec_num) + double precision , intent(in) :: asymp_jasb(2) + double precision , intent(out) :: factor_ee_deriv_e(elec_num,4,walk_num) + + integer*8 :: i, j, p, ipar, nw, ii + double precision :: x, spin_fact, y + double precision :: den, invden, invden2, invden3, xinv + double precision :: lap1, lap2, lap3, third + double precision, dimension(3) :: pow_ser_g + double precision, dimension(4) :: dx + + 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_deriv_e = 0.0d0 + third = 1.0d0 / 3.0d0 + + do nw =1, walk_num + do j = 1, elec_num + do i = 1, elec_num + x = ee_distance_rescaled(nw, i, j) + if(abs(x) < 1.0d-18) cycle + pow_ser_g = 0.0d0 + spin_fact = 1.0d0 + den = 1.0d0 + bord_vector(2) * x + invden = 1.0d0 / den + invden2 = invden * invden + invden3 = invden2 * invden + xinv = 1.0d0 / (x + 1.0d-18) + ipar = 1 + + do ii = 1, 4 + dx(ii) = ee_distance_rescaled_deriv_e(nw, ii, i, j) + end do + + if((i .LE. up_num .AND. j .LE. up_num ) .OR. & + (i .GT. up_num .AND. j .GT. up_num)) then + spin_fact = 0.5d0 + endif + + lap1 = 0.0d0 + lap2 = 0.0d0 + lap3 = 0.0d0 + do ii = 1, 3 + x = ee_distance_rescaled(nw, i, j) + if(abs(x) < 1.0d-18) cycle + do p = 2, bord_num + y = p * bord_vector(p + 1) * x + pow_ser_g(ii) = pow_ser_g(ii) + y * dx(ii) + lap1 = lap1 + (p - 1) * y * xinv * dx(ii) * dx(ii) + 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) * & + dx(ii) * invden2 + pow_ser_g(ii) + end do + + ii = 4 + lap2 = lap2 * dx(ii) * third + lap3 = lap3 + den * dx(ii) + lap3 = lap3 * (spin_fact * bord_vector(1) * invden3) + factor_ee_deriv_e( j, ii, nw) = factor_ee_deriv_e( j, ii, nw) + lap1 + lap2 + lap3 + + end do + end do + end do + +end function qmckl_compute_factor_ee_deriv_e_f + #+end_src + + #+CALL: generate_c_header(table=qmckl_factor_ee_deriv_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_factor_ee_deriv_e ( + 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* ee_distance_rescaled_deriv_e, + const double* asymp_jasb, + double* const factor_ee_deriv_e ); + #+end_src + + + #+CALL: generate_c_interface(table=qmckl_factor_ee_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_ee_deriv_e & + (context, & + walk_num, & + elec_num, & + up_num, & + bord_num, & + bord_vector, & + ee_distance_rescaled, & + ee_distance_rescaled_deriv_e, & + asymp_jasb, & + factor_ee_deriv_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 :: 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) :: ee_distance_rescaled_deriv_e(elec_num,elec_num,4,walk_num) + real (c_double ) , intent(in) :: asymp_jasb(2) + real (c_double ) , intent(out) :: factor_ee_deriv_e(elec_num,4,walk_num) + + integer(c_int32_t), external :: qmckl_compute_factor_ee_deriv_e_f + info = qmckl_compute_factor_ee_deriv_e_f & + (context, & + walk_num, & + elec_num, & + up_num, & + bord_num, & + bord_vector, & + ee_distance_rescaled, & + ee_distance_rescaled_deriv_e, & + asymp_jasb, & + factor_ee_deriv_e) + + end function qmckl_compute_factor_ee_deriv_e + #+end_src + +*** Test + #+begin_src python :results output :exports none :noweb yes +import numpy as np + +<> + +<> + +kappa = 1.0 + +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]) + +elec_dist_deriv_e = np.zeros(shape=(4,elec_num, elec_num),dtype=float) +for j in range(elec_num): + for i in range(elec_num): + rij_inv = 1.0 / elec_dist[i, j] + for ii in range(3): + elec_dist_deriv_e[ii, i, j] = (elec_coord[j][ii] - elec_coord[i][ii]) * rij_inv + elec_dist_deriv_e[3, i, j] = 2.0 * rij_inv + elec_dist_deriv_e[:, j, j] = 0.0 + +ee_distance_rescaled_deriv_e = np.zeros(shape=(4,elec_num,elec_num),dtype=float) +for j in range(elec_num): + for i in range(elec_num): + f = 1.0 - kappa * ee_distance_rescaled[i][j] + for ii in range(4): + ee_distance_rescaled_deriv_e[ii][i][j] = elec_dist_deriv_e[ii][i][j] + ee_distance_rescaled_deriv_e[3][i][j] = ee_distance_rescaled_deriv_e[3][i][j] + \ + (-kappa * ee_distance_rescaled_deriv_e[0][i][j] * ee_distance_rescaled_deriv_e[0][i][j]) + \ + (-kappa * ee_distance_rescaled_deriv_e[1][i][j] * ee_distance_rescaled_deriv_e[1][i][j]) + \ + (-kappa * ee_distance_rescaled_deriv_e[2][i][j] * ee_distance_rescaled_deriv_e[2][i][j]) + for ii in range(4): + ee_distance_rescaled_deriv_e[ii][i][j] = ee_distance_rescaled_deriv_e[ii][i][j] * f + +third = 1.0 / 3.0 +factor_ee_deriv_e = np.zeros(shape=(4,elec_num),dtype=float) +dx = np.zeros(shape=(4),dtype=float) +pow_ser_g = np.zeros(shape=(4),dtype=float) +for j in range(elec_num): + for i in range(elec_num): + x = ee_distance_rescaled[j][i] + if abs(x) < 1e-18: + continue + pow_ser_g = np.zeros(shape=(4),dtype=float) + spin_fact = 1.0 + den = 1.0 + bord_vector[1] * ee_distance_rescaled[j][i] + invden = 1.0 / den + invden2 = invden * invden + invden3 = invden2 * invden + xinv = 1.0 / (ee_distance_rescaled[j][i] + 1.0E-18) + ipar = 1 + + for ii in range(4): + dx[ii] = ee_distance_rescaled_deriv_e[ii][j][i] + + 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 + for ii in range(3): + x = ee_distance_rescaled[j][i] + if x < 1e-18: + continue + for p in range(2,bord_num+1): + y = p * bord_vector[(p-1) + 1] * x + pow_ser_g[ii] = pow_ser_g[ii] + y * dx[ii] + 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] * \ + dx[ii] * invden2 + pow_ser_g[ii] + + ii = 3 + lap2 = lap2 * dx[ii] * third + lap3 = lap3 + den * dx[ii] + lap3 = lap3 * (spin_fact * bord_vector[0] * invden3) + factor_ee_deriv_e[ii][j] = factor_ee_deriv_e[ii][j] + lap1 + lap2 + lap3 + +print("factor_ee_deriv_e[0][0]:",factor_ee_deriv_e[0][0]) +print("factor_ee_deriv_e[1][0]:",factor_ee_deriv_e[1][0]) +print("factor_ee_deriv_e[2][0]:",factor_ee_deriv_e[2][0]) +print("factor_ee_deriv_e[3][0]:",factor_ee_deriv_e[3][0]) +print(factor_ee_deriv_e) + + #+end_src + + #+RESULTS: + #+begin_example + asym_one : 0.43340325572525706 + asymp_jasb[0] : 0.5323750557252571 + asymp_jasb[1] : 0.31567342786262853 + factor_ee_deriv_e[0][0]: 0.16364894652107934 + factor_ee_deriv_e[1][0]: -0.6927548119830084 + factor_ee_deriv_e[2][0]: 0.073267755223968 + factor_ee_deriv_e[3][0]: 1.5111672803213185 + [[ 0.16364895 0.60354957 -0.19825547 0.02359797 -0.13123153 -0.18789233 + 0.07762515 -0.42459184 0.27920265 -0.2056531 ] + [-0.69275481 0.15690393 0.09831069 0.18490587 0.04361723 0.3250686 + 0.12657961 -0.01736522 -0.40149005 0.17622416] + [ 0.07326776 -0.27532276 0.22396943 0.18771633 -0.34506246 0.07298062 + 0.63302352 -0.00910198 -0.30238713 -0.25908332] + [ 1.51116728 1.5033247 0.00325003 2.89377255 0.1338393 2.15893795 + 1.74732003 0.23561147 2.67455607 0.82810434]] + #+end_example + + + #+begin_src c :tangle (eval c_test) +/* Check if Jastrow is properly initialized */ +assert(qmckl_jastrow_provided(context)); + +// calculate factor_ee_deriv_e +double factor_ee_deriv_e[walk_num][4][elec_num]; +rc = qmckl_get_jastrow_factor_ee_deriv_e(context, &(factor_ee_deriv_e[0][0][0])); + +// check factor_ee_deriv_e +assert(fabs(factor_ee_deriv_e[0][0][0]-0.16364894652107934) < 1.e-12); +assert(fabs(factor_ee_deriv_e[0][1][0]+0.6927548119830084 ) < 1.e-12); +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 + +** 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[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(nucl_num) + double precision , intent(in) :: aord_vector(aord_num + 1, type_nucl_num) + double precision , intent(in) :: en_distance_rescaled(walk_num, nucl_num, elec_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, a, i) + power_ser = 0.0d0 + + do p = 2, aord_num + x = x * en_distance_rescaled(nw, a, i) + 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, a, i) / & + (1.0d0 + aord_vector(2, type_nucl_vector(a)) * & + en_distance_rescaled(nw, a, i)) & + + 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(nucl_num) + integer (c_int64_t) , intent(in) , value :: aord_num + real (c_double ) , intent(in) :: aord_vector(aord_num + 1, type_nucl_num) + real (c_double ) , intent(in) :: en_distance_rescaled(walk_num, nucl_num, elec_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 + +** Electron-nucleus component derivative \(f'_{en}\) + Calculate the electron-electron jastrow component ~factor_en_deriv_e~ derivative + with respect to the electron coordinates using the ~en_distance_rescaled~ and + ~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); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_get_jastrow_factor_en_deriv_e(qmckl_context context, double* const factor_en_deriv_e) +{ + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + + qmckl_exit_code rc; + + rc = qmckl_provide_factor_en_deriv_e(context); + if (rc != QMCKL_SUCCESS) return rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + int64_t sze = ctx->electron.walk_num * 4 * ctx->electron.num; + memcpy(factor_en_deriv_e, ctx->jastrow.factor_en_deriv_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_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) +{ + + 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; + + /* Check if en rescaled distance derivatives is provided */ + rc = qmckl_provide_en_distance_rescaled_deriv_e(context); + if(rc != QMCKL_SUCCESS) return rc; + + /* Compute if necessary */ + if (ctx->date > ctx->jastrow.factor_en_deriv_e_date) { + + /* Allocate array */ + if (ctx->jastrow.factor_en_deriv_e == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->electron.walk_num * 4 * ctx->electron.num * sizeof(double); + double* factor_en_deriv_e = (double*) qmckl_malloc(context, mem_info); + + if (factor_en_deriv_e == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_provide_factor_en_deriv_e", + NULL); + } + 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, + 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->electron.en_distance_rescaled_deriv_e, + ctx->jastrow.factor_en_deriv_e); + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->jastrow.factor_en_deriv_e_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + #+end_src + +*** Compute + :PROPERTIES: + :Name: qmckl_compute_factor_en_deriv_e + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_factor_en_deriv_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 | nucl_num | in | Number of nucleii | + | int64_t | type_nucl_num | in | Number of unique nuclei | + | int64_t | type_nucl_vector[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 | en_distance_rescaled_deriv_e[walk_num][4][nucl_num][elec_num] | in | Electron-nucleus distance derivatives | + | double | factor_en_deriv_e[walk_num][4][elec_num] | out | Electron-nucleus jastrow | + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer function qmckl_compute_factor_en_deriv_e_f(context, walk_num, elec_num, nucl_num, type_nucl_num, & + type_nucl_vector, aord_num, aord_vector, & + en_distance_rescaled, en_distance_rescaled_deriv_e, factor_en_deriv_e) & + 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(nucl_num) + double precision , intent(in) :: aord_vector(aord_num + 1, type_nucl_num) + double precision , intent(in) :: en_distance_rescaled(walk_num, elec_num, nucl_num) + double precision , intent(in) :: en_distance_rescaled_deriv_e(walk_num, 4, elec_num, nucl_num) + double precision , intent(out) :: factor_en_deriv_e(elec_num,4,walk_num) + + integer*8 :: i, a, p, ipar, nw, ii + double precision :: x, spin_fact, den, invden, invden2, invden3, xinv + double precision :: y, lap1, lap2, lap3, third + double precision, dimension(3) :: power_ser_g + double precision, dimension(4) :: dx + + 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_deriv_e = 0.0d0 + third = 1.0d0 / 3.0d0 + + do nw =1, walk_num + do a = 1, nucl_num + do i = 1, elec_num + x = en_distance_rescaled(nw, i, a) + if(abs(x) < 1.0d-18) continue + power_ser_g = 0.0d0 + den = 1.0d0 + aord_vector(2, type_nucl_vector(a)) * x + invden = 1.0d0 / den + invden2 = invden * invden + invden3 = invden2 * invden + xinv = 1.0d0 / x + + do ii = 1, 4 + dx(ii) = en_distance_rescaled_deriv_e(nw, ii, i, a) + end do + + lap1 = 0.0d0 + lap2 = 0.0d0 + lap3 = 0.0d0 + do ii = 1, 3 + x = en_distance_rescaled(nw, i, a) + do p = 2, aord_num + y = p * aord_vector(p + 1, type_nucl_vector(a)) * x + power_ser_g(ii) = power_ser_g(ii) + y * dx(ii) + lap1 = lap1 + (p - 1) * y * xinv * dx(ii) * dx(ii) + lap2 = lap2 + y + x = x * en_distance_rescaled(nw, i, a) + end do + + lap3 = lap3 - 2.0d0 * aord_vector(2, type_nucl_vector(a)) * dx(ii) * dx(ii) + + factor_en_deriv_e(i, ii, nw) = factor_en_deriv_e(i, ii, nw) + aord_vector(1, type_nucl_vector(a)) & + * dx(ii) * invden2 & + + power_ser_g(ii) + + end do + + ii = 4 + lap2 = lap2 * dx(ii) * third + lap3 = lap3 + den * dx(ii) + lap3 = lap3 * aord_vector(1, type_nucl_vector(a)) * invden3 + factor_en_deriv_e(i, ii, nw) = factor_en_deriv_e(i, ii, nw) + lap1 + lap2 + lap3 + + end do + end do + end do + +end function qmckl_compute_factor_en_deriv_e_f + #+end_src + + #+CALL: generate_c_header(table=qmckl_factor_en_deriv_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_factor_en_deriv_e ( + 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, + const double* en_distance_rescaled_deriv_e, + double* const factor_en_deriv_e ); + #+end_src + + + #+CALL: generate_c_interface(table=qmckl_factor_en_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_en_deriv_e & + (context, & + walk_num, & + elec_num, & + nucl_num, & + type_nucl_num, & + type_nucl_vector, & + aord_num, & + aord_vector, & + en_distance_rescaled, & + en_distance_rescaled_deriv_e, & + factor_en_deriv_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 :: nucl_num + integer (c_int64_t) , intent(in) , value :: type_nucl_num + integer (c_int64_t) , intent(in) :: type_nucl_vector(nucl_num) + integer (c_int64_t) , intent(in) , value :: aord_num + real (c_double ) , intent(in) :: aord_vector(aord_num + 1, type_nucl_num) + real (c_double ) , intent(in) :: en_distance_rescaled(elec_num,nucl_num,walk_num) + real (c_double ) , intent(in) :: en_distance_rescaled_deriv_e(elec_num,nucl_num,4,walk_num) + real (c_double ) , intent(out) :: factor_en_deriv_e(elec_num,4,walk_num) + + integer(c_int32_t), external :: qmckl_compute_factor_en_deriv_e_f + info = qmckl_compute_factor_en_deriv_e_f & + (context, & + walk_num, & + elec_num, & + nucl_num, & + type_nucl_num, & + type_nucl_vector, & + aord_num, & + aord_vector, & + en_distance_rescaled, & + en_distance_rescaled_deriv_e, & + factor_en_deriv_e) + + end function qmckl_compute_factor_en_deriv_e + #+end_src + +*** Test + #+begin_src python :results output :exports none :noweb yes +import numpy as np + +<> + +kappa = 1.0 + +elec_coord = np.array(elec_coord)[0] +nucl_coord = np.array(nucl_coord) +elnuc_dist = np.zeros(shape=(elec_num, nucl_num),dtype=float) +for i in range(elec_num): + for j in range(nucl_num): + elnuc_dist[i, j] = np.linalg.norm(elec_coord[i] - nucl_coord[:,j]) + +elnuc_dist_deriv_e = np.zeros(shape=(4, elec_num, nucl_num),dtype=float) +for a in range(nucl_num): + for i in range(elec_num): + rij_inv = 1.0 / elnuc_dist[i, a] + for ii in range(3): + elnuc_dist_deriv_e[ii, i, a] = (elec_coord[i][ii] - nucl_coord[ii][a]) * rij_inv + elnuc_dist_deriv_e[3, i, a] = 2.0 * rij_inv + +en_distance_rescaled_deriv_e = np.zeros(shape=(4,elec_num,nucl_num),dtype=float) +for a in range(nucl_num): + for i in range(elec_num): + f = 1.0 - kappa * en_distance_rescaled[i][a] + for ii in range(4): + en_distance_rescaled_deriv_e[ii][i][a] = elnuc_dist_deriv_e[ii][i][a] + en_distance_rescaled_deriv_e[3][i][a] = en_distance_rescaled_deriv_e[3][i][a] + \ + (-kappa * en_distance_rescaled_deriv_e[0][i][a] * en_distance_rescaled_deriv_e[0][i][a]) + \ + (-kappa * en_distance_rescaled_deriv_e[1][i][a] * en_distance_rescaled_deriv_e[1][i][a]) + \ + (-kappa * en_distance_rescaled_deriv_e[2][i][a] * en_distance_rescaled_deriv_e[2][i][a]) + for ii in range(4): + en_distance_rescaled_deriv_e[ii][i][a] = en_distance_rescaled_deriv_e[ii][i][a] * f + +third = 1.0 / 3.0 +factor_en_deriv_e = np.zeros(shape=(4,elec_num),dtype=float) +dx = np.zeros(shape=(4),dtype=float) +pow_ser_g = np.zeros(shape=(3),dtype=float) +for a in range(nucl_num): + for i in range(elec_num): + x = en_distance_rescaled[i][a] + if abs(x) < 1e-18: + continue + pow_ser_g = np.zeros(shape=(3),dtype=float) + den = 1.0 + aord_vector[1][type_nucl_vector[a]-1] * x + invden = 1.0 / den + invden2 = invden * invden + invden3 = invden2 * invden + xinv = 1.0 / (x + 1.0E-18) + + for ii in range(4): + dx[ii] = en_distance_rescaled_deriv_e[ii][i][a] + + lap1 = 0.0 + lap2 = 0.0 + lap3 = 0.0 + for ii in range(3): + x = en_distance_rescaled[i][a] + if x < 1e-18: + continue + for p in range(2,aord_num+1): + y = p * aord_vector[(p-1) + 1][type_nucl_vector[a]-1] * x + pow_ser_g[ii] = pow_ser_g[ii] + y * dx[ii] + 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] * \ + dx[ii] * invden2 + pow_ser_g[ii] + + ii = 3 + lap2 = lap2 * dx[ii] * third + lap3 = lap3 + den * dx[ii] + lap3 = lap3 * (aord_vector[0][type_nucl_vector[a]-1] * invden3) + factor_en_deriv_e[ii][i] = factor_en_deriv_e[ii][i] + lap1 + lap2 + lap3 + +print("factor_en_deriv_e[0][0]:",factor_en_deriv_e[0][0]) +print("factor_en_deriv_e[1][0]:",factor_en_deriv_e[1][0]) +print("factor_en_deriv_e[2][0]:",factor_en_deriv_e[2][0]) +print("factor_en_deriv_e[3][0]:",factor_en_deriv_e[3][0]) + + + #+end_src + + #+RESULTS: + : factor_en_deriv_e[0][0]: 0.11609919541763383 + : 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 */ +assert(qmckl_jastrow_provided(context)); + +// calculate factor_en_deriv_e +double factor_en_deriv_e[walk_num][4][elec_num]; +rc = qmckl_get_jastrow_factor_en_deriv_e(context, &(factor_en_deriv_e[0][0][0])); + +// check factor_en_deriv_e +assert(fabs(factor_en_deriv_e[0][0][0]-0.11609919541763383) < 1.e-12); +assert(fabs(factor_en_deriv_e[0][1][0]+0.23301394780804574) < 1.e-12); +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 + 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 + + ! Prepare table of exponentiated distances raised to appropriate power + een_rescaled_e = 0.0d0 + do nw = 1, walk_num + een_rescaled_e_ij = 0.0d0 + een_rescaled_e_ij(:, 1) = 1.0d0 + + 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, 2) + end do + end do + + ! prepare the actual een table + een_rescaled_e(0, :, :, nw) = 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): + 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): + 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) +assert(fabs(een_rescaled_e[0][0][2][1]-0.08084493981483197) < 1.e-12); +assert(fabs(een_rescaled_e[0][0][3][1]-0.1066745707571846) < 1.e-12); +assert(fabs(een_rescaled_e[0][0][4][1]-0.01754273169464735) < 1.e-12); +assert(fabs(een_rescaled_e[0][1][3][2]-0.02214680362033448) < 1.e-12); +assert(fabs(een_rescaled_e[0][1][4][2]-0.0005700154999202759) < 1.e-12); +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~ stores the table 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 +qmckl_exit_code qmckl_get_jastrow_een_rescaled_e_deriv_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_deriv_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_deriv_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 * 4 * ctx->electron.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1); + memcpy(distance_rescaled, ctx->jastrow.een_rescaled_e_deriv_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_deriv_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_deriv_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_een_rescaled_e(context); + if(rc != QMCKL_SUCCESS) return rc; + + /* Compute if necessary */ + if (ctx->date > ctx->jastrow.een_rescaled_e_deriv_e_date) { + + /* Allocate array */ + if (ctx->jastrow.een_rescaled_e_deriv_e == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->electron.num * 4 * ctx->electron.num * + ctx->electron.walk_num * (ctx->jastrow.cord_num + 1) * sizeof(double); + double* een_rescaled_e_deriv_e = (double*) qmckl_malloc(context, mem_info); + + if (een_rescaled_e_deriv_e == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_een_rescaled_e_deriv_e", + NULL); + } + ctx->jastrow.een_rescaled_e_deriv_e = een_rescaled_e_deriv_e; + } + + qmckl_exit_code rc = + qmckl_compute_factor_een_rescaled_e_deriv_e(context, + ctx->electron.walk_num, + ctx->electron.num, + ctx->jastrow.cord_num, + ctx->electron.rescale_factor_kappa_ee, + ctx->electron.coord_new, + ctx->electron.ee_distance, + ctx->jastrow.een_rescaled_e, + ctx->jastrow.een_rescaled_e_deriv_e); + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->jastrow.een_rescaled_e_deriv_e_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + #+end_src + +*** Compute + :PROPERTIES: + :Name: qmckl_compute_een_rescaled_e_deriv_e + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_factor_een_rescaled_e_deriv_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 | coord_new[walk_num][3][elec_num] | in | Electron coordinates | + | 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] | in | Electron-electron distances | + | double | een_rescaled_e_deriv_e[walk_num][elec_num][4][elec_num][0:cord_num] | out | Electron-electron rescaled distances | + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer function qmckl_compute_factor_een_rescaled_e_deriv_e_f(context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, & + coord_new, ee_distance, een_rescaled_e, een_rescaled_e_deriv_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) :: coord_new(elec_num,3,walk_num) + double precision , intent(in) :: ee_distance(elec_num,elec_num,walk_num) + double precision , intent(in) :: een_rescaled_e(0:cord_num,elec_num,elec_num,walk_num) + double precision , intent(out) :: een_rescaled_e_deriv_e(0:cord_num,elec_num,4,elec_num,walk_num) + double precision,dimension(:,:,:),allocatable :: elec_dist_deriv_e + double precision :: x, rij_inv, kappa_l + integer*8 :: i, j, k, l, nw, ii + + allocate(elec_dist_deriv_e(4,elec_num,elec_num)) + + 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 + + ! Prepare table of exponentiated distances raised to appropriate power + een_rescaled_e_deriv_e = 0.0d0 + do nw = 1, walk_num + do j = 1, elec_num + do i = 1, elec_num + rij_inv = 1.0d0 / ee_distance(i, j, nw) + do ii = 1, 3 + elec_dist_deriv_e(ii, i, j) = (coord_new(i, ii, nw) - coord_new(j, ii, nw)) * rij_inv + end do + elec_dist_deriv_e(4, i, j) = 2.0d0 * rij_inv + end do + elec_dist_deriv_e(:, j, j) = 0.0d0 + end do + + ! prepare the actual een table + do l = 1, cord_num + kappa_l = - dble(l) * rescale_factor_kappa_ee + do j = 1, elec_num + do i = 1, elec_num + do ii = 1, 4 + een_rescaled_e_deriv_e(l, i, ii, j, nw) = kappa_l * elec_dist_deriv_e(ii, i, j) + end do + + een_rescaled_e_deriv_e(l, i, 4, j, nw) = een_rescaled_e_deriv_e(l, i, 4, j, nw) & + + een_rescaled_e_deriv_e(l, i, 1, j, nw) * een_rescaled_e_deriv_e(l, i, 1, j, nw) & + + een_rescaled_e_deriv_e(l, i, 2, j, nw) * een_rescaled_e_deriv_e(l, i, 2, j, nw) & + + een_rescaled_e_deriv_e(l, i, 3, j, nw) * een_rescaled_e_deriv_e(l, i, 3, j, nw) + + do ii = 1, 4 + een_rescaled_e_deriv_e(l, i, ii, j, nw) = een_rescaled_e_deriv_e(l, i, ii, j, nw) * & + een_rescaled_e(l, i, j, nw) + end do + end do + end do + end do + end do + +end function qmckl_compute_factor_een_rescaled_e_deriv_e_f + #+end_src + + #+CALL: generate_c_header(table=qmckl_factor_een_rescaled_e_deriv_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_factor_een_rescaled_e_deriv_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* coord_new, + const double* ee_distance, + const double* een_rescaled_e, + double* const een_rescaled_e_deriv_e ); + #+end_src + + + #+CALL: generate_c_interface(table=qmckl_factor_een_rescaled_e_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_e_deriv_e & + (context, & + walk_num, & + elec_num, & + cord_num, & + rescale_factor_kappa_ee, & + coord_new, & + ee_distance, & + een_rescaled_e, & + een_rescaled_e_deriv_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) :: coord_new(elec_num,3,walk_num) + real (c_double ) , intent(in) :: ee_distance(elec_num,elec_num,walk_num) + real (c_double ) , intent(in) :: een_rescaled_e(0:cord_num,elec_num,elec_num,walk_num) + real (c_double ) , intent(out) :: een_rescaled_e_deriv_e(0:cord_num,elec_num,4,elec_num,walk_num) + + integer(c_int32_t), external :: qmckl_compute_factor_een_rescaled_e_deriv_e_f + info = qmckl_compute_factor_een_rescaled_e_deriv_e_f & + (context, & + walk_num, & + elec_num, & + cord_num, & + rescale_factor_kappa_ee, & + coord_new, & + ee_distance, & + een_rescaled_e, & + een_rescaled_e_deriv_e) + + end function qmckl_compute_factor_een_rescaled_e_deriv_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): + 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): + 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)); + + #+end_src + +** Electron-nucleus rescaled distances for each order + + ~een_rescaled_n~ stores the table of the rescaled distances between + electrons and nucleii raised to the power \(p\) defined by ~cord_num~: + + \[ + C_{ia,p} = \left( 1 - \exp{-\kappa C_{ia}} \right)^p + \] + + where \(C_{ia}\) is the matrix of electron-nucleus distances. + +*** Get + + #+begin_src c :comments org :tangle (eval h_func) :noweb yes +qmckl_exit_code qmckl_get_jastrow_een_rescaled_n(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_n(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_n(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->nucleus.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1); + memcpy(distance_rescaled, ctx->jastrow.een_rescaled_n, 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_n(qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_provide_een_rescaled_n(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_en_distance(context); + if(rc != QMCKL_SUCCESS) return rc; + + /* Compute if necessary */ + if (ctx->date > ctx->jastrow.een_rescaled_n_date) { + + /* Allocate array */ + if (ctx->jastrow.een_rescaled_n == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->electron.num * ctx->nucleus.num * + ctx->electron.walk_num * (ctx->jastrow.cord_num + 1) * sizeof(double); + double* een_rescaled_n = (double*) qmckl_malloc(context, mem_info); + + if (een_rescaled_n == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_een_rescaled_n", + NULL); + } + ctx->jastrow.een_rescaled_n = een_rescaled_n; + } + + qmckl_exit_code rc = + qmckl_compute_een_rescaled_n(context, + ctx->electron.walk_num, + ctx->electron.num, + ctx->nucleus.num, + ctx->jastrow.cord_num, + ctx->electron.rescale_factor_kappa_en, + ctx->electron.en_distance, + ctx->jastrow.een_rescaled_n); + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->jastrow.een_rescaled_n_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + #+end_src + +*** Compute + :PROPERTIES: + :Name: qmckl_compute_een_rescaled_n + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_factor_een_rescaled_n_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 atoms | + | int64_t | cord_num | in | Order of polynomials | + | double | rescale_factor_kappa_en | in | Factor to rescale ee distances | + | double | en_distance[walk_num][elec_num][nucl_num] | in | Electron-nucleus distances | + | double | een_rescaled_n[walk_num][elec_num][nucl_num][0:cord_num] | out | Electron-nucleus rescaled distances | + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer function qmckl_compute_een_rescaled_n_f(context, walk_num, elec_num, nucl_num, cord_num, rescale_factor_kappa_en, & + en_distance, een_rescaled_n) & + 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) :: nucl_num + integer*8 , intent(in) :: cord_num + double precision , intent(in) :: rescale_factor_kappa_en + double precision , intent(in) :: en_distance(elec_num,nucl_num,walk_num) + double precision , intent(out) :: een_rescaled_n(0:cord_num,nucl_num,elec_num,walk_num) + double precision :: x + integer*8 :: i, a, k, l, nw + + 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 (cord_num <= 0) then + info = QMCKL_INVALID_ARG_5 + return + endif + + ! Prepare table of exponentiated distances raised to appropriate power + een_rescaled_n = 0.0d0 + do nw = 1, walk_num + + ! prepare the actual een table + een_rescaled_n(0, :, :, nw) = 1.0d0 + + do a = 1, nucl_num + do i = 1, elec_num + een_rescaled_n(1, a, i, nw) = dexp(-rescale_factor_kappa_en * en_distance(i, a, nw)) + end do + end do + + do l = 2, cord_num + do a = 1, nucl_num + do i = 1, elec_num + een_rescaled_n(l, a, i, nw) = een_rescaled_n(l - 1, a, i, nw) * een_rescaled_n(1, a, i, nw) + end do + end do + end do + end do + +end function qmckl_compute_een_rescaled_n_f + #+end_src + + #+CALL: generate_c_header(table=qmckl_factor_een_rescaled_n_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_n ( + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t nucl_num, + const int64_t cord_num, + const double rescale_factor_kappa_en, + const double* en_distance, + 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")) + + #+RESULTS: + #+begin_src f90 :tangle (eval f) :comments org :exports none + integer(c_int32_t) function qmckl_compute_een_rescaled_n & + (context, & + walk_num, & + elec_num, & + nucl_num, & + cord_num, & + rescale_factor_kappa_en, & + en_distance, & + een_rescaled_n) & + 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 :: cord_num + real (c_double ) , intent(in) , value :: rescale_factor_kappa_en + real (c_double ) , intent(in) :: en_distance(nucl_num,elec_num,walk_num) + real (c_double ) , intent(out) :: een_rescaled_n(0:cord_num,nucl_num,elec_num,walk_num) + + integer(c_int32_t), external :: qmckl_compute_een_rescaled_n_f + info = qmckl_compute_een_rescaled_n_f & + (context, & + walk_num, & + elec_num, & + nucl_num, & + cord_num, & + rescale_factor_kappa_en, & + en_distance, & + een_rescaled_n) + + end function qmckl_compute_een_rescaled_n + #+end_src + +*** Test + + #+begin_src python :results output :exports none :noweb yes +import numpy as np + +<> + +elec_coord = np.array(elec_coord)[0] +nucl_coord = np.array(nucl_coord) +elnuc_dist = np.zeros(shape=(elec_num, nucl_num),dtype=float) +for i in range(elec_num): + for a in range(nucl_num): + elnuc_dist[i, a] = np.linalg.norm(elec_coord[i] - nucl_coord[:,a]) + +kappa = 1.0 + +een_rescaled_n = np.zeros(shape=(nucl_num, elec_num, cord_num + 1), dtype=float) +een_rescaled_n[:,:,0] = 1.0 + +for a in range(nucl_num): + for i in range(elec_num): + een_rescaled_n[a, i, 1] = np.exp(-kappa * elnuc_dist[i, a]) + +for l in range(2,cord_num+1): + for a in range(nucl_num): + for i in range(elec_num): + een_rescaled_n[a, i, l] = een_rescaled_n[a, i, l - 1] * een_rescaled_n[a, i, 1] + +print(" een_rescaled_n[0, 2, 1] = ",een_rescaled_n[0, 2, 1]) +print(" een_rescaled_n[0, 3, 1] = ",een_rescaled_n[0, 3, 1]) +print(" een_rescaled_n[0, 4, 1] = ",een_rescaled_n[0, 4, 1]) +print(" een_rescaled_n[1, 3, 2] = ",een_rescaled_n[1, 3, 2]) +print(" een_rescaled_n[1, 4, 2] = ",een_rescaled_n[1, 4, 2]) +print(" een_rescaled_n[1, 5, 2] = ",een_rescaled_n[1, 5, 2]) + #+end_src + + #+RESULTS: + : een_rescaled_n[0, 2, 1] = 0.10612983920006765 + : een_rescaled_n[0, 3, 1] = 0.135652809635553 + : een_rescaled_n[0, 4, 1] = 0.023391817607642338 + : een_rescaled_n[1, 3, 2] = 0.880957224822116 + : een_rescaled_n[1, 4, 2] = 0.027185942659395074 + : een_rescaled_n[1, 5, 2] = 0.01343938025140174 + + #+begin_src c :tangle (eval c_test) +assert(qmckl_electron_provided(context)); + +double een_rescaled_n[walk_num][elec_num][nucl_num][(cord_num + 1)]; +rc = qmckl_get_jastrow_een_rescaled_n(context, &(een_rescaled_n[0][0][0][0])); + +// value of (0,2,1) +assert(fabs(een_rescaled_n[0][2][0][1]-0.10612983920006765) < 1.e-12); +assert(fabs(een_rescaled_n[0][3][0][1]-0.135652809635553) < 1.e-12); +assert(fabs(een_rescaled_n[0][4][0][1]-0.023391817607642338) < 1.e-12); +assert(fabs(een_rescaled_n[0][3][1][2]-0.880957224822116) < 1.e-12); +assert(fabs(een_rescaled_n[0][4][1][2]-0.027185942659395074) < 1.e-12); +assert(fabs(een_rescaled_n[0][5][1][2]-0.01343938025140174) < 1.e-12); + + #+end_src + +** Electron-nucleus rescaled distances for each order and derivatives + + ~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 +qmckl_exit_code qmckl_get_jastrow_een_rescaled_n_deriv_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_n_deriv_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_n_deriv_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 * 4 * ctx->nucleus.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1); + memcpy(distance_rescaled, ctx->jastrow.een_rescaled_n_deriv_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_n_deriv_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_n_deriv_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_en_distance(context); + if(rc != QMCKL_SUCCESS) return rc; + + /* Check if ee distance is provided */ + rc = qmckl_provide_een_rescaled_n(context); + if(rc != QMCKL_SUCCESS) return rc; + + /* Compute if necessary */ + if (ctx->date > ctx->jastrow.een_rescaled_n_deriv_e_date) { + + /* Allocate array */ + if (ctx->jastrow.een_rescaled_n_deriv_e == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->electron.num * 4 * ctx->nucleus.num * + ctx->electron.walk_num * (ctx->jastrow.cord_num + 1) * sizeof(double); + double* een_rescaled_n_deriv_e = (double*) qmckl_malloc(context, mem_info); + + if (een_rescaled_n_deriv_e == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_een_rescaled_n_deriv_e", + NULL); + } + ctx->jastrow.een_rescaled_n_deriv_e = een_rescaled_n_deriv_e; + } + + qmckl_exit_code rc = + qmckl_compute_factor_een_rescaled_n_deriv_e(context, + ctx->electron.walk_num, + ctx->electron.num, + ctx->nucleus.num, + ctx->jastrow.cord_num, + ctx->electron.rescale_factor_kappa_en, + ctx->electron.coord_new, + ctx->nucleus.coord, + ctx->electron.en_distance, + ctx->jastrow.een_rescaled_n, + ctx->jastrow.een_rescaled_n_deriv_e); + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->jastrow.een_rescaled_n_deriv_e_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + #+end_src + +*** Compute + :PROPERTIES: + :Name: qmckl_compute_factor_een_rescaled_n_deriv_e + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_compute_factor_een_rescaled_n_deriv_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 | nucl_num | in | Number of atoms | + | int64_t | cord_num | in | Order of polynomials | + | double | rescale_factor_kappa_en | in | Factor to rescale ee distances | + | double | coord_new[walk_num][3][elec_num] | in | Electron coordinates | + | double | coord[3][nucl_num] | in | Nuclear coordinates | + | double | en_distance[walk_num][elec_num][nucl_num] | in | Electron-nucleus distances | + | double | een_rescaled_n[walk_num][elec_num][nucl_num][0:cord_num] | in | Electron-nucleus distances | + | double | een_rescaled_n_deriv_e[walk_num][elec_num][4][nucl_num][0:cord_num] | out | Electron-nucleus rescaled distances | + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer function qmckl_compute_factor_een_rescaled_n_deriv_e_f(context, walk_num, elec_num, nucl_num, & + cord_num, rescale_factor_kappa_en, & + coord_new, coord, en_distance, een_rescaled_n, een_rescaled_n_deriv_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) :: nucl_num + integer*8 , intent(in) :: cord_num + double precision , intent(in) :: rescale_factor_kappa_en + double precision , intent(in) :: coord_new(elec_num,3,walk_num) + double precision , intent(in) :: coord(nucl_num,3) + double precision , intent(in) :: en_distance(elec_num,nucl_num,walk_num) + double precision , intent(in) :: een_rescaled_n(0:cord_num,nucl_num,elec_num,walk_num) + double precision , intent(out) :: een_rescaled_n_deriv_e(0:cord_num,nucl_num,4,elec_num,walk_num) + double precision,dimension(:,:,:),allocatable :: elnuc_dist_deriv_e + double precision :: x, ria_inv, kappa_l + integer*8 :: i, a, k, l, nw, ii + + allocate(elnuc_dist_deriv_e(4, elec_num, nucl_num)) + + 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 (cord_num <= 0) then + info = QMCKL_INVALID_ARG_5 + return + endif + + ! Prepare table of exponentiated distances raised to appropriate power + een_rescaled_n_deriv_e = 0.0d0 + do nw = 1, walk_num + + ! prepare the actual een table + do a = 1, nucl_num + do i = 1, elec_num + ria_inv = 1.0d0 / en_distance(i, a, nw) + do ii = 1, 3 + elnuc_dist_deriv_e(ii, i, a) = (coord_new(i, ii, nw) - coord(a, ii)) * ria_inv + end do + elnuc_dist_deriv_e(4, i, a) = 2.0d0 * ria_inv + end do + end do + + do l = 0, cord_num + kappa_l = - dble(l) * rescale_factor_kappa_en + do a = 1, nucl_num + do i = 1, elec_num + do ii = 1, 4 + een_rescaled_n_deriv_e(l, a, ii, i, nw) = kappa_l * elnuc_dist_deriv_e(ii, i, a) + end do + + een_rescaled_n_deriv_e(l, a, 4, i, nw) = een_rescaled_n_deriv_e(l, a, 4, i, nw) & + + een_rescaled_n_deriv_e(l, a, 1, i, nw) * een_rescaled_n_deriv_e(l, a, 1, i, nw) & + + een_rescaled_n_deriv_e(l, a, 2, i, nw) * een_rescaled_n_deriv_e(l, a, 2, i, nw) & + + een_rescaled_n_deriv_e(l, a, 3, i, nw) * een_rescaled_n_deriv_e(l, a, 3, i, nw) + + do ii = 1, 4 + een_rescaled_n_deriv_e(l, a, ii, i, nw) = een_rescaled_n_deriv_e(l, a, ii, i, nw) * & + een_rescaled_n(l, a, i, nw) + end do + end do + end do + end do + end do + +end function qmckl_compute_factor_een_rescaled_n_deriv_e_f + #+end_src + + #+CALL: generate_c_header(table=qmckl_compute_factor_een_rescaled_n_deriv_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_factor_een_rescaled_n_deriv_e ( + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t nucl_num, + const int64_t cord_num, + const double rescale_factor_kappa_en, + const double* coord_new, + const double* coord, + const double* en_distance, + const double* een_rescaled_n, + 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 & + (context, & + walk_num, & + elec_num, & + nucl_num, & + cord_num, & + rescale_factor_kappa_en, & + coord_new, & + coord, & + en_distance, & + een_rescaled_n, & + een_rescaled_n_deriv_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 :: nucl_num + integer (c_int64_t) , intent(in) , value :: cord_num + real (c_double ) , intent(in) , value :: rescale_factor_kappa_en + real (c_double ) , intent(in) :: coord_new(elec_num,3,walk_num) + real (c_double ) , intent(in) :: coord(nucl_num,3) + real (c_double ) , intent(in) :: en_distance(nucl_num,elec_num,walk_num) + real (c_double ) , intent(in) :: een_rescaled_n(0:cord_num,nucl_num,elec_num,walk_num) + real (c_double ) , intent(out) :: een_rescaled_n_deriv_e(0:cord_num,nucl_num,4,elec_num,walk_num) + + integer(c_int32_t), external :: qmckl_compute_factor_een_rescaled_n_deriv_e_f + info = qmckl_compute_factor_een_rescaled_n_deriv_e_f & + (context, & + walk_num, & + elec_num, & + nucl_num, & + cord_num, & + rescale_factor_kappa_en, & + coord_new, & + coord, & + en_distance, & + een_rescaled_n, & + een_rescaled_n_deriv_e) + + end function qmckl_compute_factor_een_rescaled_n_deriv_e + #+end_src + +*** Test + + #+begin_src python :results output :exports none :noweb yes +import numpy as np + +<> + +elec_coord = np.array(elec_coord)[0] +nucl_coord = np.array(nucl_coord) +elnuc_dist = np.zeros(shape=(elec_num, nucl_num),dtype=float) +for i in range(elec_num): + for a in range(nucl_num): + elnuc_dist[i, a] = np.linalg.norm(elec_coord[i] - nucl_coord[:,a]) + +kappa = 1.0 + +een_rescaled_n = np.zeros(shape=(nucl_num, elec_num, cord_num + 1), dtype=float) +een_rescaled_n[:,:,0] = 1.0 + +for a in range(nucl_num): + for i in range(elec_num): + een_rescaled_n[a, i, 1] = np.exp(-kappa * elnuc_dist[i, a]) + +for l in range(2,cord_num+1): + for a in range(nucl_num): + for i in range(elec_num): + een_rescaled_n[a, i, l] = een_rescaled_n[a, i, l - 1] * een_rescaled_n[a, i, 1] + +print(" een_rescaled_n[0, 2, 1] = ",een_rescaled_n[0, 2, 1]) +print(" een_rescaled_n[0, 3, 1] = ",een_rescaled_n[0, 3, 1]) +print(" een_rescaled_n[0, 4, 1] = ",een_rescaled_n[0, 4, 1]) +print(" een_rescaled_n[1, 3, 2] = ",een_rescaled_n[1, 3, 2]) +print(" een_rescaled_n[1, 4, 2] = ",een_rescaled_n[1, 4, 2]) +print(" een_rescaled_n[1, 5, 2] = ",een_rescaled_n[1, 5, 2]) + #+end_src + + #+RESULTS: + : een_rescaled_n[0, 2, 1] = 0.10612983920006765 + : een_rescaled_n[0, 3, 1] = 0.135652809635553 + : een_rescaled_n[0, 4, 1] = 0.023391817607642338 + : een_rescaled_n[1, 3, 2] = 0.880957224822116 + : een_rescaled_n[1, 4, 2] = 0.027185942659395074 + : een_rescaled_n[1, 5, 2] = 0.01343938025140174 + + #+begin_src c :tangle (eval c_test) +//assert(qmckl_electron_provided(context)); + + #+end_src + +** Prepare for electron-electron-nucleus Jastrow \(f_{een}\) + + 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~. + +*** Get + + #+begin_src c :comments org :tangle (eval h_func) :noweb yes +qmckl_exit_code qmckl_get_jastrow_dim_cord_vect(qmckl_context context, int64_t* const dim_cord_vect); +qmckl_exit_code qmckl_get_jastrow_cord_vect_full(qmckl_context context, double* const cord_vect_full); +qmckl_exit_code qmckl_get_jastrow_lkpm_combined_index(qmckl_context context, int64_t* const lkpm_combined_index); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_get_jastrow_dim_cord_vect(qmckl_context context, int64_t* const dim_cord_vect) +{ + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + + qmckl_exit_code rc; + + rc = qmckl_provide_dim_cord_vect(context); + if (rc != QMCKL_SUCCESS) return rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + *dim_cord_vect = ctx->jastrow.dim_cord_vect; + + return QMCKL_SUCCESS; +} + +qmckl_exit_code qmckl_get_jastrow_cord_vect_full(qmckl_context context, double* const cord_vect_full) +{ + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + + qmckl_exit_code rc; + + rc = qmckl_provide_dim_cord_vect(context); + if (rc != QMCKL_SUCCESS) return rc; + + rc = qmckl_provide_cord_vect_full(context); + if (rc != QMCKL_SUCCESS) return rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + size_t sze = ctx->jastrow.dim_cord_vect * ctx->nucleus.num; + memcpy(cord_vect_full, ctx->jastrow.cord_vect_full, sze * sizeof(double)); + + return QMCKL_SUCCESS; +} + +qmckl_exit_code qmckl_get_jastrow_lkpm_combined_index(qmckl_context context, int64_t* const lkpm_combined_index) +{ + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + + qmckl_exit_code rc; + + rc = qmckl_provide_dim_cord_vect(context); + if (rc != QMCKL_SUCCESS) return rc; + + rc = qmckl_provide_cord_vect_full(context); + if (rc != QMCKL_SUCCESS) return rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + size_t sze = ctx->jastrow.dim_cord_vect * 4; + memcpy(lkpm_combined_index, ctx->jastrow.lkpm_combined_index, 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_dim_cord_vect(qmckl_context context); +qmckl_exit_code qmckl_provide_cord_vect_full(qmckl_context context); +qmckl_exit_code qmckl_provide_lkpm_combined_index(qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_provide_dim_cord_vect(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); + + /* Compute if necessary */ + if (ctx->date > ctx->jastrow.dim_cord_vect_date) { + + qmckl_exit_code rc = + qmckl_compute_dim_cord_vect(context, + ctx->jastrow.cord_num, + &(ctx->jastrow.dim_cord_vect)); + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->jastrow.dim_cord_vect_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + +qmckl_exit_code qmckl_provide_cord_vect_full(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 dim_cord_vect is provided */ + qmckl_exit_code rc = qmckl_provide_dim_cord_vect(context); + if(rc != QMCKL_SUCCESS) return rc; + + /* Compute if necessary */ + if (ctx->date > ctx->jastrow.cord_vect_full_date) { + + /* Allocate array */ + if (ctx->jastrow.cord_vect_full == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->jastrow.dim_cord_vect * ctx->nucleus.num * sizeof(double); + double* cord_vect_full = (double*) qmckl_malloc(context, mem_info); + + if (cord_vect_full == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_provide_cord_vect_full", + NULL); + } + ctx->jastrow.cord_vect_full = cord_vect_full; + } + + qmckl_exit_code rc = + qmckl_compute_cord_vect_full(context, + ctx->nucleus.num, + ctx->jastrow.cord_num, + ctx->jastrow.dim_cord_vect, + ctx->jastrow.type_nucl_num, + ctx->jastrow.type_nucl_vector, + ctx->jastrow.cord_vector, + ctx->jastrow.cord_vect_full); + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->jastrow.cord_vect_full_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + +qmckl_exit_code qmckl_provide_lkpm_combined_index(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 dim_cord_vect is provided */ + qmckl_exit_code rc = qmckl_provide_dim_cord_vect(context); + if(rc != QMCKL_SUCCESS) return rc; + + /* Compute if necessary */ + if (ctx->date > ctx->jastrow.lkpm_combined_index_date) { + + /* Allocate array */ + if (ctx->jastrow.lkpm_combined_index == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = 4 * ctx->jastrow.dim_cord_vect * sizeof(int64_t); + int64_t* lkpm_combined_index = (int64_t*) qmckl_malloc(context, mem_info); + + if (lkpm_combined_index == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_provide_lkpm_combined_index", + NULL); + } + ctx->jastrow.lkpm_combined_index = lkpm_combined_index; + } + + qmckl_exit_code rc = + qmckl_compute_lkpm_combined_index(context, + ctx->jastrow.cord_num, + ctx->jastrow.dim_cord_vect, + ctx->jastrow.lkpm_combined_index); + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->jastrow.lkpm_combined_index_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + #+end_src + +*** Compute dim_cord_vect + :PROPERTIES: + :Name: qmckl_compute_dim_cord_vect + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_factor_dim_cord_vect_args + | qmckl_context | context | in | Global state | + | int64_t | cord_num | in | Order of polynomials | + | int64_t | dim_cord_vect | out | dimension of cord_vect_full table | + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer function qmckl_compute_dim_cord_vect_f(context, cord_num, dim_cord_vect) & + result(info) + use qmckl + implicit none + integer(qmckl_context), intent(in) :: context + integer*8 , intent(in) :: cord_num + integer*8 , intent(out) :: dim_cord_vect + double precision :: x + integer*8 :: i, a, k, l, p, lmax + + info = QMCKL_SUCCESS + + if (context == QMCKL_NULL_CONTEXT) then + info = QMCKL_INVALID_CONTEXT + return + endif + + if (cord_num <= 0) then + info = QMCKL_INVALID_ARG_2 + return + endif + + dim_cord_vect = 0 + + do p = 2, cord_num + do k = p - 1, 0, -1 + if (k .ne. 0) then + lmax = p - k + else + lmax = p - k - 2 + endif + do l = lmax, 0, -1 + if (iand(p - k - l, 1_8) == 1) cycle + dim_cord_vect = dim_cord_vect + 1 + end do + end do + end do + +end function qmckl_compute_dim_cord_vect_f + #+end_src + + #+CALL: generate_c_header(table=qmckl_factor_dim_cord_vect_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + + #+RESULTS: + #+begin_src c :tangle (eval h_func) :comments org + qmckl_exit_code qmckl_compute_dim_cord_vect ( + const qmckl_context context, + const int64_t cord_num, + int64_t* const dim_cord_vect ); + #+end_src + + + #+CALL: generate_c_interface(table=qmckl_factor_dim_cord_vect_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_dim_cord_vect & + (context, cord_num, dim_cord_vect) & + 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 :: cord_num + integer (c_int64_t) , intent(out) :: dim_cord_vect + + integer(c_int32_t), external :: qmckl_compute_dim_cord_vect_f + info = qmckl_compute_dim_cord_vect_f & + (context, cord_num, dim_cord_vect) + + end function qmckl_compute_dim_cord_vect + #+end_src + +*** Compute cord_vect_full + :PROPERTIES: + :Name: qmckl_compute_cord_vect_full + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_factor_cord_vect_full_args + | qmckl_context | context | in | Global state | + | int64_t | nucl_num | in | Number of atoms | + | int64_t | cord_num | in | Order of polynomials | + | int64_t | dim_cord_vect | in | dimension of cord full table | + | int64_t | type_nucl_num | in | dimension of cord full table | + | int64_t | type_nucl_vector[nucl_num] | in | dimension of cord full table | + | double | cord_vector[cord_num][type_nucl_num] | in | dimension of cord full table | + | double | cord_vect_full[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, cord_num, dim_cord_vect, type_nucl_num, & + type_nucl_vector, cord_vector, cord_vect_full) & + result(info) + use qmckl + implicit none + integer(qmckl_context), intent(in) :: context + integer*8 , intent(in) :: nucl_num + integer*8 , intent(in) :: cord_num + integer*8 , intent(in) :: dim_cord_vect + integer*8 , intent(in) :: type_nucl_num + integer*8 , intent(in) :: type_nucl_vector(nucl_num) + double precision , intent(in) :: cord_vector(cord_num, type_nucl_num) + double precision , intent(out) :: cord_vect_full(nucl_num,dim_cord_vect) + double precision :: x + integer*8 :: i, a, k, l, nw + + info = QMCKL_SUCCESS + + if (context == QMCKL_NULL_CONTEXT) then + info = QMCKL_INVALID_CONTEXT + return + endif + + if (nucl_num <= 0) then + info = QMCKL_INVALID_ARG_2 + return + endif + + if (cord_num <= 0) then + info = QMCKL_INVALID_ARG_3 + return + endif + + if (type_nucl_num <= 0) then + info = QMCKL_INVALID_ARG_4 + return + endif + + if (dim_cord_vect <= 0) then + info = QMCKL_INVALID_ARG_5 + return + endif + + + do a = 1, nucl_num + cord_vect_full(1:dim_cord_vect,a) = cord_vector(1:dim_cord_vect,type_nucl_vector(a)) + end do + +end function qmckl_compute_cord_vect_full_f + #+end_src + + #+CALL: generate_c_header(table=qmckl_factor_cord_vect_full_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + + #+RESULTS: + #+begin_src c :tangle (eval h_func) :comments org + qmckl_exit_code qmckl_compute_cord_vect_full ( + const qmckl_context context, + const int64_t nucl_num, + const int64_t cord_num, + const int64_t dim_cord_vect, + const int64_t type_nucl_num, + const int64_t* type_nucl_vector, + const double* cord_vector, + double* const cord_vect_full ); + #+end_src + + + #+CALL: generate_c_interface(table=qmckl_factor_cord_vect_full_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_cord_vect_full & + (context, & + nucl_num, & + cord_num, & + dim_cord_vect, & + type_nucl_num, & + type_nucl_vector, & + cord_vector, & + cord_vect_full) & + 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 :: nucl_num + integer (c_int64_t) , intent(in) , value :: cord_num + integer (c_int64_t) , intent(in) , value :: dim_cord_vect + integer (c_int64_t) , intent(in) , value :: type_nucl_num + integer (c_int64_t) , intent(in) :: type_nucl_vector(nucl_num) + real (c_double ) , intent(in) :: cord_vector(type_nucl_num,cord_num) + real (c_double ) , intent(out) :: cord_vect_full(nucl_num,dim_cord_vect) + + integer(c_int32_t), external :: qmckl_compute_cord_vect_full_f + info = qmckl_compute_cord_vect_full_f & + (context, & + nucl_num, & + cord_num, & + dim_cord_vect, & + type_nucl_num, & + type_nucl_vector, & + cord_vector, & + cord_vect_full) + + end function qmckl_compute_cord_vect_full + #+end_src + +*** Compute lkpm_combined_index + :PROPERTIES: + :Name: qmckl_compute_lkpm_combined_index + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_factor_lkpm_combined_index_args + | qmckl_context | context | in | Global state | + | int64_t | cord_num | in | Order of polynomials | + | int64_t | dim_cord_vect | in | dimension of cord full table | + | int64_t | lpkm_combined_index[4][dim_cord_vect] | out | Full list of combined indices | + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer function qmckl_compute_lkpm_combined_index_f(context, cord_num, dim_cord_vect, & + lkpm_combined_index) & + result(info) + use qmckl + implicit none + integer(qmckl_context), intent(in) :: context + integer*8 , intent(in) :: cord_num + integer*8 , intent(in) :: dim_cord_vect + integer*8 , intent(out) :: lkpm_combined_index(dim_cord_vect, 4) + double precision :: x + integer*8 :: i, a, k, l, kk, p, lmax, m + + info = QMCKL_SUCCESS + + if (context == QMCKL_NULL_CONTEXT) then + info = QMCKL_INVALID_CONTEXT + return + endif + + if (cord_num <= 0) then + info = QMCKL_INVALID_ARG_2 + return + endif + + if (dim_cord_vect <= 0) then + info = QMCKL_INVALID_ARG_3 + return + endif + + + kk = 0 + do p = 2, cord_num + do k = p - 1, 0, -1 + if (k .ne. 0) then + lmax = p - k + else + lmax = p - k - 2 + end if + do l = lmax, 0, -1 + if (iand(p - k - l, 1_8) .eq. 1) cycle + m = (p - k - l)/2 + kk = kk + 1 + lkpm_combined_index(kk, 1) = l + lkpm_combined_index(kk, 2) = k + lkpm_combined_index(kk, 3) = p + lkpm_combined_index(kk, 4) = m + end do + end do + end do + +end function qmckl_compute_lkpm_combined_index_f + #+end_src + + #+CALL: generate_c_header(table=qmckl_factor_lkpm_combined_index_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + + #+RESULTS: + #+begin_src c :tangle (eval h_func) :comments org + qmckl_exit_code qmckl_compute_lkpm_combined_index ( + const qmckl_context context, + const int64_t cord_num, + const int64_t dim_cord_vect, + int64_t* const lpkm_combined_index ); + #+end_src + + + #+CALL: generate_c_interface(table=qmckl_factor_lkpm_combined_index_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_lkpm_combined_index & + (context, cord_num, dim_cord_vect, lpkm_combined_index) & + 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 :: cord_num + integer (c_int64_t) , intent(in) , value :: dim_cord_vect + integer (c_int64_t) , intent(out) :: lpkm_combined_index(dim_cord_vect,4) + + integer(c_int32_t), external :: qmckl_compute_lkpm_combined_index_f + info = qmckl_compute_lkpm_combined_index_f & + (context, cord_num, dim_cord_vect, lpkm_combined_index) + + end function qmckl_compute_lkpm_combined_index + #+end_src + + +*** Test + + #+begin_src python :results output :exports none :noweb yes +import numpy as np + +<> + +elec_coord = np.array(elec_coord)[0] +nucl_coord = np.array(nucl_coord) +elnuc_dist = np.zeros(shape=(elec_num, nucl_num),dtype=float) +for i in range(elec_num): + for a in range(nucl_num): + elnuc_dist[i, a] = np.linalg.norm(elec_coord[i] - nucl_coord[:,a]) + +kappa = 1.0 + +een_rescaled_n = np.zeros(shape=(nucl_num, elec_num, cord_num + 1), dtype=float) +een_rescaled_n[:,:,0] = 1.0 + +for a in range(nucl_num): + for i in range(elec_num): + een_rescaled_n[a, i, 1] = np.exp(-kappa * elnuc_dist[i, a]) + +for l in range(2,cord_num+1): + for a in range(nucl_num): + for i in range(elec_num): + een_rescaled_n[a, i, l] = een_rescaled_n[a, i, l - 1] * een_rescaled_n[a, i, 1] + +print(" een_rescaled_n[0, 2, 1] = ",een_rescaled_n[0, 2, 1]) +print(" een_rescaled_n[0, 3, 1] = ",een_rescaled_n[0, 3, 1]) +print(" een_rescaled_n[0, 4, 1] = ",een_rescaled_n[0, 4, 1]) +print(" een_rescaled_n[1, 3, 2] = ",een_rescaled_n[1, 3, 2]) +print(" een_rescaled_n[1, 4, 2] = ",een_rescaled_n[1, 4, 2]) +print(" een_rescaled_n[1, 5, 2] = ",een_rescaled_n[1, 5, 2]) + #+end_src + + #+RESULTS: + : een_rescaled_n[0, 2, 1] = 0.10612983920006765 + : een_rescaled_n[0, 3, 1] = 0.135652809635553 + : een_rescaled_n[0, 4, 1] = 0.023391817607642338 + : een_rescaled_n[1, 3, 2] = 0.880957224822116 + : een_rescaled_n[1, 4, 2] = 0.027185942659395074 + : een_rescaled_n[1, 5, 2] = 0.01343938025140174 + + #+begin_src c :tangle (eval c_test) +//assert(qmckl_electron_provided(context)); +// + + #+end_src + +** Electron-electron-nucleus Jastrow \(f_{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); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_get_jastrow_factor_een(qmckl_context context, double* const factor_een) +{ + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + + qmckl_exit_code rc; + + rc = qmckl_provide_factor_een(context); + if (rc != QMCKL_SUCCESS) return rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + int64_t sze = ctx->electron.walk_num * ctx->electron.num; + memcpy(factor_een, ctx->jastrow.factor_een, 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_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) +{ + + 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_een_rescaled_e(context); + if(rc != QMCKL_SUCCESS) return rc; + + /* Check if en rescaled distance derivatives is provided */ + rc = qmckl_provide_een_rescaled_n(context); + if(rc != QMCKL_SUCCESS) return rc; + + /* Check if en rescaled distance derivatives is provided */ + rc = qmckl_provide_cord_vect_full(context); + if(rc != QMCKL_SUCCESS) return rc; + + /* Check if en rescaled distance derivatives is provided */ + rc = qmckl_provide_lkpm_combined_index(context); + if(rc != QMCKL_SUCCESS) return rc; + + /* Compute if necessary */ + if (ctx->date > ctx->jastrow.factor_een_date) { + + /* Allocate array */ + if (ctx->jastrow.factor_een == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->electron.walk_num * sizeof(double); + double* factor_een = (double*) qmckl_malloc(context, mem_info); + + if (factor_een == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_provide_factor_een", + NULL); + } + ctx->jastrow.factor_een = factor_een; + } + + qmckl_exit_code rc = + qmckl_compute_factor_een(context, + ctx->electron.walk_num, + ctx->electron.num, + ctx->nucleus.num, + ctx->jastrow.cord_num, + ctx->jastrow.dim_cord_vect, + ctx->jastrow.cord_vect_full, + ctx->jastrow.lkpm_combined_index, + ctx->jastrow.een_rescaled_e, + ctx->jastrow.een_rescaled_n, + ctx->jastrow.factor_een); + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->jastrow.factor_een_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + #+end_src + +*** Compute + :PROPERTIES: + :Name: qmckl_compute_factor_een + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_factor_een_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 | cord_num | in | order of polynomials | + | int64_t | dim_cord_vect | in | dimension of full coefficient vector | + | double | cord_vect_full[dim_cord_vect][nucl_num] | in | full coefficient vector | + | int64_t | lkpm_combined_index[4][dim_cord_vect] | in | combined indices | + | double | een_rescaled_e[walk_num][elec_num][elec_num][0:cord_num] | in | Electron-nucleus rescaled | + | double | een_rescaled_n[walk_num][elec_num][nucl_num][0:cord_num] | in | Electron-nucleus rescaled factor | + | double | factor_een[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, & + een_rescaled_e, een_rescaled_n, factor_een) & + result(info) + use qmckl + implicit none + integer(qmckl_context), intent(in) :: context + integer*8 , intent(in) :: walk_num, elec_num, cord_num, nucl_num, dim_cord_vect + integer*8 , intent(in) :: lkpm_combined_index(4,dim_cord_vect) + double precision , intent(in) :: cord_vect_full(dim_cord_vect, nucl_num) + double precision , intent(in) :: een_rescaled_e(walk_num, elec_num, elec_num, 0:cord_num) + double precision , intent(in) :: een_rescaled_n(walk_num, elec_num, nucl_num, 0:cord_num) + double precision , intent(out) :: factor_een(walk_num) + + integer*8 :: i, a, j, l, k, p, m, n, nw + double precision :: accu, accu2, cn + + 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 (cord_num <= 0) then + info = QMCKL_INVALID_ARG_5 + return + endif + + factor_een = 0.0d0 + + do nw =1, walk_num + do n = 1, dim_cord_vect + l = lkpm_combined_index(1, n) + k = lkpm_combined_index(2, n) + p = lkpm_combined_index(3, n) + m = lkpm_combined_index(4, n) + + do a = 1, nucl_num + accu2 = 0.0d0 + cn = cord_vect_full(n, a) + do j = 1, elec_num + accu = 0.0d0 + do i = 1, elec_num + accu = accu + een_rescaled_e(nw, i, j, k) * & + een_rescaled_n(nw, i, a, m) + end do + accu2 = accu2 + accu * een_rescaled_n(nw, j, a, m + l) + end do + factor_een(nw) = factor_een(nw) + accu2 + cn + end do + end do + end do + +end function qmckl_compute_factor_een_f + #+end_src + + #+CALL: generate_c_header(table=qmckl_factor_een_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_een ( + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t nucl_num, + const int64_t cord_num, + const int64_t dim_cord_vect, + const double* cord_vect_full, + const int64_t* lkpm_combined_index, + const double* een_rescaled_e, + const double* een_rescaled_n, + double* const factor_een ); + #+end_src + + + #+CALL: generate_c_interface(table=qmckl_factor_een_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 & + (context, & + walk_num, & + elec_num, & + nucl_num, & + cord_num, & + dim_cord_vect, & + cord_vect_full, & + lkpm_combined_index, & + een_rescaled_e, & + een_rescaled_n, & + factor_een) & + 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 :: cord_num + integer (c_int64_t) , intent(in) , value :: dim_cord_vect + real (c_double ) , intent(in) :: cord_vect_full(nucl_num,dim_cord_vect) + integer (c_int64_t) , intent(in) :: lkpm_combined_index(dim_cord_vect,4) + real (c_double ) , intent(in) :: een_rescaled_e(0:cord_num,elec_num,elec_num,walk_num) + real (c_double ) , intent(in) :: een_rescaled_n(0:cord_num,nucl_num,elec_num,walk_num) + real (c_double ) , intent(out) :: factor_een(walk_num) + + integer(c_int32_t), external :: qmckl_compute_factor_een_f + info = qmckl_compute_factor_een_f & + (context, & + walk_num, & + elec_num, & + nucl_num, & + cord_num, & + dim_cord_vect, & + cord_vect_full, & + lkpm_combined_index, & + een_rescaled_e, & + een_rescaled_n, & + factor_een) + + end function qmckl_compute_factor_een + #+end_src + +*** Test + #+begin_src python :results output :exports none :noweb yes +import numpy as np + +<> + +kappa = 1.0 + +elec_coord = np.array(elec_coord)[0] +nucl_coord = np.array(nucl_coord) +elnuc_dist = np.zeros(shape=(elec_num, nucl_num),dtype=float) +for i in range(elec_num): + for j in range(nucl_num): + elnuc_dist[i, j] = np.linalg.norm(elec_coord[i] - nucl_coord[:,j]) + +elnuc_dist_deriv_e = np.zeros(shape=(4, elec_num, nucl_num),dtype=float) +for a in range(nucl_num): + for i in range(elec_num): + rij_inv = 1.0 / elnuc_dist[i, a] + for ii in range(3): + elnuc_dist_deriv_e[ii, i, a] = (elec_coord[i][ii] - nucl_coord[ii][a]) * rij_inv + elnuc_dist_deriv_e[3, i, a] = 2.0 * rij_inv + +en_distance_rescaled_deriv_e = np.zeros(shape=(4,elec_num,nucl_num),dtype=float) +for a in range(nucl_num): + for i in range(elec_num): + f = 1.0 - kappa * en_distance_rescaled[i][a] + for ii in range(4): + en_distance_rescaled_deriv_e[ii][i][a] = elnuc_dist_deriv_e[ii][i][a] + en_distance_rescaled_deriv_e[3][i][a] = en_distance_rescaled_deriv_e[3][i][a] + \ + (-kappa * en_distance_rescaled_deriv_e[0][i][a] * en_distance_rescaled_deriv_e[0][i][a]) + \ + (-kappa * en_distance_rescaled_deriv_e[1][i][a] * en_distance_rescaled_deriv_e[1][i][a]) + \ + (-kappa * en_distance_rescaled_deriv_e[2][i][a] * en_distance_rescaled_deriv_e[2][i][a]) + for ii in range(4): + en_distance_rescaled_deriv_e[ii][i][a] = en_distance_rescaled_deriv_e[ii][i][a] * f + +third = 1.0 / 3.0 +factor_en_deriv_e = np.zeros(shape=(4,elec_num),dtype=float) +dx = np.zeros(shape=(4),dtype=float) +pow_ser_g = np.zeros(shape=(3),dtype=float) +for a in range(nucl_num): + for i in range(elec_num): + x = en_distance_rescaled[i][a] + if abs(x) < 1e-18: + continue + pow_ser_g = np.zeros(shape=(3),dtype=float) + den = 1.0 + aord_vector[1][type_nucl_vector[a]-1] * x + invden = 1.0 / den + invden2 = invden * invden + invden3 = invden2 * invden + xinv = 1.0 / (x + 1.0E-18) + + for ii in range(4): + dx[ii] = en_distance_rescaled_deriv_e[ii][i][a] + + lap1 = 0.0 + lap2 = 0.0 + lap3 = 0.0 + for ii in range(3): + x = en_distance_rescaled[i][a] + if x < 1e-18: + continue + for p in range(2,aord_num+1): + y = p * aord_vector[(p-1) + 1][type_nucl_vector[a]-1] * x + pow_ser_g[ii] = pow_ser_g[ii] + y * dx[ii] + 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] * \ + dx[ii] * invden2 + pow_ser_g[ii] + + ii = 3 + lap2 = lap2 * dx[ii] * third + lap3 = lap3 + den * dx[ii] + lap3 = lap3 * (aord_vector[0][type_nucl_vector[a]-1] * invden3) + factor_en_deriv_e[ii][i] = factor_en_deriv_e[ii][i] + lap1 + lap2 + lap3 + +print("factor_en_deriv_e[0][0]:",factor_en_deriv_e[0][0]) +print("factor_en_deriv_e[1][0]:",factor_en_deriv_e[1][0]) +print("factor_en_deriv_e[2][0]:",factor_en_deriv_e[2][0]) +print("factor_en_deriv_e[3][0]:",factor_en_deriv_e[3][0]) + + + #+end_src + + #+RESULTS: + : factor_en_deriv_e[0][0]: 0.11609919541763383 + : 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 */ +//assert(qmckl_jastrow_provided(context)); +// + + #+end_src + +** Electron-electron-nucleus Jastrow \(f_{een}\) derivative + + 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); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_get_jastrow_factor_een_deriv_e(qmckl_context context, double* const factor_een_deriv_e) +{ + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + + qmckl_exit_code rc; + + rc = qmckl_provide_factor_een_deriv_e(context); + if (rc != QMCKL_SUCCESS) return rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + int64_t sze = ctx->electron.walk_num * ctx->electron.num; + memcpy(factor_een_deriv_e, ctx->jastrow.factor_een_deriv_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_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) +{ + + 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_een_rescaled_e(context); + if(rc != QMCKL_SUCCESS) return rc; + + /* Check if en rescaled distance derivatives is provided */ + rc = qmckl_provide_een_rescaled_n(context); + if(rc != QMCKL_SUCCESS) return rc; + + /* Check if en rescaled distance is provided */ + rc = qmckl_provide_een_rescaled_e_deriv_e(context); + if(rc != QMCKL_SUCCESS) return rc; + + /* Check if en rescaled distance derivatives is provided */ + rc = qmckl_provide_een_rescaled_n_deriv_e(context); + if(rc != QMCKL_SUCCESS) return rc; + + /* Check if en rescaled distance derivatives is provided */ + rc = qmckl_provide_cord_vect_full(context); + if(rc != QMCKL_SUCCESS) return rc; + + /* Check if en rescaled distance derivatives is provided */ + rc = qmckl_provide_lkpm_combined_index(context); + if(rc != QMCKL_SUCCESS) return rc; + + /* Compute if necessary */ + if (ctx->date > ctx->jastrow.factor_een_deriv_e_date) { + + /* Allocate array */ + if (ctx->jastrow.factor_een_deriv_e == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = 4 * ctx->electron.num * ctx->electron.walk_num * sizeof(double); + double* factor_een_deriv_e = (double*) qmckl_malloc(context, mem_info); + + if (factor_een_deriv_e == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_provide_factor_een_deriv_e", + NULL); + } + 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, + ctx->electron.num, + ctx->nucleus.num, + ctx->jastrow.cord_num, + ctx->jastrow.dim_cord_vect, + ctx->jastrow.cord_vect_full, + ctx->jastrow.lkpm_combined_index, + ctx->jastrow.een_rescaled_e, + ctx->jastrow.een_rescaled_n, + ctx->jastrow.een_rescaled_e_deriv_e, + ctx->jastrow.een_rescaled_n_deriv_e, + ctx->jastrow.factor_een_deriv_e); + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->jastrow.factor_een_deriv_e_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + #+end_src + +*** Compute + :PROPERTIES: + :Name: qmckl_compute_factor_een_deriv_e + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_factor_een_deriv_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 | nucl_num | in | Number of nucleii | + | int64_t | cord_num | in | order of polynomials | + | int64_t | dim_cord_vect | in | dimension of full coefficient vector | + | double | cord_vect_full[dim_cord_vect][nucl_num] | in | full coefficient vector | + | int64_t | lkpm_combined_index[4][dim_cord_vect] | in | combined indices | + | double | een_rescaled_e[walk_num][elec_num][elec_num][0:cord_num] | in | Electron-nucleus rescaled | + | double | een_rescaled_n[walk_num][elec_num][nucl_num][0:cord_num] | in | Electron-nucleus rescaled factor | + | double | een_rescaled_e_deriv_e[walk_num][elec_num][4][elec_num][0:cord_num] | in | Electron-nucleus rescaled | + | double | een_rescaled_n_deriv_e[walk_num][elec_num][4][nucl_num][0:cord_num] | in | Electron-nucleus rescaled factor | + | double | factor_een_deriv_e[walk_num][4][elec_num] | out | 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, & + een_rescaled_e, een_rescaled_n, & + een_rescaled_e_deriv_e, een_rescaled_n_deriv_e, factor_een_deriv_e) & + result(info) + use qmckl + implicit none + integer(qmckl_context), intent(in) :: context + integer*8 , intent(in) :: walk_num, elec_num, cord_num, nucl_num, dim_cord_vect + integer*8 , intent(in) :: lkpm_combined_index(4,dim_cord_vect) + double precision , intent(in) :: cord_vect_full(dim_cord_vect, nucl_num) + double precision , intent(in) :: een_rescaled_e(walk_num, elec_num, elec_num, 0:cord_num) + double precision , intent(in) :: een_rescaled_n(walk_num, elec_num, nucl_num, 0:cord_num) + double precision , intent(in) :: een_rescaled_e_deriv_e(walk_num, elec_num, 4, elec_num, 0:cord_num) + double precision , intent(in) :: een_rescaled_n_deriv_e(walk_num, elec_num, 4, nucl_num, 0:cord_num) + double precision , intent(out) :: factor_een_deriv_e(elec_num, 4, walk_num) + + integer*8 :: i, a, j, l, k, p, m, n, nw + double precision :: accu, accu2, cn + double precision :: daccu(1:4), daccu2(1:4) + + 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 (cord_num <= 0) then + info = QMCKL_INVALID_ARG_5 + return + endif + + factor_een_deriv_e = 0.0d0 + + do nw =1, walk_num + do n = 1, dim_cord_vect + l = lkpm_combined_index(1, n) + k = lkpm_combined_index(2, n) + p = lkpm_combined_index(3, n) + m = lkpm_combined_index(4, n) + + do a = 1, nucl_num + cn = cord_vect_full(n, a) + do j = 1, elec_num + accu = 0.0d0 + accu2 = 0.0d0 + daccu = 0.0d0 + daccu2 = 0.0d0 + do i = 1, elec_num + accu = accu + een_rescaled_e(nw, i, j, k) * & + een_rescaled_n(nw, i, a, m) + accu2 = accu2 + een_rescaled_e(nw, i, j, k) * & + een_rescaled_n(nw, i, a, m + l) + daccu(1:4) = daccu(1:4) + een_rescaled_e_deriv_e(nw, j, 1:4, i, k) * & + een_rescaled_n(nw, i, a, m) + daccu2(1:4) = daccu2(1:4) + een_rescaled_e_deriv_e(nw, j, 1:4, i, k) * & + een_rescaled_n(nw, i, a, m + l) + end do + factor_een_deriv_e(j, 1:4, nw) = factor_een_deriv_e(j, 1:4, nw) + & + (accu * een_rescaled_n_deriv_e(nw, j, 1:4, a, m + l) & + + daccu(1:4) * een_rescaled_n(nw, j, a, m + l) & + + daccu2(1:4) * een_rescaled_n(nw, j, a, m) & + + accu2 * een_rescaled_n_deriv_e(nw, j, 1:4, a, m)) * cn + + factor_een_deriv_e(j, 4, nw) = factor_een_deriv_e(j, 4, nw) + 2.0d0 * ( & + daccu (1) * een_rescaled_n_deriv_e(nw, j, 1, a, m + l) + & + daccu (2) * een_rescaled_n_deriv_e(nw, j, 2, a, m + l) + & + daccu (3) * een_rescaled_n_deriv_e(nw, j, 3, a, m + l) + & + daccu2(1) * een_rescaled_n_deriv_e(nw, j, 1, a, m ) + & + daccu2(2) * een_rescaled_n_deriv_e(nw, j, 2, a, m ) + & + daccu2(3) * een_rescaled_n_deriv_e(nw, j, 3, a, m ) ) * cn + + end do + end do + end do + end do + +end function qmckl_compute_factor_een_deriv_e_f + #+end_src + + #+CALL: generate_c_header(table=qmckl_factor_een_deriv_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_factor_een_deriv_e ( + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t nucl_num, + const int64_t cord_num, + const int64_t dim_cord_vect, + const double* cord_vect_full, + const int64_t* lkpm_combined_index, + const double* een_rescaled_e, + 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 ); + #+end_src + + + #+CALL: generate_c_interface(table=qmckl_factor_een_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_deriv_e & + (context, & + walk_num, & + elec_num, & + nucl_num, & + cord_num, & + dim_cord_vect, & + cord_vect_full, & + lkpm_combined_index, & + een_rescaled_e, & + een_rescaled_n, & + een_rescaled_e_deriv_e, & + een_rescaled_n_deriv_e, & + factor_een_deriv_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 :: nucl_num + integer (c_int64_t) , intent(in) , value :: cord_num + integer (c_int64_t) , intent(in) , value :: dim_cord_vect + real (c_double ) , intent(in) :: cord_vect_full(nucl_num,dim_cord_vect) + integer (c_int64_t) , intent(in) :: lkpm_combined_index(dim_cord_vect,4) + real (c_double ) , intent(in) :: een_rescaled_e(0:cord_num,elec_num,elec_num,walk_num) + real (c_double ) , intent(in) :: een_rescaled_n(0:cord_num,nucl_num,elec_num,walk_num) + real (c_double ) , intent(in) :: een_rescaled_e_deriv_e(0:cord_num,elec_num,4,elec_num,walk_num) + real (c_double ) , intent(in) :: een_rescaled_n_deriv_e(0:cord_num,nucl_num,4,elec_num,walk_num) + real (c_double ) , intent(out) :: factor_een_deriv_e(elec_num,4,walk_num) + + integer(c_int32_t), external :: qmckl_compute_factor_een_deriv_e_f + info = 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, & + een_rescaled_e, & + een_rescaled_n, & + een_rescaled_e_deriv_e, & + een_rescaled_n_deriv_e, & + factor_een_deriv_e) + + end function qmckl_compute_factor_een_deriv_e + #+end_src + +*** Test + #+begin_src python :results output :exports none :noweb yes +import numpy as np + +<> + +kappa = 1.0 + +elec_coord = np.array(elec_coord)[0] +nucl_coord = np.array(nucl_coord) +elnuc_dist = np.zeros(shape=(elec_num, nucl_num),dtype=float) +for i in range(elec_num): + for j in range(nucl_num): + elnuc_dist[i, j] = np.linalg.norm(elec_coord[i] - nucl_coord[:,j]) + +elnuc_dist_deriv_e = np.zeros(shape=(4, elec_num, nucl_num),dtype=float) +for a in range(nucl_num): + for i in range(elec_num): + rij_inv = 1.0 / elnuc_dist[i, a] + for ii in range(3): + elnuc_dist_deriv_e[ii, i, a] = (elec_coord[i][ii] - nucl_coord[ii][a]) * rij_inv + elnuc_dist_deriv_e[3, i, a] = 2.0 * rij_inv + +en_distance_rescaled_deriv_e = np.zeros(shape=(4,elec_num,nucl_num),dtype=float) +for a in range(nucl_num): + for i in range(elec_num): + f = 1.0 - kappa * en_distance_rescaled[i][a] + for ii in range(4): + en_distance_rescaled_deriv_e[ii][i][a] = elnuc_dist_deriv_e[ii][i][a] + en_distance_rescaled_deriv_e[3][i][a] = en_distance_rescaled_deriv_e[3][i][a] + \ + (-kappa * en_distance_rescaled_deriv_e[0][i][a] * en_distance_rescaled_deriv_e[0][i][a]) + \ + (-kappa * en_distance_rescaled_deriv_e[1][i][a] * en_distance_rescaled_deriv_e[1][i][a]) + \ + (-kappa * en_distance_rescaled_deriv_e[2][i][a] * en_distance_rescaled_deriv_e[2][i][a]) + for ii in range(4): + en_distance_rescaled_deriv_e[ii][i][a] = en_distance_rescaled_deriv_e[ii][i][a] * f + +third = 1.0 / 3.0 +factor_en_deriv_e = np.zeros(shape=(4,elec_num),dtype=float) +dx = np.zeros(shape=(4),dtype=float) +pow_ser_g = np.zeros(shape=(3),dtype=float) +for a in range(nucl_num): + for i in range(elec_num): + x = en_distance_rescaled[i][a] + if abs(x) < 1e-18: + continue + pow_ser_g = np.zeros(shape=(3),dtype=float) + den = 1.0 + aord_vector[1][type_nucl_vector[a]-1] * x + invden = 1.0 / den + invden2 = invden * invden + invden3 = invden2 * invden + xinv = 1.0 / (x + 1.0E-18) + + for ii in range(4): + dx[ii] = en_distance_rescaled_deriv_e[ii][i][a] + + lap1 = 0.0 + lap2 = 0.0 + lap3 = 0.0 + for ii in range(3): + x = en_distance_rescaled[i][a] + if x < 1e-18: + continue + for p in range(2,aord_num+1): + y = p * aord_vector[(p-1) + 1][type_nucl_vector[a]-1] * x + pow_ser_g[ii] = pow_ser_g[ii] + y * dx[ii] + 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] * \ + dx[ii] * invden2 + pow_ser_g[ii] + + ii = 3 + lap2 = lap2 * dx[ii] * third + lap3 = lap3 + den * dx[ii] + lap3 = lap3 * (aord_vector[0][type_nucl_vector[a]-1] * invden3) + factor_en_deriv_e[ii][i] = factor_en_deriv_e[ii][i] + lap1 + lap2 + lap3 + +print("factor_en_deriv_e[0][0]:",factor_en_deriv_e[0][0]) +print("factor_en_deriv_e[1][0]:",factor_en_deriv_e[1][0]) +print("factor_en_deriv_e[2][0]:",factor_en_deriv_e[2][0]) +print("factor_en_deriv_e[3][0]:",factor_en_deriv_e[3][0]) + + + #+end_src + + #+RESULTS: + : factor_en_deriv_e[0][0]: 0.11609919541763383 + : 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 */ + + #+end_src + +* End of files :noexport: + + #+begin_src c :tangle (eval h_private_type) +#endif + #+end_src + +*** Test + #+begin_src c :tangle (eval c_test) + rc = qmckl_context_destroy(context); + assert (rc == QMCKL_SUCCESS); + + return 0; +} + #+end_src + + +# -*- mode: org -*- +# vim: syntax=c + + + diff --git a/org/qmckl_tests.org b/org/qmckl_tests.org index a2fbbfe..0fb684b 100644 --- a/org/qmckl_tests.org +++ b/org/qmckl_tests.org @@ -1013,3 +1013,159 @@ double chbrclf_elec_coord[chbrclf_walk_num][chbrclf_elec_num][3] = { { #+END_src + + +* N2 + + This test is mainly for the Jastrow factor and was supplied by + Ramon Panades Baruetta. The coordinates and Jastrow coefficients + have been taken from his fork of IRPJast. The core electrons are + treated by pseudopotentials thus excluded from the actual calculation. + + | Number of atoms | 2 | + | Number of alpha electrons | 5 | + | Number of beta electrons | 5 | + | Number of core electrons | 4 | + +** XYZ coordinates + +#+BEGIN_example + 2 +N2 + N 0.000000 0.000000 0.000000 + N 0.000000 0.000000 2.059801 +#+END_example + + Nuclear coordinates are stored in atomic units in transposed format. + +#+begin_src c :tangle ../tests/n2.h +#define n2_nucl_num ((int64_t) 2) + +double n2_charge[n2_nucl_num] = { 5., 5.}; + +double n2_nucl_coord[3][n2_nucl_num] = +{ {0.000000, 0.000000 }, + {0.000000, 0.000000 }, + {0.000000, 2.059801 } }; +#+end_src + +** Electron coordinates + + + Electron coordinates are stored in atomic units in normal format. + +#+begin_src c :tangle ../tests/n2.h +#define n2_elec_up_num ((int64_t) 5) +#define n2_elec_dn_num ((int64_t) 5) +#define n2_elec_num ((int64_t) 10) +#define n2_walk_num ((int64_t) 1) + +double n2_elec_coord[n2_walk_num][n2_elec_num][3] = { { + {-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}}}; + +#+end_src + +** Jastrow related data + + This test is mainly for the Jastrow factor and was supplied by + Ramon Panades Baruetta. + +#+begin_src c :tangle ../tests/n2.h +/* Jastrow related */ + +#define n2_type_nucl_num ((int64_t) 1) +#define n2_aord_num ((int64_t) 5) +#define n2_bord_num ((int64_t) 5) +#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. }, + {-0.380512}, + {-0.157996}, + {-0.031558}, + { 0.021512}}; + +double n2_bord_vector[n2_bord_num + 1] = { + 0.5 , + 0.15366 , + 0.0672262 , + 0.02157 , + 0.0073096 , + 0.002866 }; + +double n2_cord_vector[n2_cord_num][n2_type_nucl_num] = { + { 5.717020e-01}, + {-5.142530e-01}, + {-5.130430e-01}, + { 9.486000e-03}, + {-4.205000e-03}, + { 4.263258e-01}, + { 8.288150e-02}, + { 5.118600e-03}, + {-2.997800e-03}, + {-5.270400e-03}, + {-7.500000e-05}, + {-8.301650e-02}, + { 1.454340e-02}, + { 5.143510e-02}, + { 9.250000e-04}, + {-4.099100e-03}, + { 4.327600e-03}, + {-1.654470e-03}, + { 2.614000e-03}, + {-1.477000e-03}, + {-1.137000e-03}, + {-4.010475e-02}, + { 6.106710e-03}}; + +double n2_cord_vector_full[n2_dim_cord_vec][n2_nucl_num] = { + { 5.717020e-01, 5.717020e-01}, + {-5.142530e-01, -5.142530e-01}, + {-5.130430e-01, -5.130430e-01}, + { 9.486000e-03, 9.486000e-03}, + {-4.205000e-03, -4.205000e-03}, + { 4.263258e-01, 4.263258e-01}, + { 8.288150e-02, 8.288150e-02}, + { 5.118600e-03, 5.118600e-03}, + {-2.997800e-03, -2.997800e-03}, + {-5.270400e-03, -5.270400e-03}, + {-7.500000e-05, -7.500000e-05}, + {-8.301650e-02, -8.301650e-02}, + { 1.454340e-02, 1.454340e-02}, + { 5.143510e-02, 5.143510e-02}, + { 9.250000e-04, 9.250000e-04}, + {-4.099100e-03, -4.099100e-03}, + { 4.327600e-03, 4.327600e-03}, + {-1.654470e-03, -1.654470e-03}, + { 2.614000e-03, 2.614000e-03}, + {-1.477000e-03, -1.477000e-03}, + {-1.137000e-03, -1.137000e-03}, + {-4.010475e-02, -4.010475e-02}, + { 6.106710e-03, 6.106710e-03}}; + +double n2_lkpm_of_cindex[4][n2_dim_cord_vec] = { + {1, 1, 2, 0, 0, 0, 2, 1, 1, 2, 3, 0, 2, 1, 3, 0, 0, 1, 3, 1, 1, 0, 3}, + {1, 1, 3, 4, 0, 2, 2, 4, 0, 0, 2, 4, 1, 3, 1, 4, 0, 1, 1, 4, 1, 2, 0}, + {4, 1, 0, 0, 4, 2, 1, 4, 5, 0, 2, 3, 5, 0, 0, 3, 5, 1, 3, 2, 5, 0, 1}, + {2, 5, 1, 4, 1, 5, 0, 2, 1, 5, 1, 0, 1, 5, 2, 3, 0, 5, 1, 1, 0, 5, 2}}; + +#+end_src + + +# -*- mode: org -*- +# vim: syntax=c diff --git a/org/table_of_contents b/org/table_of_contents index e9e2190..50a6159 100644 --- a/org/table_of_contents +++ b/org/table_of_contents @@ -6,6 +6,7 @@ qmckl_numprec.org qmckl_nucleus.org qmckl_electron.org qmckl_ao.org +qmckl_jastrow.org qmckl_distance.org qmckl_utils.org qmckl_tests.org diff --git a/tools/build_qmckl_h.sh b/tools/build_qmckl_h.sh index fdec31e..6365689 100755 --- a/tools/build_qmckl_h.sh +++ b/tools/build_qmckl_h.sh @@ -43,29 +43,29 @@ cat << EOF > ${OUTPUT} * ------------------------------------------ * QMCkl - Quantum Monte Carlo kernel library * ------------------------------------------ - * + * * Documentation : https://trex-coe.github.io/qmckl * Issues : https://github.com/trex-coe/qmckl/issues - * + * * BSD 3-Clause License - * + * * Copyright (c) 2020, TREX Center of Excellence * All rights reserved. - * + * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions are met: - * + * * 1. Redistributions of source code must retain the above copyright notice, this * list of conditions and the following disclaimer. - * + * * 2. Redistributions in binary form must reproduce the above copyright notice, * this list of conditions and the following disclaimer in the documentation * and/or other materials provided with the distribution. - * + * * 3. Neither the name of the copyright holder nor the names of its * contributors may be used to endorse or promote products derived from * this software without specific prior written permission. - * + * * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE @@ -76,10 +76,10 @@ cat << EOF > ${OUTPUT} * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - * - * - * - * + * + * + * + * */ #ifndef __QMCKL_H__ @@ -92,7 +92,7 @@ EOF for i in ${HEADERS} do - header=${srcdir}/src/$i + header=${srcdir}/src/$i if [[ -f $header ]] ; then echo "/* $header */" >> ${OUTPUT} cat $header >> ${OUTPUT} diff --git a/tools/lib.org b/tools/lib.org index 3536c32..01d86fb 100644 --- a/tools/lib.org +++ b/tools/lib.org @@ -124,8 +124,7 @@ for d in parse_table(table): results += [ f" {const}{c_type} {name}" ] results=',\n'.join(results) -template = f"""{rettyp} {fname} ( -{results} ); """ +template = f"""{rettyp} {fname} (\n{results} ); """ return template #+END_SRC @@ -259,3 +258,4 @@ return results #+begin_src f90 :tangle (eval fh_func) :comments org :exports none #+end_src +