From 36293d892e6961a65006b193719dc9fcf213693c Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 3 Jun 2021 01:32:50 +0200 Subject: [PATCH 1/9] Clean configure --- .gitignore | 40 ++++++++++++++++++++-------------------- Makefile.am | 2 +- configure.ac | 21 ++++++++------------- org/.gitignore | 1 + org/qmckl_ao.org | 2 +- pkgconfig/qmckl.pc.in | 6 +++--- tools/build_doc.sh | 5 +++-- tools/build_makefile.py | 4 ++-- tools/tangle.sh | 5 +++-- tools/theme.setup | 2 +- 10 files changed, 43 insertions(+), 45 deletions(-) diff --git a/.gitignore b/.gitignore index 96b7a1b..61f6fbb 100644 --- a/.gitignore +++ b/.gitignore @@ -1,17 +1,10 @@ +*.la +*.lo .deps/ -compile -depcomp -generated.mk -stamp-h1 -libtool -ltmain.sh -missing -m4/libtool.m4 -m4/ltoptions.m4 -m4/ltsugar.m4 -m4/ltversion.m4 -m4/lt~obsolete.m4 +.dirstamp +.libs Makefile +Makefile.in aclocal.m4 autom4te.cache/ config.guess @@ -19,13 +12,20 @@ config.log config.status config.sub configure -install-sh -qmckl.pc -Makefile.in -test-driver -.libs -.dirstamp -*.la -*.lo +generated.mk +m4/libtool.m4 +m4/ltoptions.m4 +m4/ltsugar.m4 +m4/ltversion.m4 +m4/lt~obsolete.m4 qmckl-*.tar.gz qmckl.mod +qmckl.pc +stamp-h1 +tools/compile +tools/depcomp +tools/install-sh +tools/libtool +tools/ltmain.sh +tools/missing +tools/test-driver diff --git a/Makefile.am b/Makefile.am index b911ad8..110ff8a 100644 --- a/Makefile.am +++ b/Makefile.am @@ -159,7 +159,7 @@ tests/chbrclf.h: $(qmckl_h) generated.mk: $(ORG_FILES) - python $(srcdir)/tools/build_makefile.py + $(PYTHON) $(srcdir)/tools/build_makefile.py cppcheck: cppcheck.out diff --git a/configure.ac b/configure.ac index 7a6c365..c313f56 100644 --- a/configure.ac +++ b/configure.ac @@ -38,13 +38,13 @@ AC_PREREQ([2.69]) AC_INIT([qmckl],[0.1.1], [https://github.com/TREX-CoE/qmckl/issues], [], [https://trex-coe.github.io/qmckl/index.html]) +AC_CONFIG_AUX_DIR(tools) AM_INIT_AUTOMAKE([subdir-objects color-tests parallel-tests silent-rules 1.11]) AM_MAINTAINER_MODE() LT_INIT AC_CONFIG_SRCDIR([Makefile.in]) AC_CONFIG_HEADERS([include/config.h]) AC_CONFIG_MACRO_DIR([m4]) -#AM_MAINTAINER_MODE VERSION_MAJOR=`echo ${PACKAGE_VERSION} | cut -d. -f1` VERSION_MINOR=`echo ${PACKAGE_VERSION} | cut -d. -f2` @@ -64,6 +64,9 @@ AC_LANG(C) # Checks for programs. AC_PROG_CC +# Make sure the c compiler supports C99 +m4_version_prereq([2.70],[], [AC_PROG_CC_C99]) +AS_IF([test "$ac_cv_prog_cc_c99" = "no"], [AC_MSG_ERROR([The compiler does not support C99])]) AC_PROG_CC_C_O AC_PROG_FC AC_PROG_FC_C_O @@ -72,13 +75,10 @@ AC_FC_FREEFORM AC_PROG_LIBTOOL AC_PROG_INSTALL AC_PROG_LN_S -AC_PROG_CC_C99 PKG_PROG_PKG_CONFIG([]) PKG_LIBS="" PKG_CFLAGS="" -AC_SUBST([HAS_CPPCHECK]) - # Checks for libraries. AC_FC_LIBRARY_LDFLAGS @@ -178,15 +178,8 @@ if test "x${QMCKL_DEVEL}" != "x"; then QMCKL_DEVEL=" -- Developer mode" AC_PROG_AWK - AC_CHECK_PROGS([PYTHON],[python python3 python2],[no]) - if test x${PYTHON} == xno ; then - AC_MSG_ERROR([ - -------------------------------------------- - Error: Python is required to build makefiles - -------------------------------------------- - ]) - fi - python ${srcdir}/tools/build_makefile.py + AM_PATH_PYTHON + ${PYTHON} ${srcdir}/tools/build_makefile.py AC_CHECK_PROGS([EMACS],[emacs26 emacs],[no]) if test x${EMACS} == xno ; then @@ -211,6 +204,8 @@ LIBS="${LAPACK_LIBS} ${BLAS_LIBS} ${LIBS}" PKG_LIBS="${PKG_LIBS} ${LIBS}" AC_SUBST([PKG_LIBS]) AC_SUBST([PKG_CFLAGS]) +AC_SUBST([HAS_CPPCHECK]) + AC_CONFIG_FILES([Makefile pkgconfig/qmckl.pc diff --git a/org/.gitignore b/org/.gitignore index e0963f0..931ace3 100644 --- a/org/.gitignore +++ b/org/.gitignore @@ -1,2 +1,3 @@ +ltximg/ *.tangled *.exported diff --git a/org/qmckl_ao.org b/org/qmckl_ao.org index 906bf17..06565aa 100644 --- a/org/qmckl_ao.org +++ b/org/qmckl_ao.org @@ -94,6 +94,7 @@ int main() { | ~shell_factor~ | ~[shell_num]~ | Normalization factor for each shell | | ~exponent~ | ~[prim_num]~ | Array of exponents | | ~coefficient~ | ~[prim_num]~ | Array of coefficients | + | ~prim_factor~ | ~[prim_num]~ | Normalization factors of the primtives | For H_2 with the following basis set, @@ -1661,7 +1662,6 @@ end function test_qmckl_ao_gaussian_vgl #+end_src ** TODO Slater basis functions - ** TODO Radial functions on a grid * Combining radial and polynomial parts * End of files :noexport: diff --git a/pkgconfig/qmckl.pc.in b/pkgconfig/qmckl.pc.in index 0d82f6b..a629932 100644 --- a/pkgconfig/qmckl.pc.in +++ b/pkgconfig/qmckl.pc.in @@ -7,7 +7,7 @@ Name: @PACKAGE_NAME@ Description: Quantum Monte Carlo kernel library URL: https://github.com/trex-coe/qmckl Version: @PACKAGE_VERSION@ -Cflags: -I${includedir} @PKG_CFLAGS@ -Libs: -L${libdir} -lqmckl @PKG_LIBS@ -Requires: @PKG_BLAS@ +Cflags: -I${includedir} +Libs: -L${libdir} -lqmckl +Libs.private: @PKG_LIBS@ diff --git a/tools/build_doc.sh b/tools/build_doc.sh index 389e3f3..37935dd 100755 --- a/tools/build_doc.sh +++ b/tools/build_doc.sh @@ -19,7 +19,7 @@ function extract_doc() html=${DOCS}/html/$(basename ${org%.org}.html) text=${DOCS}/text/$(basename ${org%.org}.txt) - ./missing emacs --batch \ + ./tools/missing emacs --batch \ --load ${HTMLIZE} \ --load ${CONFIG_DOC} \ ${org} \ @@ -34,8 +34,9 @@ function extract_doc() for i in $@ do exported=${i%.org}.exported + exported=$(dirname $exported)/.$(basename $exported) NOW=$(date +"%m%d%H%M.%S") - extract_doc ${i} &> $exported + extract_doc ${i} > $exported # Make log file older than the exported files touch -t ${NOW} $exported diff --git a/tools/build_makefile.py b/tools/build_makefile.py index 5290b3f..6f1f30f 100755 --- a/tools/build_makefile.py +++ b/tools/build_makefile.py @@ -33,8 +33,8 @@ def main(): for org in glob("org/*.org"): i = org.split('/')[-1].rsplit(".",1)[0] - tangled = "org/"+i+".tangled" - exported = "org/"+i+".exported" + tangled = "org/."+i+".tangled" + exported = "org/."+i+".exported" c_test_x = "tests/test_"+i c_test_o = "tests/test_"+i+".$(OBJEXT)" f_test_o = "tests/test_"+i+"_f.$(OBJEXT)" diff --git a/tools/tangle.sh b/tools/tangle.sh index a94189d..0f28038 100755 --- a/tools/tangle.sh +++ b/tools/tangle.sh @@ -24,7 +24,7 @@ function tangle() elif [[ ${org_file} -ot ${f_file} ]] ; then return fi - ./missing \ + ./tools/missing \ emacs --batch ${org_file} \ --load=${PWD}/tools/config_tangle.el \ -f org-babel-tangle @@ -33,8 +33,9 @@ function tangle() for i in $@ do tangled=${i%.org}.tangled + tangled=$(dirname $tangled)/.$(basename $tangled) NOW=$(date +"%m%d%H%M.%S") - tangle ${i} &> $tangled + tangle ${i} > $tangled # Make log file older than the tangled files touch -t ${NOW} $tangled diff --git a/tools/theme.setup b/tools/theme.setup index 3c46c18..156662b 100644 --- a/tools/theme.setup +++ b/tools/theme.setup @@ -6,7 +6,7 @@ #+INFOJS_OPT: toc:t mouse:underline path:org-info.js #+HTML_HEAD: -#+STARTUP: align fold nodlcheck hidestars oddeven lognotestate +#+STARTUP: align fold nodlcheck hidestars oddeven lognotestate latexpreview #+AUTHOR: TREX CoE #+LANGUAGE: en From d88ade486b405b8c189f13f05262f9c1a86f6e41 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 3 Jun 2021 01:48:26 +0200 Subject: [PATCH 2/9] Clean cppcheck --- Makefile.am | 4 +++- org/qmckl_nucleus.org | 17 +++++++---------- 2 files changed, 10 insertions(+), 11 deletions(-) diff --git a/Makefile.am b/Makefile.am index 110ff8a..57ad20b 100644 --- a/Makefile.am +++ b/Makefile.am @@ -166,7 +166,9 @@ cppcheck: cppcheck.out cppcheck.out: $(qmckl_h) cd src/ && \ cppcheck --addon=cert -q --error-exitcode=0 \ - --enable=all \ + --enable=all --suppress="unusedStructMember"\ + --suppress="unusedFunction" \ + --suppress="missingIncludeSystem" \ --language=c --std=c99 -rp --platform=unix64 \ -I../include *.c *.h 2>../$@ diff --git a/org/qmckl_nucleus.org b/org/qmckl_nucleus.org index 8bd2893..004560f 100644 --- a/org/qmckl_nucleus.org +++ b/org/qmckl_nucleus.org @@ -392,18 +392,15 @@ qmckl_set_nucleus_kappa(qmckl_context context, const double rescale_factor_kappa <> //TODO: Check for small values of kappa - //if (rescale_factor_kappa == 0) { - // return qmckl_failwith( context, - // QMCKL_INVALID_ARG_2, - // "qmckl_set_nucleus_kappa", - // "rescale_factor_kappa cannot be 0"); - //} + if (rescale_factor_kappa <= 0.0) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_set_nucleus_kappa", + "rescale_factor_kappa cannot be 0"); + } int32_t mask = 1 << 2; - int64_t num; - qmckl_exit_code rc; - ctx->nucleus.rescale_factor_kappa = rescale_factor_kappa; <> @@ -836,7 +833,7 @@ end function qmckl_compute_nn_distance_rescaled_f #+begin_src c :tangle (eval h_private_func) :comments org :exports none qmckl_exit_code qmckl_compute_nn_distance_rescaled ( - const qmckl_context context, + const qmckl_context context, const int64_t nucl_num, const double rescale_factor_kappa, const double* coord, From eb09d47bc0064618ba52818e11673e10402cb58a Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 3 Jun 2021 01:50:51 +0200 Subject: [PATCH 3/9] Fixing Actions --- tools/install_htmlize.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tools/install_htmlize.sh b/tools/install_htmlize.sh index 29f8208..b069777 100755 --- a/tools/install_htmlize.sh +++ b/tools/install_htmlize.sh @@ -2,7 +2,7 @@ # # Installs the htmlize Emacs plugin -./missing git clone "https://github.com/hniksic/emacs-htmlize" +./tools/missing git clone "https://github.com/hniksic/emacs-htmlize" mv emacs-htmlize/htmlize.el $1 rm -rf emacs-htmlize From 58cd3f009d92d56c74078d4e0acad7e69e3cb690 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 3 Jun 2021 18:26:00 +0200 Subject: [PATCH 4/9] Added test data for ch3clbrf --- org/qmckl_ao.org | 101 +++++++++++-- org/qmckl_tests.org | 339 ++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 426 insertions(+), 14 deletions(-) diff --git a/org/qmckl_ao.org b/org/qmckl_ao.org index 06565aa..2e40ada 100644 --- a/org/qmckl_ao.org +++ b/org/qmckl_ao.org @@ -10,27 +10,30 @@ Gaussian ($p=2$): \[ R_s(\mathbf{r}) = \mathcal{N}_s |\mathbf{r}-\mathbf{R}_A|^{n_s} - \sum_{k=1}^{N_{\text{prim}}} a_{ks} + \sum_{k=1}^{N_{\text{prim}}} a_{ks}\, f_{ks} \exp \left( - \gamma_{ks} | \mathbf{r}-\mathbf{R}_A | ^p \right). \] In the case of Gaussian functions, $n_s$ is always zero. The normalization factor $\mathcal{N}_s$ ensures that all the functions -of the shell are normalized to unity. As this normalization requires -the ability to compute overlap integrals, it should be written in the -file to ensure that the file is self-contained and does not require -the client program to have the ability to compute such integrals. +of the shell are normalized to unity. Usually, basis sets are given +a combination of normalized primitives, so the normalization +coefficients of the primitives, $f_{ks}$, need also to be provided. Atomic orbitals (AOs) are defined as \[ -\chi_i (\mathbf{r}) = P_{\eta(i)}(\mathbf{r})\, R_{\theta(i)} (\mathbf{r}) +\chi_i (\mathbf{r}) = \mathcal{M}_i\, P_{\eta(i)}(\mathbf{r})\, R_{\theta(i)} (\mathbf{r}) \] where $\theta(i)$ returns the shell on which the AO is expanded, and $\eta(i)$ denotes which angular function is chosen. +Here, the parameter $\mathcal{M}_i$ is an extra parameter which allows +the normalization of the different functions of the same shell to be +different, as in GAMESS for example. -In this section we describe the kernels used to compute the values, +In this section we describe first how the basis set is stored in the +context, and then we present the kernels used to compute the values, gradients and Laplacian of the atomic basis functions. * Headers :noexport: @@ -125,16 +128,23 @@ type = 'G' shell_num = 12 prim_num = 20 shell_center = [1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2] -shell_ang_mom = ['S', 'S', 'S', 'P', 'P', 'D', 'S', 'S', 'S', 'P', 'P', 'D'] +shell_ang_mom = [0, 0, 0, 1, 1, 2, 0, 0, 0, 1, 1, 2] shell_factor = [ 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.] shell_prim_num = [5, 1, 1, 1, 1, 1, 5, 1, 1, 1, 1, 1] shell_prim_index = [1, 6, 7, 8, 9, 10, 11, 16, 17, 18, 19, 20] -exponent = [ 33.87, 5.095, 1.159, 0.3258, 0.1027, 0.3258, 0.1027, - 1.407, 0.388, 1.057, 33.87, 5.095, 1.159, 0.3258, 0.1027, - 0.3258, 0.1027, 1.407, 0.388, 1.057] -coefficient = [ 0.006068, 0.045308, 0.202822, 0.503903, 0.383421, - 1.0, 1.0, 1.0, 1.0, 1.0, 0.006068, 0.045308, 0.202822, - 0.503903, 0.383421, 1.0, 1.0, 1.0, 1.0, 1.0] +exponent = [ 33.87, 5.095, 1.159, 0.3258, 0.1027, 0.3258, 0.1027, 1.407, + 0.388, 1.057, 33.87, 5.095, 1.159, 0.3258, 0.1027, 0.3258, 0.1027, 1.407, + 0.388, 1.057] +coefficient = [ 0.006068, 0.045308, 0.202822, 0.503903, 0.383421, 1.0, 1.0, + 1.0, 1.0, 1.0, 0.006068, 0.045308, 0.202822, 0.503903, 0.383421, 1.0, 1.0, + 1.0, 1.0, 1.0] +prim_factor = [ 1.0006253235944540e+01, 2.4169531573445120e+00, 7.9610924849766440e-01 + 3.0734305383061117e-01, 1.2929684417481876e-01, 3.0734305383061117e-01, + 1.2929684417481876e-01, 2.1842769845268308e+00, 4.3649547399719840e-01, + 1.8135965626177861e+00, 1.0006253235944540e+01, 2.4169531573445120e+00, + 7.9610924849766440e-01, 3.0734305383061117e-01, 1.2929684417481876e-01, + 3.0734305383061117e-01, 1.2929684417481876e-01, 2.1842769845268308e+00, + 4.3649547399719840e-01, 1.8135965626177861e+00 ] #+END_EXAMPLE ** Data structure @@ -151,6 +161,7 @@ typedef struct qmckl_ao_basis_struct { double * shell_factor; double * exponent ; double * coefficient ; + double * prim_factor ; bool provided; char type; } qmckl_ao_basis_struct; @@ -174,6 +185,7 @@ int64_t* qmckl_get_ao_basis_shell_prim_index (const qmckl_context context); double* qmckl_get_ao_basis_shell_factor (const qmckl_context context); double* qmckl_get_ao_basis_exponent (const qmckl_context context); double* qmckl_get_ao_basis_coefficient (const qmckl_context context); +double* qmckl_get_ao_basis_prim_factor (const qmckl_context context); #+end_src When all the data for the AOs have been provided, the following @@ -383,6 +395,25 @@ double* qmckl_get_ao_basis_coefficient (const qmckl_context context) { } +double* qmckl_get_ao_basis_prim_factor (const qmckl_context context) { + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return NULL; + } + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + int32_t mask = 1 << 8; + + if ( (ctx->ao_basis.uninitialized & mask) != 0) { + return NULL; + } + + assert (ctx->ao_basis.prim_factor != NULL); + return ctx->ao_basis.prim_factor; +} + + bool qmckl_ao_basis_provided(const qmckl_context context) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { @@ -412,6 +443,7 @@ qmckl_exit_code qmckl_set_ao_basis_shell_prim_num (qmckl_context context, con qmckl_exit_code qmckl_set_ao_basis_shell_factor (qmckl_context context, const double * shell_factor); qmckl_exit_code qmckl_set_ao_basis_exponent (qmckl_context context, const double * exponent); qmckl_exit_code qmckl_set_ao_basis_coefficient (qmckl_context context, const double * coefficient); +qmckl_exit_code qmckl_set_ao_basis_prim_factor (qmckl_context context, const double * prim_factor); #+end_src #+NAME:pre2 @@ -789,6 +821,47 @@ qmckl_exit_code qmckl_set_ao_basis_coefficient(qmckl_context context, const dou <> } +qmckl_exit_code qmckl_set_ao_basis_prim_factor(qmckl_context context, const double* prim_factor) { + <> + + int32_t mask = 1 << 7; + + const int64_t prim_num = qmckl_get_ao_basis_prim_num(context); + if (prim_num == 0L) { + return qmckl_failwith( context, + QMCKL_FAILURE, + "qmckl_set_ao_basis_prim_factor", + "prim_num is not set"); + } + + + if (ctx->ao_basis.prim_factor != NULL) { + qmckl_exit_code rc = qmckl_free(context, ctx->ao_basis.prim_factor); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_set_ao_basis_prim_factor", + NULL); + } + } + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = prim_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_ao_basis_prim_factor", + NULL); + } + + memcpy(new_array, prim_factor, mem_info.size); + + ctx->ao_basis.prim_factor = new_array; + + <> +} + #+end_src ** TODO Fortran interfaces diff --git a/org/qmckl_tests.org b/org/qmckl_tests.org index 77eedd4..304c9e4 100644 --- a/org/qmckl_tests.org +++ b/org/qmckl_tests.org @@ -523,6 +523,345 @@ F 1 #+END_example +#+begin_src c :tangle ../tests/chbrclf.h +#define chbrclf_shell_num 72 +#define chbrclf_prim_num 297 + +chbrclf_basis_shell_center[chbrclf_shell_num] = + {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, + 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5}; + +chbrclf_basis_shell_ang_mom[chbrclf_shell_num] = + {0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 3, 3, 0, 0, 0, 0, 1, 1, 1, 2, 2, 0, + 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 3, 3, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, + 2, 2, 2, 3, 3, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3}; + +chbrclf_basis_shell_prim_index[chbrclf_shell_num] = + {1, 11, 21, 22, 23, 24, 29, 30, 31, 32, 33, 34, 35, 36, 37, 42, 43, 44, 45, + 46, 47, 48, 49, 50, 60, 70, 71, 72, 73, 78, 79, 80, 81, 82, 83, 84, 85, 86, + 101, 116, 131, 132, 133, 134, 143, 152, 153, 154, 155, 156, 157, 158, 159, + 160, 180, 200, 220, 240, 241, 242, 243, 256, 269, 282, 283, 284, 285, 293, + 294, 295, 296, 297}; + +chbrclf_basis_shell_prim_num[chbrclf_shell_num] = + {10, 10, 1, 1, 1, 5, 1, 1, 1, 1, 1, 1, 1, 1, 5, 1, 1, 1, 1, 1, 1, 1, 1, 10, + 10, 1, 1, 1, 5, 1, 1, 1, 1, 1, 1, 1, 1, 15, 15, 15, 1, 1, 1, 9, 9, 1, 1, 1, + 1, 1, 1, 1, 1, 20, 20, 20, 20, 1, 1, 1, 13, 13, 13, 1, 1, 1, 8, 1, 1, 1, 1, + 1}; + +chbrclf_basis_shell_factor[chbrclf_shell_num] = + {1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., + 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., + 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., + 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.}; + +chbrclf_basis_exponent[chbrclf_prim_num] = + {8.2360000000000000e+03, 1.2350000000000000e+03, 2.8080000000000001e+02, + 7.9269999999999996e+01, 2.5590000000000000e+01, 8.9969999999999999e+00, + 3.3190000000000000e+00, 9.0590000000000004e-01, 3.6430000000000001e-01, + 1.2850000000000000e-01, 8.2360000000000000e+03, 1.2350000000000000e+03, + 2.8080000000000001e+02, 7.9269999999999996e+01, 2.5590000000000000e+01, + 8.9969999999999999e+00, 3.3190000000000000e+00, 9.0590000000000004e-01, + 3.6430000000000001e-01, 1.2850000000000000e-01, 9.0590000000000004e-01, + 1.2850000000000000e-01, 4.4019999999999997e-02, 1.8710000000000001e+01, + 4.1330000000000000e+00, 1.2000000000000000e+00, 3.8269999999999998e-01, + 1.2089999999999999e-01, 3.8269999999999998e-01, 1.2089999999999999e-01, + 3.5690000000000000e-02, 1.0970000000000000e+00, 3.1800000000000000e-01, + 1.0000000000000001e-01, 7.6100000000000001e-01, 2.6800000000000002e-01, + 3.3869999999999997e+01, 5.0949999999999998e+00, 1.1590000000000000e+00, + 3.2579999999999998e-01, 1.0270000000000000e-01, 3.2579999999999998e-01, + 1.0270000000000000e-01, 2.5260000000000001e-02, 1.4070000000000000e+00, + 3.8800000000000001e-01, 1.0199999999999999e-01, 1.0569999999999999e+00, + 2.4700000000000000e-01, 1.9500000000000000e+04, 2.9230000000000000e+03, + 6.6450000000000000e+02, 1.8750000000000000e+02, 6.0619999999999997e+01, + 2.1420000000000002e+01, 7.9500000000000002e+00, 2.2570000000000001e+00, + 8.8149999999999995e-01, 3.0409999999999998e-01, 1.9500000000000000e+04, + 2.9230000000000000e+03, 6.6450000000000000e+02, 1.8750000000000000e+02, + 6.0619999999999997e+01, 2.1420000000000002e+01, 7.9500000000000002e+00, + 2.2570000000000001e+00, 8.8149999999999995e-01, 3.0409999999999998e-01, + 2.2570000000000001e+00, 3.0409999999999998e-01, 9.1579999999999995e-02, + 4.3880000000000003e+01, 9.9260000000000002e+00, 2.9300000000000002e+00, + 9.1320000000000001e-01, 2.6719999999999999e-01, 9.1320000000000001e-01, + 2.6719999999999999e-01, 7.3609999999999995e-02, 3.1070000000000002e+00, + 8.5499999999999998e-01, 2.9199999999999998e-01, 1.9170000000000000e+00, + 7.2399999999999998e-01, 4.5610000000000000e+05, 6.8330000000000000e+04, + 1.5550000000000000e+04, 4.4050000000000000e+03, 1.4390000000000000e+03, + 5.2039999999999998e+02, 2.0309999999999999e+02, 8.3959999999999994e+01, + 3.6200000000000003e+01, 1.5830000000000000e+01, 6.3339999999999996e+00, + 2.6940000000000000e+00, 9.7680000000000000e-01, 4.3130000000000002e-01, + 1.6250000000000001e-01, 4.5610000000000000e+05, 6.8330000000000000e+04, + 1.5550000000000000e+04, 4.4050000000000000e+03, 1.4390000000000000e+03, + 5.2039999999999998e+02, 2.0309999999999999e+02, 8.3959999999999994e+01, + 3.6200000000000003e+01, 1.5830000000000000e+01, 6.3339999999999996e+00, + 2.6940000000000000e+00, 9.7680000000000000e-01, 4.3130000000000002e-01, + 1.6250000000000001e-01, 4.5610000000000000e+05, 6.8330000000000000e+04, + 1.5550000000000000e+04, 4.4050000000000000e+03, 1.4390000000000000e+03, + 5.2039999999999998e+02, 2.0309999999999999e+02, 8.3959999999999994e+01, + 3.6200000000000003e+01, 1.5830000000000000e+01, 6.3339999999999996e+00, + 2.6940000000000000e+00, 9.7680000000000000e-01, 4.3130000000000002e-01, + 1.6250000000000001e-01, 9.7680000000000000e-01, 1.6250000000000001e-01, + 5.9100000000000000e-02, 6.6329999999999995e+02, 1.5680000000000001e+02, + 4.9979999999999997e+01, 1.8420000000000002e+01, 7.2400000000000002e+00, + 2.9220000000000002e+00, 1.0220000000000000e+00, 3.8179999999999997e-01, + 1.3009999999999999e-01, 6.6329999999999995e+02, 1.5680000000000001e+02, + 4.9979999999999997e+01, 1.8420000000000002e+01, 7.2400000000000002e+00, + 2.9220000000000002e+00, 1.0220000000000000e+00, 3.8179999999999997e-01, + 1.3009999999999999e-01, 1.0220000000000000e+00, 1.3009999999999999e-01, + 4.1900000000000000e-02, 1.0460000000000000e+00, 3.4399999999999997e-01, + 1.3500000000000001e-01, 7.0599999999999996e-01, 3.1200000000000000e-01, + 1.0639000000000000e+07, 1.5934000000000000e+06, 3.6261000000000000e+05, + 1.0270000000000000e+05, 3.3501000000000000e+04, 1.2093000000000000e+04, + 4.7158999999999996e+03, 1.9555999999999999e+03, 8.5261000000000001e+02, + 3.8767000000000002e+02, 1.8268000000000001e+02, 8.8245000000000005e+01, + 3.9262999999999998e+01, 1.9234000000000002e+01, 9.4056999999999995e+00, + 4.1600999999999999e+00, 1.8995000000000000e+00, 6.0472000000000004e-01, + 3.0114000000000002e-01, 1.2515000000000001e-01, 1.0639000000000000e+07, + 1.5934000000000000e+06, 3.6261000000000000e+05, 1.0270000000000000e+05, + 3.3501000000000000e+04, 1.2093000000000000e+04, 4.7158999999999996e+03, + 1.9555999999999999e+03, 8.5261000000000001e+02, 3.8767000000000002e+02, + 1.8268000000000001e+02, 8.8245000000000005e+01, 3.9262999999999998e+01, + 1.9234000000000002e+01, 9.4056999999999995e+00, 4.1600999999999999e+00, + 1.8995000000000000e+00, 6.0472000000000004e-01, 3.0114000000000002e-01, + 1.2515000000000001e-01, 1.0639000000000000e+07, 1.5934000000000000e+06, + 3.6261000000000000e+05, 1.0270000000000000e+05, 3.3501000000000000e+04, + 1.2093000000000000e+04, 4.7158999999999996e+03, 1.9555999999999999e+03, + 8.5261000000000001e+02, 3.8767000000000002e+02, 1.8268000000000001e+02, + 8.8245000000000005e+01, 3.9262999999999998e+01, 1.9234000000000002e+01, + 9.4056999999999995e+00, 4.1600999999999999e+00, 1.8995000000000000e+00, + 6.0472000000000004e-01, 3.0114000000000002e-01, 1.2515000000000001e-01, + 1.0639000000000000e+07, 1.5934000000000000e+06, 3.6261000000000000e+05, + 1.0270000000000000e+05, 3.3501000000000000e+04, 1.2093000000000000e+04, + 4.7158999999999996e+03, 1.9555999999999999e+03, 8.5261000000000001e+02, + 3.8767000000000002e+02, 1.8268000000000001e+02, 8.8245000000000005e+01, + 3.9262999999999998e+01, 1.9234000000000002e+01, 9.4056999999999995e+00, + 4.1600999999999999e+00, 1.8995000000000000e+00, 6.0472000000000004e-01, + 3.0114000000000002e-01, 1.2515000000000001e-01, 6.0472000000000004e-01, + 1.2515000000000001e-01, 4.5593000000000002e-02, 8.6765000000000000e+03, + 2.0559000000000001e+03, 6.6623000000000002e+02, 2.5309999999999999e+02, + 1.0612000000000000e+02, 4.7241999999999997e+01, 2.1824999999999999e+01, + 9.9684000000000008e+00, 4.5171000000000001e+00, 1.9982000000000000e+00, + 7.0987999999999996e-01, 2.8144999999999998e-01, 1.0204000000000001e-01, + 8.6765000000000000e+03, 2.0559000000000001e+03, 6.6623000000000002e+02, + 2.5309999999999999e+02, 1.0612000000000000e+02, 4.7241999999999997e+01, + 2.1824999999999999e+01, 9.9684000000000008e+00, 4.5171000000000001e+00, + 1.9982000000000000e+00, 7.0987999999999996e-01, 2.8144999999999998e-01, + 1.0204000000000001e-01, 8.6765000000000000e+03, 2.0559000000000001e+03, + 6.6623000000000002e+02, 2.5309999999999999e+02, 1.0612000000000000e+02, + 4.7241999999999997e+01, 2.1824999999999999e+01, 9.9684000000000008e+00, + 4.5171000000000001e+00, 1.9982000000000000e+00, 7.0987999999999996e-01, + 2.8144999999999998e-01, 1.0204000000000001e-01, 7.0987999999999996e-01, + 1.0204000000000001e-01, 3.5142000000000000e-02, 4.0382999999999998e+02, + 1.2117000000000000e+02, 4.6344999999999999e+01, 1.9721000000000000e+01, + 8.8623999999999992e+00, 3.9962000000000000e+00, 1.7636000000000001e+00, + 7.0618999999999998e-01, 7.0618999999999998e-01, 2.6390000000000002e-01, + 1.0470000000000000e-01, 5.5149999999999999e-01, 2.5800000000000001e-01}; + +chbrclf_basis_coefficient[chbrclf_prim_num] = + {5.3100000000000000e-04, 4.1079999999999997e-03, 2.1087000000000002e-02, + 8.1852999999999995e-02, 2.3481700000000000e-01, 4.3440099999999998e-01, + 3.4612900000000002e-01, 3.9378000000000003e-02, -8.9829999999999997e-03, + 2.3850000000000000e-03, -1.1300000000000000e-04, -8.7799999999999998e-04, + -4.5399999999999998e-03, -1.8133000000000000e-02, -5.5759999999999997e-02, + -1.2689500000000001e-01, -1.7035200000000000e-01, 1.4038200000000001e-01, + 5.9868399999999999e-01, 3.9538899999999999e-01, 1.0000000000000000e+00, + 1.0000000000000000e+00, 1.0000000000000000e+00, 1.4031000000000000e-02, + 8.6865999999999999e-02, 2.9021599999999997e-01, 5.0100800000000001e-01, + 3.4340599999999999e-01, 1.0000000000000000e+00, 1.0000000000000000e+00, + 1.0000000000000000e+00, 1.0000000000000000e+00, 1.0000000000000000e+00, + 1.0000000000000000e+00, 1.0000000000000000e+00, 1.0000000000000000e+00, + 6.0679999999999996e-03, 4.5308000000000001e-02, 2.0282200000000000e-01, + 5.0390299999999999e-01, 3.8342100000000001e-01, 1.0000000000000000e+00, + 1.0000000000000000e+00, 1.0000000000000000e+00, 1.0000000000000000e+00, + 1.0000000000000000e+00, 1.0000000000000000e+00, 1.0000000000000000e+00, + 1.0000000000000000e+00, 5.0699999999999996e-04, 3.9230000000000003e-03, + 2.0199999999999999e-02, 7.9009999999999997e-02, 2.3043900000000000e-01, + 4.3287199999999998e-01, 3.4996400000000000e-01, 4.3233000000000001e-02, + -7.8919999999999997e-03, 2.3839999999999998e-03, -1.1700000000000000e-04, + -9.1200000000000005e-04, -4.7169999999999998e-03, -1.9085999999999999e-02, + -5.9655000000000000e-02, -1.4001000000000000e-01, -1.7678199999999999e-01, + 1.7162500000000000e-01, 6.0504300000000000e-01, 3.6951200000000001e-01, + 1.0000000000000000e+00, 1.0000000000000000e+00, 1.0000000000000000e+00, + 1.6664999999999999e-02, 1.0447200000000000e-01, 3.1725999999999999e-01, + 4.8734300000000003e-01, 3.3460400000000001e-01, 1.0000000000000000e+00, + 1.0000000000000000e+00, 1.0000000000000000e+00, 1.0000000000000000e+00, + 1.0000000000000000e+00, 1.0000000000000000e+00, 1.0000000000000000e+00, + 1.0000000000000000e+00, 4.9296999999999999e-05, 3.8302900000000001e-04, + 2.0085400000000001e-03, 8.3855800000000001e-03, 2.9470300000000001e-02, + 8.7832499999999994e-02, 2.1147299999999999e-01, 3.6536400000000002e-01, + 3.4088400000000002e-01, 1.0213300000000000e-01, 3.1167500000000002e-03, + 1.0575100000000000e-03, -3.7800000000000003e-04, 1.5613600000000000e-04, + -5.1412600000000003e-05, -1.3830400000000001e-05, -1.0727900000000000e-04, + -5.6508299999999997e-04, -2.3613499999999999e-03, -8.4588600000000003e-03, + -2.5963799999999999e-02, -6.8636199999999994e-02, -1.4187400000000000e-01, + -1.9931900000000000e-01, -1.9566199999999999e-02, 4.9974099999999999e-01, + 5.6373600000000001e-01, 7.9032500000000006e-02, -8.3509099999999996e-03, + 2.3245599999999998e-03, 4.1854599999999997e-06, 3.2439500000000000e-05, + 1.7110500000000001e-04, 7.1417599999999996e-04, 2.5670500000000000e-03, + 7.8855200000000000e-03, 2.1086700000000000e-02, 4.4226399999999999e-02, + 6.5167000000000003e-02, 6.0301199999999999e-03, -2.0649500000000001e-01, + -4.0587099999999998e-01, 7.5955800000000004e-02, 7.2566100000000000e-01, + 3.9442300000000002e-01, 1.0000000000000000e+00, 1.0000000000000000e+00, + 1.0000000000000000e+00, 2.4044800000000001e-03, 1.9214800000000001e-02, + 8.8509699999999997e-02, 2.5602000000000003e-01, 4.3692700000000001e-01, + 3.5033399999999998e-01, 5.8549499999999997e-02, -4.5842299999999999e-03, + 2.2696999999999999e-03, -6.5214500000000003e-04, -5.1944499999999998e-03, + -2.4693799999999998e-02, -7.2816699999999998e-02, -1.3403000000000001e-01, + -9.4774200000000003e-02, 2.6228899999999999e-01, 5.6466700000000003e-01, + 3.4125000000000000e-01, 1.0000000000000000e+00, 1.0000000000000000e+00, + 1.0000000000000000e+00, 1.0000000000000000e+00, 1.0000000000000000e+00, + 1.0000000000000000e+00, 1.0000000000000000e+00, 1.0000000000000000e+00, + 5.9000000000000003e-06, 4.6100000000000002e-05, 2.4220000000000001e-04, + 1.0226000000000000e-03, 3.7112999999999998e-03, 1.1978500000000000e-02, + 3.4692700000000000e-02, 8.9123900000000006e-02, 1.9345570000000001e-01, + 3.2090190000000002e-01, 3.2992329999999997e-01, 1.4941209999999999e-01, + 1.4993800000000000e-02, -9.1650000000000000e-04, 4.3800000000000002e-04, + -2.3980000000000000e-04, 7.3600000000000000e-05, -3.6699999999999998e-05, + 2.3900000000000002e-05, -5.5999999999999997e-06, -1.9000000000000000e-06, + -1.4500000000000000e-05, -7.6100000000000007e-05, -3.2100000000000000e-04, + -1.1708999999999999e-03, -3.7967999999999999e-03, -1.1230700000000000e-02, + -2.9927700000000002e-02, -7.1270600000000003e-02, -1.4031360000000001e-01, + -2.0307629999999999e-01, -9.6098500000000003e-02, 3.5580859999999997e-01, + 5.9217920000000002e-01, 2.2159770000000001e-01, 1.3764800000000001e-02, + 8.3949999999999997e-04, -4.5099999999999998e-05, -8.4999999999999999e-06, + -1.2400000000000000e-05, 6.9999999999999997e-07, 5.6999999999999996e-06, + 3.0300000000000001e-05, 1.2750000000000001e-04, 4.6589999999999999e-04, + 1.5096000000000000e-03, 4.4852000000000000e-03, 1.1983499999999999e-02, + 2.8957100000000000e-02, 5.8156600000000003e-02, 8.8813299999999998e-02, + 4.4524399999999999e-02, -2.0603869999999999e-01, -5.1270170000000004e-01, + -1.5093490000000001e-01, 6.7892030000000003e-01, 5.8176969999999995e-01, + 4.6755499999999998e-02, -1.1182500000000000e-02, 2.4402000000000000e-03, + -1.9999999999999999e-07, -1.7999999999999999e-06, -9.3000000000000007e-06, + -3.9100000000000002e-05, -1.4280000000000000e-04, -4.6279999999999997e-04, + -1.3749999999999999e-03, -3.6784000000000001e-03, -8.8981000000000008e-03, + -1.7952900000000001e-02, -2.7573199999999999e-02, -1.4095300000000000e-02, + 6.7256099999999999e-02, 1.7669280000000001e-01, 5.2886099999999998e-02, + -3.0759550000000002e-01, -4.7006579999999998e-01, 2.5587610000000000e-01, + 6.9803409999999999e-01, 2.9672559999999998e-01, 1.0000000000000000e+00, + 1.0000000000000000e+00, 1.0000000000000000e+00, 4.3570000000000002e-04, + 3.7815000000000001e-03, 2.0478199999999998e-02, 7.9283400000000004e-02, + 2.1784729999999999e-01, 3.8785849999999999e-01, 3.5943500000000000e-01, + 1.1219949999999999e-01, 4.3873999999999996e-03, 1.7809000000000000e-03, + -4.5760000000000001e-04, 2.1220000000000001e-04, -7.3399999999999995e-05, + -1.7479999999999999e-04, -1.5263000000000000e-03, -8.3399000000000008e-03, + -3.3220300000000001e-02, -9.5418000000000003e-02, -1.8240260000000000e-01, + -1.5583079999999999e-01, 1.8678990000000001e-01, 5.4277330000000001e-01, + 3.8733089999999998e-01, 4.5306899999999997e-02, -4.3784000000000002e-03, + 1.8110999999999999e-03, 4.5099999999999998e-05, 3.9639999999999999e-04, + 2.1554999999999999e-03, 8.6719999999999992e-03, 2.4868000000000001e-02, + 4.8547199999999999e-02, 3.9615600000000001e-02, -6.0574900000000001e-02, + -1.8716990000000000e-01, -1.3777570000000000e-01, 2.9280210000000001e-01, + 5.7608959999999998e-01, 3.0786170000000002e-01, 1.0000000000000000e+00, + 1.0000000000000000e+00, 1.0000000000000000e+00, 1.4732000000000000e-03, + 1.2672500000000000e-02, 5.8045100000000002e-02, 1.7051030000000000e-01, + 3.1859579999999998e-01, 3.8450230000000002e-01, 2.7377370000000001e-01, + 7.4396699999999996e-02, 1.0000000000000000e+00, 1.0000000000000000e+00, + 1.0000000000000000e+00, 1.0000000000000000e+00, 1.0000000000000000e+00}; + +chbrclf_basis_prim_factor[chbrclf_prim_num] = + {6.1616545431994848e+02, 1.4847738511079908e+02, 4.8888635917437597e+01, + 1.8933972232608955e+01, 8.1089160941724145e+00, 3.7024003863155635e+00, + 1.7525302846177560e+00, 6.6179013183966806e-01, 3.3419848027174592e-01, + 1.5296336817449557e-01, 6.1616545431994848e+02, 1.4847738511079908e+02, + 4.8888635917437597e+01, 1.8933972232608955e+01, 8.1089160941724145e+00, + 3.7024003863155635e+00, 1.7525302846177560e+00, 6.6179013183966806e-01, + 3.3419848027174592e-01, 1.5296336817449557e-01, 6.6179013183966806e-01, + 1.5296336817449557e-01, 6.8493225861981921e-02, 5.5466699238441954e+01, + 8.3998560685400019e+00, 1.7902622143452276e+00, 4.2905519588435126e-01, + 1.0161854305479753e-01, 4.2905519588435126e-01, 1.0161854305479753e-01, + 2.2111758010684022e-02, 1.9354014159719681e+00, 2.2164447815916102e-01, + 2.9269105913429974e-02, 7.9628755341813429e-01, 7.6077763383300537e-02, + 1.0006253235944540e+01, 2.4169531573445120e+00, 7.9610924849766440e-01, + 3.0734305383061117e-01, 1.2929684417481876e-01, 3.0734305383061117e-01, + 1.2929684417481876e-01, 4.5158041868216925e-02, 2.1842769845268308e+00, + 4.3649547399719840e-01, 8.2165651391863506e-02, 1.8135965626177861e+00, + 1.4243906834168285e-01, 1.1760777961352585e+03, 2.8332291650568584e+02, + 9.3278452222064189e+01, 3.6112790320330610e+01, 1.5483603491420400e+01, + 7.0961811262237955e+00, 3.3743135718917610e+00, 1.3123774561862465e+00, + 6.4837584903810197e-01, 2.9185854115641796e-01, 1.1760777961352585e+03, + 2.8332291650568584e+02, 9.3278452222064189e+01, 3.6112790320330610e+01, + 1.5483603491420400e+01, 7.0961811262237955e+00, 3.3743135718917610e+00, + 1.3123774561862465e+00, 6.4837584903810197e-01, 2.9185854115641796e-01, + 1.3123774561862465e+00, 2.9185854115641796e-01, 1.1864804090515012e-01, + 1.6098053319659394e+02, 2.5113539500925100e+01, 5.4641751683240054e+00, + 1.2724697488890255e+00, 2.7383291850797253e-01, 1.2724697488890255e+00, + 2.7383291850797253e-01, 5.4652635549581594e-02, 1.1967588544907814e+01, + 1.2512674783228661e+00, 1.9091033786232964e-01, 6.3658177914201666e+00, + 7.1181350527019893e-01, 1.2508497509090121e+04, 3.0120949486045815e+03, + 9.9244828879919498e+02, 3.8536256746518035e+02, 1.6651589181438584e+02, + 7.7653861740552273e+01, 3.8343569039456206e+01, 1.9768076572700657e+01, + 1.0518202315565855e+01, 5.6561481770276556e+00, 2.8455658258813377e+00, + 1.4986778401686554e+00, 7.0026807980299932e-01, 3.7931034358525295e-01, + 1.8241061862759339e-01, 1.2508497509090121e+04, 3.0120949486045815e+03, + 9.9244828879919498e+02, 3.8536256746518035e+02, 1.6651589181438584e+02, + 7.7653861740552273e+01, 3.8343569039456206e+01, 1.9768076572700657e+01, + 1.0518202315565855e+01, 5.6561481770276556e+00, 2.8455658258813377e+00, + 1.4986778401686554e+00, 7.0026807980299932e-01, 3.7931034358525295e-01, + 1.8241061862759339e-01, 1.2508497509090121e+04, 3.0120949486045815e+03, + 9.9244828879919498e+02, 3.8536256746518035e+02, 1.6651589181438584e+02, + 7.7653861740552273e+01, 3.8343569039456206e+01, 1.9768076572700657e+01, + 1.0518202315565855e+01, 5.6561481770276556e+00, 2.8455658258813377e+00, + 1.4986778401686554e+00, 7.0026807980299932e-01, 3.7931034358525295e-01, + 1.8241061862759339e-01, 7.0026807980299932e-01, 1.8241061862759339e-01, + 8.5428091252337218e-02, 4.7981915824835833e+03, 7.9090197643220097e+02, + 1.8942417420993877e+02, 5.4394140688380837e+01, 1.6928298001679121e+01, + 5.4455324755334713e+00, 1.4647169510384077e+00, 4.2779429930932966e-01, + 1.1137474138247395e-01, 4.7981915824835833e+03, 7.9090197643220097e+02, + 1.8942417420993877e+02, 5.4394140688380837e+01, 1.6928298001679121e+01, + 5.4455324755334713e+00, 1.4647169510384077e+00, 4.2779429930932966e-01, + 1.1137474138247395e-01, 1.4647169510384077e+00, 1.1137474138247395e-01, + 2.7021385701525968e-02, 1.7806964960637739e+00, 2.5432363995130330e-01, + 4.9487276238674341e-02, 6.7261264398159915e-01, 1.0710316041250582e-01, + 1.3276564169936487e+05, 3.1963453883259335e+04, 1.0531492549479324e+04, + 4.0887252186430137e+03, 1.7648322054568350e+03, 8.2188420053589937e+02, + 4.0558653897685440e+02, 2.0958978573126603e+02, 1.1245348234327012e+02, + 6.2266827352587455e+01, 3.5414274108483511e+01, 2.0520014948183928e+01, + 1.1178864762092555e+01, 6.5457935075374944e+00, 3.8278369835025656e+00, + 2.0760505305870112e+00, 1.1531599898261422e+00, 4.8873770086696849e-01, + 2.8972530369835303e-01, 1.4996269365355971e-01, 1.3276564169936487e+05, + 3.1963453883259335e+04, 1.0531492549479324e+04, 4.0887252186430137e+03, + 1.7648322054568350e+03, 8.2188420053589937e+02, 4.0558653897685440e+02, + 2.0958978573126603e+02, 1.1245348234327012e+02, 6.2266827352587455e+01, + 3.5414274108483511e+01, 2.0520014948183928e+01, 1.1178864762092555e+01, + 6.5457935075374944e+00, 3.8278369835025656e+00, 2.0760505305870112e+00, + 1.1531599898261422e+00, 4.8873770086696849e-01, 2.8972530369835303e-01, + 1.4996269365355971e-01, 1.3276564169936487e+05, 3.1963453883259335e+04, + 1.0531492549479324e+04, 4.0887252186430137e+03, 1.7648322054568350e+03, + 8.2188420053589937e+02, 4.0558653897685440e+02, 2.0958978573126603e+02, + 1.1245348234327012e+02, 6.2266827352587455e+01, 3.5414274108483511e+01, + 2.0520014948183928e+01, 1.1178864762092555e+01, 6.5457935075374944e+00, + 3.8278369835025656e+00, 2.0760505305870112e+00, 1.1531599898261422e+00, + 4.8873770086696849e-01, 2.8972530369835303e-01, 1.4996269365355971e-01, + 1.3276564169936487e+05, 3.1963453883259335e+04, 1.0531492549479324e+04, + 4.0887252186430137e+03, 1.7648322054568350e+03, 8.2188420053589937e+02, + 4.0558653897685440e+02, 2.0958978573126603e+02, 1.1245348234327012e+02, + 6.2266827352587455e+01, 3.5414274108483511e+01, 2.0520014948183928e+01, + 1.1178864762092555e+01, 6.5457935075374944e+00, 3.8278369835025656e+00, + 2.0760505305870112e+00, 1.1531599898261422e+00, 4.8873770086696849e-01, + 2.8972530369835303e-01, 1.4996269365355971e-01, 4.8873770086696849e-01, + 1.4996269365355971e-01, 7.0320786489653203e-02, 1.1936329579215313e+05, + 1.9732975244933248e+04, 4.8247000542937221e+03, 1.4389816948051262e+03, + 4.8549709936386239e+02, 1.7654297142185436e+02, 6.7240804881705529e+01, + 2.5247705079657806e+01, 9.3867385006594475e+00, 3.3864040992879496e+00, + 9.2879798315626561e-01, 2.9220769881703862e-01, 8.2205930646140513e-02, + 1.1936329579215313e+05, 1.9732975244933248e+04, 4.8247000542937221e+03, + 1.4389816948051262e+03, 4.8549709936386239e+02, 1.7654297142185436e+02, + 6.7240804881705529e+01, 2.5247705079657806e+01, 9.3867385006594475e+00, + 3.3864040992879496e+00, 9.2879798315626561e-01, 2.9220769881703862e-01, + 8.2205930646140513e-02, 1.1936329579215313e+05, 1.9732975244933248e+04, + 4.8247000542937221e+03, 1.4389816948051262e+03, 4.8549709936386239e+02, + 1.7654297142185436e+02, 6.7240804881705529e+01, 2.5247705079657806e+01, + 9.3867385006594475e+00, 3.3864040992879496e+00, 9.2879798315626561e-01, + 2.9220769881703862e-01, 8.2205930646140513e-02, 9.2879798315626561e-01, + 8.2205930646140513e-02, 2.1688183591227813e-02, 5.9876577632594533e+04, + 7.2836806319891484e+03, 1.3549226646722386e+03, 3.0376315094739988e+02, + 7.4924579607137730e+01, 1.8590543353806009e+01, 4.4423176930919421e+00, + 8.9541051939952665e-01, 8.9541051939952665e-01, 1.5992942988584680e-01, + 3.1718756222897104e-02, 3.8586186799894789e-01, 6.9839124768946298e-02}; + +#+end_src + +** TODO Molecular orbitals ** Electron coordinates Electron coordinates are stored in atomic units in normal format. From 8090edc771bef34339d84ed6e849df9401d885be Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 3 Jun 2021 22:34:23 +0200 Subject: [PATCH 5/9] Fix AO tests --- org/qmckl_ao.org | 136 ++++++++++++++++++----------------------- org/qmckl_context.org | 2 +- org/qmckl_electron.org | 16 ++--- org/qmckl_tests.org | 26 ++++---- 4 files changed, 83 insertions(+), 97 deletions(-) diff --git a/org/qmckl_ao.org b/org/qmckl_ao.org index 2e40ada..64cb745 100644 --- a/org/qmckl_ao.org +++ b/org/qmckl_ao.org @@ -56,6 +56,8 @@ gradients and Laplacian of the atomic basis functions. #include "config.h" #endif +#include "chbrclf.h" + int main() { qmckl_context context; context = qmckl_context_create(); @@ -155,7 +157,7 @@ typedef struct qmckl_ao_basis_struct { int64_t shell_num; int64_t prim_num; int64_t * shell_center; - char * shell_ang_mom; + int32_t * shell_ang_mom; int64_t * shell_prim_num; int64_t * shell_prim_index; double * shell_factor; @@ -179,7 +181,7 @@ char qmckl_get_ao_basis_type (const qmckl_context context); int64_t qmckl_get_ao_basis_shell_num (const qmckl_context context); int64_t qmckl_get_ao_basis_prim_num (const qmckl_context context); int64_t* qmckl_get_ao_basis_shell_center (const qmckl_context context); -char* qmckl_get_ao_basis_shell_ang_mom (const qmckl_context context); +int32_t* qmckl_get_ao_basis_shell_ang_mom (const qmckl_context context); int64_t* qmckl_get_ao_basis_shell_prim_num (const qmckl_context context); int64_t* qmckl_get_ao_basis_shell_prim_index (const qmckl_context context); double* qmckl_get_ao_basis_shell_factor (const qmckl_context context); @@ -280,7 +282,7 @@ int64_t* qmckl_get_ao_basis_shell_center (const qmckl_context context) { } -char* qmckl_get_ao_basis_shell_ang_mom (const qmckl_context context) { +int32_t* qmckl_get_ao_basis_shell_ang_mom (const qmckl_context context) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return NULL; } @@ -403,7 +405,7 @@ double* qmckl_get_ao_basis_prim_factor (const qmckl_context context) { qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - int32_t mask = 1 << 8; + int32_t mask = 1 << 10; if ( (ctx->ao_basis.uninitialized & mask) != 0) { return NULL; @@ -438,7 +440,7 @@ qmckl_exit_code qmckl_set_ao_basis_shell_num (qmckl_context context, con qmckl_exit_code qmckl_set_ao_basis_prim_num (qmckl_context context, const int64_t prim_num); qmckl_exit_code qmckl_set_ao_basis_shell_prim_index (qmckl_context context, const int64_t * shell_prim_index); qmckl_exit_code qmckl_set_ao_basis_shell_center (qmckl_context context, const int64_t * shell_center); -qmckl_exit_code qmckl_set_ao_basis_shell_ang_mom (qmckl_context context, const char * shell_ang_mom); +qmckl_exit_code qmckl_set_ao_basis_shell_ang_mom (qmckl_context context, const int32_t * shell_ang_mom); qmckl_exit_code qmckl_set_ao_basis_shell_prim_num (qmckl_context context, const int64_t * shell_prim_num); qmckl_exit_code qmckl_set_ao_basis_shell_factor (qmckl_context context, const double * shell_factor); qmckl_exit_code qmckl_set_ao_basis_exponent (qmckl_context context, const double * exponent); @@ -463,7 +465,6 @@ ctx->ao_basis.provided = (ctx->ao_basis.uninitialized == 0); return QMCKL_SUCCESS; #+end_src - #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_set_ao_basis_type(qmckl_context context, const char t) { <> @@ -575,7 +576,7 @@ qmckl_exit_code qmckl_set_ao_basis_shell_center(qmckl_context context, const in } -qmckl_exit_code qmckl_set_ao_basis_shell_ang_mom(qmckl_context context, const char* shell_ang_mom) { +qmckl_exit_code qmckl_set_ao_basis_shell_ang_mom(qmckl_context context, const int32_t* shell_ang_mom) { <> int32_t mask = 1 << 4; @@ -600,7 +601,7 @@ qmckl_exit_code qmckl_set_ao_basis_shell_ang_mom(qmckl_context context, const c qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = shell_num * sizeof(char); - char* new_array = (char*) qmckl_malloc(context, mem_info); + int32_t * new_array = (int32_t*) qmckl_malloc(context, mem_info); if (new_array == NULL) { return qmckl_failwith( context, @@ -824,7 +825,7 @@ qmckl_exit_code qmckl_set_ao_basis_coefficient(qmckl_context context, const dou qmckl_exit_code qmckl_set_ao_basis_prim_factor(qmckl_context context, const double* prim_factor) { <> - int32_t mask = 1 << 7; + int32_t mask = 1 << 10; const int64_t prim_num = qmckl_get_ao_basis_prim_num(context); if (prim_num == 0L) { @@ -868,88 +869,73 @@ qmckl_exit_code qmckl_set_ao_basis_prim_factor(qmckl_context context, const dou ** Test :noexport: - #+begin_src c :tangle (eval c_test) :exports none :exports none - /* Reference input data */ + #+begin_src c :tangle (eval c_test) :exports none :exports none +const int64_t shell_num = chbrclf_shell_num; +const int64_t prim_num = chbrclf_prim_num; +const int64_t * shell_center = &(chbrclf_basis_shell_center[0]); +const int32_t * shell_ang_mom = &(chbrclf_basis_shell_ang_mom[0]); +const int64_t * shell_prim_num = &(chbrclf_basis_shell_prim_num[0]); +const int64_t * shell_prim_index = &(chbrclf_basis_shell_prim_index[0]); +const double * shell_factor = &(chbrclf_basis_shell_factor[0]); +const double * exponent = &(chbrclf_basis_exponent[0]); +const double * coefficient = &(chbrclf_basis_coefficient[0]); +const double * prim_factor = &(chbrclf_basis_prim_factor[0]); - char typ = 'G'; - #define shell_num ((int64_t) 12) - #define prim_num ((int64_t) 20) +char typ = 'G'; - int64_t shell_center [shell_num] = - { 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2 }; +qmckl_exit_code rc; - char shell_ang_mom [shell_num] = - { 'S', 'S', 'S', 'P', 'P', 'D', 'S', 'S', 'S', 'P', 'P', 'D' }; +assert(!qmckl_ao_basis_provided(context)); - double shell_factor [shell_num] = - { 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1. }; +rc = qmckl_set_ao_basis_type (context, typ); +assert(rc == QMCKL_SUCCESS); +assert(!qmckl_ao_basis_provided(context)); - int64_t shell_prim_num [shell_num] = - {5, 1, 1, 1, 1, 1, 5, 1, 1, 1, 1, 1}; +rc = qmckl_set_ao_basis_shell_num (context, shell_num); +assert(rc == QMCKL_SUCCESS); +assert(!qmckl_ao_basis_provided(context)); - int64_t shell_prim_index [shell_num] = - {1, 6, 7, 8, 9, 10, 11, 16, 17, 18, 19, 20}; +rc = qmckl_set_ao_basis_prim_num (context, prim_num); +assert(rc == QMCKL_SUCCESS); +assert(!qmckl_ao_basis_provided(context)); - double exponent [prim_num] = - { 33.87, 5.095, 1.159, 0.3258, 0.1027, 0.3258, 0.1027, - 1.407, 0.388, 1.057, 33.87, 5.095, 1.159, 0.3258, 0.1027, - 0.3258, 0.1027, 1.407, 0.388, 1.057 }; +rc = qmckl_set_ao_basis_shell_center (context, shell_center); +assert(rc == QMCKL_SUCCESS); +assert(!qmckl_ao_basis_provided(context)); - double coefficient [prim_num] = - { 0.006068, 0.045308, 0.202822, 0.503903, 0.383421, - 1.0, 1.0, 1.0, 1.0, 1.0, 0.006068, 0.045308, 0.202822, - 0.503903, 0.383421, 1.0, 1.0, 1.0, 1.0, 1.0 }; - /* --- */ +rc = qmckl_set_ao_basis_shell_ang_mom (context, shell_ang_mom); +assert(rc == QMCKL_SUCCESS); +assert(!qmckl_ao_basis_provided(context)); - qmckl_exit_code rc; +rc = qmckl_set_ao_basis_shell_factor (context, shell_factor); +assert(rc == QMCKL_SUCCESS); +assert(!qmckl_ao_basis_provided(context)); - assert(!qmckl_ao_basis_provided(context)); +rc = qmckl_set_ao_basis_shell_center (context, shell_prim_num); +assert(rc == QMCKL_SUCCESS); +assert(!qmckl_ao_basis_provided(context)); - rc = qmckl_set_ao_basis_type (context, typ); - assert(rc == QMCKL_SUCCESS); - assert(!qmckl_ao_basis_provided(context)); +rc = qmckl_set_ao_basis_shell_prim_num (context, shell_prim_num); +assert(rc == QMCKL_SUCCESS); +assert(!qmckl_ao_basis_provided(context)); - rc = qmckl_set_ao_basis_shell_num (context, shell_num); - assert(rc == QMCKL_SUCCESS); - assert(!qmckl_ao_basis_provided(context)); +rc = qmckl_set_ao_basis_shell_prim_index (context, shell_prim_index); +assert(rc == QMCKL_SUCCESS); +assert(!qmckl_ao_basis_provided(context)); - rc = qmckl_set_ao_basis_prim_num (context, prim_num); - assert(rc == QMCKL_SUCCESS); - assert(!qmckl_ao_basis_provided(context)); +rc = qmckl_set_ao_basis_exponent (context, exponent); +assert(rc == QMCKL_SUCCESS); +assert(!qmckl_ao_basis_provided(context)); - rc = qmckl_set_ao_basis_shell_center (context, shell_center); - assert(rc == QMCKL_SUCCESS); - assert(!qmckl_ao_basis_provided(context)); +rc = qmckl_set_ao_basis_coefficient (context, coefficient); +assert(rc == QMCKL_SUCCESS); +assert(!qmckl_ao_basis_provided(context)); - rc = qmckl_set_ao_basis_shell_ang_mom (context, shell_ang_mom); - assert(rc == QMCKL_SUCCESS); - assert(!qmckl_ao_basis_provided(context)); +rc = qmckl_set_ao_basis_prim_factor (context, prim_factor); +assert(rc == QMCKL_SUCCESS); +assert(qmckl_ao_basis_provided(context)); - rc = qmckl_set_ao_basis_shell_factor (context, shell_factor); - assert(rc == QMCKL_SUCCESS); - assert(!qmckl_ao_basis_provided(context)); - - rc = qmckl_set_ao_basis_shell_center (context, shell_prim_num); - assert(rc == QMCKL_SUCCESS); - assert(!qmckl_ao_basis_provided(context)); - - rc = qmckl_set_ao_basis_shell_prim_num (context, shell_prim_num); - assert(rc == QMCKL_SUCCESS); - assert(!qmckl_ao_basis_provided(context)); - - rc = qmckl_set_ao_basis_shell_prim_index (context, shell_prim_index); - assert(rc == QMCKL_SUCCESS); - assert(!qmckl_ao_basis_provided(context)); - - rc = qmckl_set_ao_basis_exponent (context, exponent); - assert(rc == QMCKL_SUCCESS); - assert(!qmckl_ao_basis_provided(context)); - - rc = qmckl_set_ao_basis_coefficient (context, coefficient); - assert(rc == QMCKL_SUCCESS); - assert(qmckl_ao_basis_provided(context)); - - #+end_src + #+end_src * Polynomial part ** Powers of $x-X_i$ diff --git a/org/qmckl_context.org b/org/qmckl_context.org index 6472365..618120d 100644 --- a/org/qmckl_context.org +++ b/org/qmckl_context.org @@ -221,7 +221,7 @@ qmckl_context qmckl_context_create() { ctx->numprec.precision = QMCKL_DEFAULT_PRECISION; ctx->numprec.range = QMCKL_DEFAULT_RANGE; - ctx->ao_basis.uninitialized = (1 << 10) - 1; + ctx->ao_basis.uninitialized = (1 << 11) - 1; ctx->nucleus.uninitialized = (1 << 4) - 1; ctx->electron.uninitialized = (1 << 3) - 1; diff --git a/org/qmckl_electron.org b/org/qmckl_electron.org index 0f47ff9..bc56cda 100644 --- a/org/qmckl_electron.org +++ b/org/qmckl_electron.org @@ -1615,7 +1615,7 @@ qmckl_exit_code qmckl_provide_en_distance_rescaled(qmckl_context context) | qmckl_context | context | in | Global state | | int64_t | elec_num | in | Number of electrons | | int64_t | nucl_num | in | Number of nuclei | - | double | rescale_factor_kappa_en | in | The factor for rescaled distances | + | double | rescale_factor_kappa_en | in | The factor for rescaled distances | | int64_t | walk_num | in | Number of walkers | | double | elec_coord[walk_num][3][elec_num] | in | Electron coordinates | | double | nucl_coord[3][elec_num] | in | Nuclear coordinates | @@ -1766,28 +1766,28 @@ assert(qmckl_nucleus_provided(context)); double en_distance_rescaled[walk_num][nucl_num][elec_num]; -rc = qmckl_get_electron_en_distance_rescaled(context, &(en_distance[0][0][0])); +rc = qmckl_get_electron_en_distance_rescaled(context, &(en_distance_rescaled[0][0][0])); assert (rc == QMCKL_SUCCESS); // TODO: check exact values //// (e,n,w) in Fortran notation //// (1,1,1) -//assert(fabs(en_distance[0][0][0] - 7.546738741619978) < 1.e-12); +//assert(fabs(en_distance_rescaled[0][0][0] - 7.546738741619978) < 1.e-12); // //// (1,2,1) -//assert(fabs(en_distance[0][1][0] - 8.77102435246984) < 1.e-12); +//assert(fabs(en_distance_rescaled[0][1][0] - 8.77102435246984) < 1.e-12); // //// (2,1,1) -//assert(fabs(en_distance[0][0][1] - 3.698922010513608) < 1.e-12); +//assert(fabs(en_distance_rescaled[0][0][1] - 3.698922010513608) < 1.e-12); // //// (1,1,2) -//assert(fabs(en_distance[1][0][0] - 5.824059436060509) < 1.e-12); +//assert(fabs(en_distance_rescaled[1][0][0] - 5.824059436060509) < 1.e-12); // //// (1,2,2) -//assert(fabs(en_distance[1][1][0] - 7.080482110317645) < 1.e-12); +//assert(fabs(en_distance_rescaled[1][1][0] - 7.080482110317645) < 1.e-12); // //// (2,1,2) -//assert(fabs(en_distance[1][0][1] - 3.1804527583077356) < 1.e-12); +//assert(fabs(en_distance_rescaled[1][0][1] - 3.1804527583077356) < 1.e-12); #+end_src diff --git a/org/qmckl_tests.org b/org/qmckl_tests.org index 304c9e4..3c090d1 100644 --- a/org/qmckl_tests.org +++ b/org/qmckl_tests.org @@ -527,36 +527,36 @@ F 1 #define chbrclf_shell_num 72 #define chbrclf_prim_num 297 -chbrclf_basis_shell_center[chbrclf_shell_num] = +int64_t chbrclf_basis_shell_center[chbrclf_shell_num] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5}; -chbrclf_basis_shell_ang_mom[chbrclf_shell_num] = +int32_t chbrclf_basis_shell_ang_mom[chbrclf_shell_num] = {0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 3, 3, 0, 0, 0, 0, 1, 1, 1, 2, 2, 0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 3, 3, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3}; -chbrclf_basis_shell_prim_index[chbrclf_shell_num] = +int64_t chbrclf_basis_shell_prim_num[chbrclf_shell_num] = + {10, 10, 1, 1, 1, 5, 1, 1, 1, 1, 1, 1, 1, 1, 5, 1, 1, 1, 1, 1, 1, 1, 1, 10, + 10, 1, 1, 1, 5, 1, 1, 1, 1, 1, 1, 1, 1, 15, 15, 15, 1, 1, 1, 9, 9, 1, 1, 1, + 1, 1, 1, 1, 1, 20, 20, 20, 20, 1, 1, 1, 13, 13, 13, 1, 1, 1, 8, 1, 1, 1, 1, + 1}; + +int64_t chbrclf_basis_shell_prim_index[chbrclf_shell_num] = {1, 11, 21, 22, 23, 24, 29, 30, 31, 32, 33, 34, 35, 36, 37, 42, 43, 44, 45, 46, 47, 48, 49, 50, 60, 70, 71, 72, 73, 78, 79, 80, 81, 82, 83, 84, 85, 86, 101, 116, 131, 132, 133, 134, 143, 152, 153, 154, 155, 156, 157, 158, 159, 160, 180, 200, 220, 240, 241, 242, 243, 256, 269, 282, 283, 284, 285, 293, 294, 295, 296, 297}; -chbrclf_basis_shell_prim_num[chbrclf_shell_num] = - {10, 10, 1, 1, 1, 5, 1, 1, 1, 1, 1, 1, 1, 1, 5, 1, 1, 1, 1, 1, 1, 1, 1, 10, - 10, 1, 1, 1, 5, 1, 1, 1, 1, 1, 1, 1, 1, 15, 15, 15, 1, 1, 1, 9, 9, 1, 1, 1, - 1, 1, 1, 1, 1, 20, 20, 20, 20, 1, 1, 1, 13, 13, 13, 1, 1, 1, 8, 1, 1, 1, 1, - 1}; - -chbrclf_basis_shell_factor[chbrclf_shell_num] = +double chbrclf_basis_shell_factor[chbrclf_shell_num] = {1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.}; -chbrclf_basis_exponent[chbrclf_prim_num] = +double chbrclf_basis_exponent[chbrclf_prim_num] = {8.2360000000000000e+03, 1.2350000000000000e+03, 2.8080000000000001e+02, 7.9269999999999996e+01, 2.5590000000000000e+01, 8.9969999999999999e+00, 3.3190000000000000e+00, 9.0590000000000004e-01, 3.6430000000000001e-01, @@ -657,7 +657,7 @@ chbrclf_basis_exponent[chbrclf_prim_num] = 7.0618999999999998e-01, 7.0618999999999998e-01, 2.6390000000000002e-01, 1.0470000000000000e-01, 5.5149999999999999e-01, 2.5800000000000001e-01}; -chbrclf_basis_coefficient[chbrclf_prim_num] = +double chbrclf_basis_coefficient[chbrclf_prim_num] = {5.3100000000000000e-04, 4.1079999999999997e-03, 2.1087000000000002e-02, 8.1852999999999995e-02, 2.3481700000000000e-01, 4.3440099999999998e-01, 3.4612900000000002e-01, 3.9378000000000003e-02, -8.9829999999999997e-03, @@ -758,7 +758,7 @@ chbrclf_basis_coefficient[chbrclf_prim_num] = 7.4396699999999996e-02, 1.0000000000000000e+00, 1.0000000000000000e+00, 1.0000000000000000e+00, 1.0000000000000000e+00, 1.0000000000000000e+00}; -chbrclf_basis_prim_factor[chbrclf_prim_num] = +double chbrclf_basis_prim_factor[chbrclf_prim_num] = {6.1616545431994848e+02, 1.4847738511079908e+02, 4.8888635917437597e+01, 1.8933972232608955e+01, 8.1089160941724145e+00, 3.7024003863155635e+00, 1.7525302846177560e+00, 6.6179013183966806e-01, 3.3419848027174592e-01, From 4ae55176410e91ba9f639240b2799f95fc8fd4ac Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 10 Jun 2021 00:10:19 +0200 Subject: [PATCH 6/9] Working on AOs --- org/qmckl_ao.org | 204 +++++++++++++++++++++++++++++++---------- org/qmckl_context.org | 2 +- org/qmckl_electron.org | 2 +- org/qmckl_error.org | 9 +- org/qmckl_tests.org | 17 ++-- 5 files changed, 171 insertions(+), 63 deletions(-) diff --git a/org/qmckl_ao.org b/org/qmckl_ao.org index 64cb745..f360291 100644 --- a/org/qmckl_ao.org +++ b/org/qmckl_ao.org @@ -89,17 +89,25 @@ int main() { The following arrays are stored in the context: - | ~type~ | | Gaussian (~'G'~) or Slater (~'S'~) | - | ~shell_num~ | | Number of shells | - | ~prim_num~ | | Total number of primitives | - | ~shell_center~ | ~[shell_num]~ | Id of the nucleus on which each shell is centered | - | ~shell_ang_mom~ | ~[shell_num]~ | Angular momentum of each shell | - | ~shell_prim_num~ | ~[shell_num]~ | Number of primitives in each shell | - | ~shell_prim_index~ | ~[shell_num]~ | Address of the first primitive of each shell in the ~EXPONENT~ array | - | ~shell_factor~ | ~[shell_num]~ | Normalization factor for each shell | - | ~exponent~ | ~[prim_num]~ | Array of exponents | - | ~coefficient~ | ~[prim_num]~ | Array of coefficients | - | ~prim_factor~ | ~[prim_num]~ | Normalization factors of the primtives | + |----------------------+---------------+----------------------------------------------------------------------| + | ~type~ | | Gaussian (~'G'~) or Slater (~'S'~) | + | ~shell_num~ | | Number of shells | + | ~prim_num~ | | Total number of primitives | + | ~nucleus_index~ | ~[nucl_num]~ | Index of the first shell of each nucleus | + | ~shell_ang_mom~ | ~[shell_num]~ | Angular momentum of each shell | + | ~shell_prim_num~ | ~[shell_num]~ | Number of primitives in each shell | + | ~shell_prim_index~ | ~[shell_num]~ | Address of the first primitive of each shell in the ~EXPONENT~ array | + | ~shell_factor~ | ~[shell_num]~ | Normalization factor for each shell | + | ~exponent~ | ~[prim_num]~ | Array of exponents | + | ~coefficient~ | ~[prim_num]~ | Array of coefficients | + | ~prim_factor~ | ~[prim_num]~ | Normalization factors of the primtives | + |----------------------+---------------+----------------------------------------------------------------------| + | ~nucl_shell_index~ | ~[nucl_num]~ | Index of the first shell for each nucleus | + | ~exponent_sorted~ | ~[prim_num]~ | Array of exponents for sorted primitives | + | ~coeff_norm_sorted~ | ~[prim_num]~ | Array of normalized coefficients for sorted primitives | + | ~prim_factor_sorted~ | ~[prim_num]~ | Normalization factors of the sorted primtives | + | ~nuclear_radius~ | ~[nucl_num]~ | Distance beyond which all the AOs are zero | + |----------------------+---------------+----------------------------------------------------------------------| For H_2 with the following basis set, @@ -129,11 +137,11 @@ D 1 type = 'G' shell_num = 12 prim_num = 20 -shell_center = [1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2] +nucleus_index = [0 , 6] shell_ang_mom = [0, 0, 0, 1, 1, 2, 0, 0, 0, 1, 1, 2] shell_factor = [ 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.] shell_prim_num = [5, 1, 1, 1, 1, 1, 5, 1, 1, 1, 1, 1] -shell_prim_index = [1, 6, 7, 8, 9, 10, 11, 16, 17, 18, 19, 20] +shell_prim_index = [0 , 5 , 6 , 7 , 8 , 9 , 10, 15, 16, 17, 18, 19] exponent = [ 33.87, 5.095, 1.159, 0.3258, 0.1027, 0.3258, 0.1027, 1.407, 0.388, 1.057, 33.87, 5.095, 1.159, 0.3258, 0.1027, 0.3258, 0.1027, 1.407, 0.388, 1.057] @@ -156,7 +164,8 @@ typedef struct qmckl_ao_basis_struct { int32_t uninitialized; int64_t shell_num; int64_t prim_num; - int64_t * shell_center; + int64_t * nucleus_index; + int64_t * nucleus_shell_num; int32_t * shell_ang_mom; int64_t * shell_prim_num; int64_t * shell_prim_index; @@ -180,7 +189,7 @@ typedef struct qmckl_ao_basis_struct { char qmckl_get_ao_basis_type (const qmckl_context context); int64_t qmckl_get_ao_basis_shell_num (const qmckl_context context); int64_t qmckl_get_ao_basis_prim_num (const qmckl_context context); -int64_t* qmckl_get_ao_basis_shell_center (const qmckl_context context); +int64_t* qmckl_get_ao_basis_nucleus_index (const qmckl_context context); int32_t* qmckl_get_ao_basis_shell_ang_mom (const qmckl_context context); int64_t* qmckl_get_ao_basis_shell_prim_num (const qmckl_context context); int64_t* qmckl_get_ao_basis_shell_prim_index (const qmckl_context context); @@ -263,7 +272,7 @@ int64_t qmckl_get_ao_basis_prim_num (const qmckl_context context) { } -int64_t* qmckl_get_ao_basis_shell_center (const qmckl_context context) { +int64_t* qmckl_get_ao_basis_nucleus_shell_num (const qmckl_context context) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return NULL; } @@ -277,8 +286,26 @@ int64_t* qmckl_get_ao_basis_shell_center (const qmckl_context context) { return NULL; } - assert (ctx->ao_basis.shell_center != NULL); - return ctx->ao_basis.shell_center; + assert (ctx->ao_basis.nucleus_shell_num != NULL); + return ctx->ao_basis.nucleus_shell_num ; +} + +int64_t* qmckl_get_ao_basis_nucleus_index (const qmckl_context context) { + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return NULL; + } + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + int32_t mask = 1 << 4; + + if ( (ctx->ao_basis.uninitialized & mask) != 0) { + return NULL; + } + + assert (ctx->ao_basis.nucleus_index != NULL); + return ctx->ao_basis.nucleus_index ; } @@ -290,7 +317,7 @@ int32_t* qmckl_get_ao_basis_shell_ang_mom (const qmckl_context context) { qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - int32_t mask = 1 << 4; + int32_t mask = 1 << 5; if ( (ctx->ao_basis.uninitialized & mask) != 0) { return NULL; @@ -309,7 +336,7 @@ int64_t* qmckl_get_ao_basis_shell_prim_num (const qmckl_context context) { qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - int32_t mask = 1 << 5; + int32_t mask = 1 << 6; if ( (ctx->ao_basis.uninitialized & mask) != 0) { return NULL; @@ -328,7 +355,7 @@ int64_t* qmckl_get_ao_basis_shell_prim_index (const qmckl_context context) { qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - int32_t mask = 1 << 6; + int32_t mask = 1 << 7; if ( (ctx->ao_basis.uninitialized & mask) != 0) { return NULL; @@ -347,7 +374,7 @@ double* qmckl_get_ao_basis_shell_factor (const qmckl_context context) { qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - int32_t mask = 1 << 7; + int32_t mask = 1 << 8; if ( (ctx->ao_basis.uninitialized & mask) != 0) { return NULL; @@ -367,7 +394,7 @@ double* qmckl_get_ao_basis_exponent (const qmckl_context context) { assert (ctx != NULL); - int32_t mask = 1 << 8; + int32_t mask = 1 << 9; if ( (ctx->ao_basis.uninitialized & mask) != 0) { return NULL; @@ -386,7 +413,7 @@ double* qmckl_get_ao_basis_coefficient (const qmckl_context context) { qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - int32_t mask = 1 << 9; + int32_t mask = 1 << 10; if ( (ctx->ao_basis.uninitialized & mask) != 0) { return NULL; @@ -405,7 +432,7 @@ double* qmckl_get_ao_basis_prim_factor (const qmckl_context context) { qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - int32_t mask = 1 << 10; + int32_t mask = 1 << 11; if ( (ctx->ao_basis.uninitialized & mask) != 0) { return NULL; @@ -432,16 +459,17 @@ bool qmckl_ao_basis_provided(const qmckl_context context) { ** Initialization functions To set the basis set, all the following functions need to be - called. When + called. #+begin_src c :comments org :tangle (eval h_func) qmckl_exit_code qmckl_set_ao_basis_type (qmckl_context context, const char t); qmckl_exit_code qmckl_set_ao_basis_shell_num (qmckl_context context, const int64_t shell_num); qmckl_exit_code qmckl_set_ao_basis_prim_num (qmckl_context context, const int64_t prim_num); -qmckl_exit_code qmckl_set_ao_basis_shell_prim_index (qmckl_context context, const int64_t * shell_prim_index); -qmckl_exit_code qmckl_set_ao_basis_shell_center (qmckl_context context, const int64_t * shell_center); +qmckl_exit_code qmckl_set_ao_basis_nucleus_index (qmckl_context context, const int64_t * nucleus_index); +qmckl_exit_code qmckl_set_ao_basis_nucleus_shell_num(qmckl_context context, const int64_t * nucleus_shell_num); qmckl_exit_code qmckl_set_ao_basis_shell_ang_mom (qmckl_context context, const int32_t * shell_ang_mom); qmckl_exit_code qmckl_set_ao_basis_shell_prim_num (qmckl_context context, const int64_t * shell_prim_num); +qmckl_exit_code qmckl_set_ao_basis_shell_prim_index (qmckl_context context, const int64_t * shell_prim_index); qmckl_exit_code qmckl_set_ao_basis_shell_factor (qmckl_context context, const double * shell_factor); qmckl_exit_code qmckl_set_ao_basis_exponent (qmckl_context context, const double * exponent); qmckl_exit_code qmckl_set_ao_basis_coefficient (qmckl_context context, const double * coefficient); @@ -535,7 +563,7 @@ qmckl_exit_code qmckl_set_ao_basis_prim_num(qmckl_context context, const int64_ } -qmckl_exit_code qmckl_set_ao_basis_shell_center(qmckl_context context, const int64_t* shell_center) { +qmckl_exit_code qmckl_set_ao_basis_nucleus_shell_num(qmckl_context context, const int64_t* nucleus_shell_num) { <> int32_t mask = 1 << 3; @@ -544,15 +572,15 @@ qmckl_exit_code qmckl_set_ao_basis_shell_center(qmckl_context context, const in if (shell_num == 0L) { return qmckl_failwith( context, QMCKL_FAILURE, - "qmckl_set_ao_basis_shell_center", + "qmckl_set_ao_basis_nucleus_shell_num", "shell_num is not set"); } - if (ctx->ao_basis.shell_center != NULL) { - qmckl_exit_code rc = qmckl_free(context, ctx->ao_basis.shell_center); + if (ctx->ao_basis.nucleus_shell_num != NULL) { + qmckl_exit_code rc = qmckl_free(context, ctx->ao_basis.nucleus_shell_num); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, - "qmckl_set_ao_basis_shell_center", + "qmckl_set_ao_basis_nucleus_shell_num", NULL); } } @@ -564,13 +592,53 @@ qmckl_exit_code qmckl_set_ao_basis_shell_center(qmckl_context context, const in if (new_array == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, - "qmckl_set_ao_basis_shell_center", + "qmckl_set_ao_basis_nucleus_shell_num", NULL); } - memcpy(new_array, shell_center, mem_info.size); + memcpy(new_array, nucleus_shell_num, mem_info.size); - ctx->ao_basis.shell_center = new_array; + ctx->ao_basis.nucleus_shell_num = new_array; + + <> +} + +qmckl_exit_code qmckl_set_ao_basis_nucleus_index(qmckl_context context, const int64_t* nucleus_index) { + <> + + int32_t mask = 1 << 4; + + const int64_t shell_num = qmckl_get_ao_basis_shell_num(context); + if (shell_num == 0L) { + return qmckl_failwith( context, + QMCKL_FAILURE, + "qmckl_set_ao_basis_nucleus_index", + "shell_num is not set"); + } + + if (ctx->ao_basis.nucleus_index != NULL) { + qmckl_exit_code rc = qmckl_free(context, ctx->ao_basis.nucleus_index); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_set_ao_basis_nucleus_index", + NULL); + } + } + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = shell_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_ao_basis_nucleus_index", + NULL); + } + + memcpy(new_array, nucleus_index, mem_info.size); + + ctx->ao_basis.nucleus_index = new_array; <> } @@ -579,7 +647,7 @@ qmckl_exit_code qmckl_set_ao_basis_shell_center(qmckl_context context, const in qmckl_exit_code qmckl_set_ao_basis_shell_ang_mom(qmckl_context context, const int32_t* shell_ang_mom) { <> - int32_t mask = 1 << 4; + int32_t mask = 1 << 5; const int64_t shell_num = qmckl_get_ao_basis_shell_num(context); if (shell_num == 0L) { @@ -621,7 +689,7 @@ qmckl_exit_code qmckl_set_ao_basis_shell_ang_mom(qmckl_context context, const i qmckl_exit_code qmckl_set_ao_basis_shell_prim_num(qmckl_context context, const int64_t* shell_prim_num) { <> - int32_t mask = 1 << 5; + int32_t mask = 1 << 6; const int64_t shell_num = qmckl_get_ao_basis_shell_num(context); if (shell_num == 0L) { @@ -663,7 +731,7 @@ qmckl_exit_code qmckl_set_ao_basis_shell_prim_num(qmckl_context context, const qmckl_exit_code qmckl_set_ao_basis_shell_prim_index(qmckl_context context, const int64_t* shell_prim_index) { <> - int32_t mask = 1 << 6; + int32_t mask = 1 << 7; const int64_t shell_num = qmckl_get_ao_basis_shell_num(context); if (shell_num == 0L) { @@ -704,7 +772,7 @@ qmckl_exit_code qmckl_set_ao_basis_shell_prim_index(qmckl_context context, cons qmckl_exit_code qmckl_set_ao_basis_shell_factor(qmckl_context context, const double* shell_factor) { <> - int32_t mask = 1 << 7; + int32_t mask = 1 << 8; const int64_t shell_num = qmckl_get_ao_basis_shell_num(context); if (shell_num == 0L) { @@ -745,7 +813,7 @@ qmckl_exit_code qmckl_set_ao_basis_shell_factor(qmckl_context context, const do qmckl_exit_code qmckl_set_ao_basis_exponent(qmckl_context context, const double* exponent) { <> - int32_t mask = 1 << 8; + int32_t mask = 1 << 9; const int64_t prim_num = qmckl_get_ao_basis_prim_num(context); if (prim_num == 0L) { @@ -785,7 +853,7 @@ qmckl_exit_code qmckl_set_ao_basis_exponent(qmckl_context context, const double qmckl_exit_code qmckl_set_ao_basis_coefficient(qmckl_context context, const double* coefficient) { <> - int32_t mask = 1 << 9; + int32_t mask = 1 << 10; const int64_t prim_num = qmckl_get_ao_basis_prim_num(context); if (prim_num == 0L) { @@ -825,7 +893,7 @@ qmckl_exit_code qmckl_set_ao_basis_coefficient(qmckl_context context, const dou qmckl_exit_code qmckl_set_ao_basis_prim_factor(qmckl_context context, const double* prim_factor) { <> - int32_t mask = 1 << 10; + int32_t mask = 1 << 11; const int64_t prim_num = qmckl_get_ao_basis_prim_num(context); if (prim_num == 0L) { @@ -865,6 +933,41 @@ qmckl_exit_code qmckl_set_ao_basis_prim_factor(qmckl_context context, const dou #+end_src + When the basis set is completely entered, other data structures are + computed to accelerate the calculations. The primitives within each + contraction are sorted in ascending order of their exponents, such + that as soon as a primitive is zero all the following functions + vanish. Also, it is possible to compute a nuclear radius beyond which + all the primitives are zero up to the numerical accuracy defined in + the context. + + # TODO : sort the basis set here + + #+begin_src c :comments org :tangle (eval h_private_type) :noweb yes :exports none +qmckl_exit_code qmckl_finalize_basis(qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_finalize_basis(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); + + int64_t nucl_num = 0; + qmckl_exit_code rc = QMCKL_FAILURE; + + rc = qmckl_get_nucleus_num(context, &nucl_num); + if (rc != QMCKL_SUCCESS) return rc; + + /* TODO : sort the basis set here */ + return QMCKL_SUCCESS; +} + #+end_src + ** TODO Fortran interfaces ** Test :noexport: @@ -872,7 +975,8 @@ qmckl_exit_code qmckl_set_ao_basis_prim_factor(qmckl_context context, const dou #+begin_src c :tangle (eval c_test) :exports none :exports none const int64_t shell_num = chbrclf_shell_num; const int64_t prim_num = chbrclf_prim_num; -const int64_t * shell_center = &(chbrclf_basis_shell_center[0]); +const int64_t * nucleus_index = &(chbrclf_basis_nucleus_index[0]); +const int64_t * nucleus_shell_num = &(chbrclf_basis_nucleus_shell_num[0]); const int32_t * shell_ang_mom = &(chbrclf_basis_shell_ang_mom[0]); const int64_t * shell_prim_num = &(chbrclf_basis_shell_prim_num[0]); const int64_t * shell_prim_index = &(chbrclf_basis_shell_prim_index[0]); @@ -899,7 +1003,11 @@ rc = qmckl_set_ao_basis_prim_num (context, prim_num); assert(rc == QMCKL_SUCCESS); assert(!qmckl_ao_basis_provided(context)); -rc = qmckl_set_ao_basis_shell_center (context, shell_center); +rc = qmckl_set_ao_basis_nucleus_index (context, nucleus_index); +assert(rc == QMCKL_SUCCESS); +assert(!qmckl_ao_basis_provided(context)); + +rc = qmckl_set_ao_basis_nucleus_shell_num (context, nucleus_shell_num); assert(rc == QMCKL_SUCCESS); assert(!qmckl_ao_basis_provided(context)); @@ -911,10 +1019,6 @@ rc = qmckl_set_ao_basis_shell_factor (context, shell_factor); assert(rc == QMCKL_SUCCESS); assert(!qmckl_ao_basis_provided(context)); -rc = qmckl_set_ao_basis_shell_center (context, shell_prim_num); -assert(rc == QMCKL_SUCCESS); -assert(!qmckl_ao_basis_provided(context)); - rc = qmckl_set_ao_basis_shell_prim_num (context, shell_prim_num); assert(rc == QMCKL_SUCCESS); assert(!qmckl_ao_basis_provided(context)); @@ -1523,6 +1627,7 @@ end function test_qmckl_ao_polynomial_vgl \[ \nabla_z v_i = -2 a_i (X_z - R_z) v_i \] \[ \Delta v_i = a_i (4 |X-R|^2 a_i - 6) v_i \] + |--------------+--------+------------------------------------------------------| | ~context~ | input | Global state | | ~X(3)~ | input | Array containing the coordinates of the points | | ~R(3)~ | input | Array containing the x,y,z coordinates of the center | @@ -1530,6 +1635,7 @@ end function test_qmckl_ao_polynomial_vgl | ~A(n)~ | input | Exponents of the Gaussians | | ~VGL(ldv,5)~ | output | Value, gradients and Laplacian of the Gaussians | | ~ldv~ | input | Leading dimension of array ~VGL~ | + |--------------+--------+------------------------------------------------------| Requirements diff --git a/org/qmckl_context.org b/org/qmckl_context.org index 618120d..6c7a3ac 100644 --- a/org/qmckl_context.org +++ b/org/qmckl_context.org @@ -221,7 +221,7 @@ qmckl_context qmckl_context_create() { ctx->numprec.precision = QMCKL_DEFAULT_PRECISION; ctx->numprec.range = QMCKL_DEFAULT_RANGE; - ctx->ao_basis.uninitialized = (1 << 11) - 1; + ctx->ao_basis.uninitialized = (1 << 12) - 1; ctx->nucleus.uninitialized = (1 << 4) - 1; ctx->electron.uninitialized = (1 << 3) - 1; diff --git a/org/qmckl_electron.org b/org/qmckl_electron.org index bc56cda..11f025f 100644 --- a/org/qmckl_electron.org +++ b/org/qmckl_electron.org @@ -1221,7 +1221,7 @@ assert(qmckl_electron_provided(context)); double ee_distance_rescaled[walk_num * elec_num * elec_num]; -rc = qmckl_get_electron_ee_distance_rescaled(context, ee_distance); +rc = qmckl_get_electron_ee_distance_rescaled(context, ee_distance_rescaled); // TODO: Get exact values //// (e1,e2,w) diff --git a/org/qmckl_error.org b/org/qmckl_error.org index 0d50308..0849de2 100644 --- a/org/qmckl_error.org +++ b/org/qmckl_error.org @@ -104,7 +104,8 @@ typedef int32_t qmckl_exit_code; | ~QMCKL_ALLOCATION_FAILED~ | 104 | 'Allocation failed' | | ~QMCKL_DEALLOCATION_FAILED~ | 105 | 'De-allocation failed' | | ~QMCKL_NOT_PROVIDED~ | 106 | 'Not provided' | - | ~QMCKL_INVALID_EXIT_CODE~ | 107 | 'Invalid exit code' | + | ~QMCKL_OUT_OF_BOUNDS~ | 107 | 'Index out of bounds' | + | ~QMCKL_INVALID_EXIT_CODE~ | 108 | 'Invalid exit code' | # We need to force Emacs not to indent the Python code: # -*- org-src-preserve-indentation: t @@ -162,7 +163,8 @@ return '\n'.join(result) #define QMCKL_ALLOCATION_FAILED ((qmckl_exit_code) 104) #define QMCKL_DEALLOCATION_FAILED ((qmckl_exit_code) 105) #define QMCKL_NOT_PROVIDED ((qmckl_exit_code) 106) - #define QMCKL_INVALID_EXIT_CODE ((qmckl_exit_code) 107) + #define QMCKL_OUT_OF_BOUNDS ((qmckl_exit_code) 107) + #define QMCKL_INVALID_EXIT_CODE ((qmckl_exit_code) 108) #+end_src #+begin_src f90 :comments org :tangle (eval fh_type) :exports none @@ -193,7 +195,8 @@ return '\n'.join(result) integer(qmckl_exit_code), parameter :: QMCKL_ALLOCATION_FAILED = 104 integer(qmckl_exit_code), parameter :: QMCKL_DEALLOCATION_FAILED = 105 integer(qmckl_exit_code), parameter :: QMCKL_NOT_PROVIDED = 106 - integer(qmckl_exit_code), parameter :: QMCKL_INVALID_EXIT_CODE = 107 + integer(qmckl_exit_code), parameter :: QMCKL_OUT_OF_BOUNDS = 107 + integer(qmckl_exit_code), parameter :: QMCKL_INVALID_EXIT_CODE = 108 #+end_src :end: diff --git a/org/qmckl_tests.org b/org/qmckl_tests.org index 3c090d1..a2fbbfe 100644 --- a/org/qmckl_tests.org +++ b/org/qmckl_tests.org @@ -527,10 +527,9 @@ F 1 #define chbrclf_shell_num 72 #define chbrclf_prim_num 297 -int64_t chbrclf_basis_shell_center[chbrclf_shell_num] = - {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, - 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5}; +int64_t chbrclf_basis_nucleus_index[chbrclf_nucl_num] = {0, 14, 23, 27, 53}; + +int64_t chbrclf_basis_nucleus_shell_num[chbrclf_nucl_num] = {14, 9, 14, 16, 19}; int32_t chbrclf_basis_shell_ang_mom[chbrclf_shell_num] = {0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 3, 3, 0, 0, 0, 0, 1, 1, 1, 2, 2, 0, @@ -544,11 +543,11 @@ int64_t chbrclf_basis_shell_prim_num[chbrclf_shell_num] = 1}; int64_t chbrclf_basis_shell_prim_index[chbrclf_shell_num] = - {1, 11, 21, 22, 23, 24, 29, 30, 31, 32, 33, 34, 35, 36, 37, 42, 43, 44, 45, - 46, 47, 48, 49, 50, 60, 70, 71, 72, 73, 78, 79, 80, 81, 82, 83, 84, 85, 86, - 101, 116, 131, 132, 133, 134, 143, 152, 153, 154, 155, 156, 157, 158, 159, - 160, 180, 200, 220, 240, 241, 242, 243, 256, 269, 282, 283, 284, 285, 293, - 294, 295, 296, 297}; + {0, 10, 20, 21, 22, 23, 28, 29, 30, 31, 32, 33, 34, 35, 36, 41, 42, 43, 44, + 45, 46, 47, 48, 49, 59, 69, 70, 71, 72, 77, 78, 79, 80, 81, 82, 83, 84, 85, + 100, 115, 130, 131, 132, 133, 142, 151, 152, 153, 154, 155, 156, 157, 158, + 159, 179, 199, 219, 239, 240, 241, 242, 255, 268, 281, 282, 283, 284, 292, + 293, 294, 295, 296}; double chbrclf_basis_shell_factor[chbrclf_shell_num] = {1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., From 4bcb9b980cd0523f2fa3c3f32903eb5859076c99 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 10 Jun 2021 22:57:59 +0200 Subject: [PATCH 7/9] Implemented AO primitives --- org/qmckl_ao.org | 902 ++++++++++++++++++++++++++++++----------- org/qmckl_context.org | 32 +- org/qmckl_electron.org | 199 +++++---- org/qmckl_nucleus.org | 95 +++-- 4 files changed, 855 insertions(+), 373 deletions(-) diff --git a/org/qmckl_ao.org b/org/qmckl_ao.org index f360291..abdb227 100644 --- a/org/qmckl_ao.org +++ b/org/qmckl_ao.org @@ -56,6 +56,7 @@ gradients and Laplacian of the atomic basis functions. #include "config.h" #endif +#include #include "chbrclf.h" int main() { @@ -83,31 +84,38 @@ int main() { #include "qmckl_context_private_type.h" #include "qmckl_memory_private_type.h" #include "qmckl_memory_private_func.h" +#include "qmckl_ao_private_type.h" +#include "qmckl_ao_private_func.h" #+end_src * Context The following arrays are stored in the context: - |----------------------+---------------+----------------------------------------------------------------------| - | ~type~ | | Gaussian (~'G'~) or Slater (~'S'~) | - | ~shell_num~ | | Number of shells | - | ~prim_num~ | | Total number of primitives | - | ~nucleus_index~ | ~[nucl_num]~ | Index of the first shell of each nucleus | - | ~shell_ang_mom~ | ~[shell_num]~ | Angular momentum of each shell | - | ~shell_prim_num~ | ~[shell_num]~ | Number of primitives in each shell | - | ~shell_prim_index~ | ~[shell_num]~ | Address of the first primitive of each shell in the ~EXPONENT~ array | - | ~shell_factor~ | ~[shell_num]~ | Normalization factor for each shell | - | ~exponent~ | ~[prim_num]~ | Array of exponents | - | ~coefficient~ | ~[prim_num]~ | Array of coefficients | - | ~prim_factor~ | ~[prim_num]~ | Normalization factors of the primtives | - |----------------------+---------------+----------------------------------------------------------------------| - | ~nucl_shell_index~ | ~[nucl_num]~ | Index of the first shell for each nucleus | - | ~exponent_sorted~ | ~[prim_num]~ | Array of exponents for sorted primitives | - | ~coeff_norm_sorted~ | ~[prim_num]~ | Array of normalized coefficients for sorted primitives | - | ~prim_factor_sorted~ | ~[prim_num]~ | Normalization factors of the sorted primtives | - | ~nuclear_radius~ | ~[nucl_num]~ | Distance beyond which all the AOs are zero | - |----------------------+---------------+----------------------------------------------------------------------| + |--------------------+---------------+----------------------------------------------------------------------| + | ~type~ | | Gaussian (~'G'~) or Slater (~'S'~) | + | ~shell_num~ | | Number of shells | + | ~prim_num~ | | Total number of primitives | + | ~nucleus_index~ | ~[nucl_num]~ | Index of the first shell of each nucleus | + | ~shell_ang_mom~ | ~[shell_num]~ | Angular momentum of each shell | + | ~shell_prim_num~ | ~[shell_num]~ | Number of primitives in each shell | + | ~shell_prim_index~ | ~[shell_num]~ | Address of the first primitive of each shell in the ~EXPONENT~ array | + | ~shell_factor~ | ~[shell_num]~ | Normalization factor for each shell | + | ~exponent~ | ~[prim_num]~ | Array of exponents | + | ~coefficient~ | ~[prim_num]~ | Array of coefficients | + | ~prim_factor~ | ~[prim_num]~ | Normalization factors of the primtives | + + Computed data + |----------------------+-------------------------------------+-----------------------------------------------------------------------------------------------| + | ~nucleus_prim_index~ | ~[nucl_num]~ | Index of the first primitive for each nucleus | + | ~primitive_vgl~ | ~[prim_num][5][walk_num][elec_num]~ | Value, gradients, Laplacian of the primitives at electron positions | + | ~primitive_vgl_date~ | ~uint64_t~ | Late modification date of Value, gradients, Laplacian of the primitives at electron positions | + |----------------------+-------------------------------------+-----------------------------------------------------------------------------------------------| + | ~nucl_shell_index~ | ~[nucl_num]~ | Index of the first shell for each nucleus | + | ~exponent_sorted~ | ~[prim_num]~ | Array of exponents for sorted primitives | + | ~coeff_norm_sorted~ | ~[prim_num]~ | Array of normalized coefficients for sorted primitives | + | ~prim_factor_sorted~ | ~[prim_num]~ | Normalization factors of the sorted primtives | + | ~nuclear_radius~ | ~[nucl_num]~ | Distance beyond which all the AOs are zero | For H_2 with the following basis set, @@ -168,21 +176,51 @@ typedef struct qmckl_ao_basis_struct { int64_t * nucleus_shell_num; int32_t * shell_ang_mom; int64_t * shell_prim_num; + int64_t * nucleus_prim_index; int64_t * shell_prim_index; double * shell_factor; double * exponent ; double * coefficient ; double * prim_factor ; + double * primitive_vgl; + int64_t primitive_vgl_date; bool provided; char type; } qmckl_ao_basis_struct; #+end_src The ~uninitialized~ integer contains one bit set to one for each - initialization function which has not bee called. It becomes equal + 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_private_func) +qmckl_exit_code qmckl_init_ao_basis(qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) +qmckl_exit_code qmckl_init_ao_basis(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->ao_basis.uninitialized = (1 << 12) - 1; + + /* Default values */ + /* ctx->ao_basis. + ,*/ + + return QMCKL_SUCCESS; +} + #+end_src + + ** Access functions #+begin_src c :comments org :tangle (eval h_private_func) :exports none @@ -489,6 +527,10 @@ qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; #+begin_src c :exports none ctx->ao_basis.uninitialized &= ~mask; ctx->ao_basis.provided = (ctx->ao_basis.uninitialized == 0); +if (ctx->ao_basis.provided) { + qmckl_exit_code rc_ = qmckl_finalize_basis(context); + if (rc_ != QMCKL_SUCCESS) return rc_; + } return QMCKL_SUCCESS; #+end_src @@ -941,9 +983,7 @@ qmckl_exit_code qmckl_set_ao_basis_prim_factor(qmckl_context context, const dou all the primitives are zero up to the numerical accuracy defined in the context. - # TODO : sort the basis set here - - #+begin_src c :comments org :tangle (eval h_private_type) :noweb yes :exports none + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_finalize_basis(qmckl_context context); #+end_src @@ -962,7 +1002,27 @@ qmckl_exit_code qmckl_finalize_basis(qmckl_context context) { rc = qmckl_get_nucleus_num(context, &nucl_num); if (rc != QMCKL_SUCCESS) return rc; + + /* nucleus_prim_index */ + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = (ctx->nucleus.num + (int64_t) 1) * sizeof(int64_t); + + ctx->ao_basis.nucleus_prim_index = (int64_t*) qmckl_malloc(context, mem_info); + + if (ctx->ao_basis.nucleus_prim_index == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_nucleus_prim_index", + NULL); + } + + for (int64_t i=0 ; iao_basis.nucleus_index[i]; + ctx->ao_basis.nucleus_prim_index[i] = ctx->ao_basis.shell_prim_index[shell_idx]; + } + ctx->ao_basis.nucleus_prim_index[nucl_num] = ctx->ao_basis.prim_num; + /* TODO : sort the basis set here */ return QMCKL_SUCCESS; } @@ -973,6 +1033,23 @@ qmckl_exit_code qmckl_finalize_basis(qmckl_context context) { ** Test :noexport: #+begin_src c :tangle (eval c_test) :exports none :exports none +const int64_t nucl_num = chbrclf_nucl_num; +const double* nucl_charge = chbrclf_charge; +const double* nucl_coord = &(chbrclf_nucl_coord[0][0]); + +qmckl_exit_code rc; +rc = qmckl_set_nucleus_num (context, nucl_num); +assert(rc == QMCKL_SUCCESS); + +rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0])); +assert(rc == QMCKL_SUCCESS); + +rc = qmckl_set_nucleus_charge(context, nucl_charge); +assert(rc == QMCKL_SUCCESS); + +assert(qmckl_nucleus_provided(context)); + + const int64_t shell_num = chbrclf_shell_num; const int64_t prim_num = chbrclf_prim_num; const int64_t * nucleus_index = &(chbrclf_basis_nucleus_index[0]); @@ -987,8 +1064,6 @@ const double * prim_factor = &(chbrclf_basis_prim_factor[0]); char typ = 'G'; -qmckl_exit_code rc; - assert(!qmckl_ao_basis_provided(context)); rc = qmckl_set_ao_basis_type (context, typ); @@ -1040,9 +1115,568 @@ assert(rc == QMCKL_SUCCESS); assert(qmckl_ao_basis_provided(context)); #+end_src + +* Radial part +** General functions for Gaussian basis functions + ~qmckl_ao_gaussian_vgl~ computes the values, gradients and + Laplacians at a given point of ~n~ Gaussian functions centered at + the same point: + + \[ v_i = \exp(-a_i |X-R|^2) \] + \[ \nabla_x v_i = -2 a_i (X_x - R_x) v_i \] + \[ \nabla_y v_i = -2 a_i (X_y - R_y) v_i \] + \[ \nabla_z v_i = -2 a_i (X_z - R_z) v_i \] + \[ \Delta v_i = a_i (4 |X-R|^2 a_i - 6) v_i \] + + |--------------+--------+------------------------------------------------------| + | ~context~ | input | Global state | + | ~X(3)~ | input | Array containing the coordinates of the points | + | ~R(3)~ | input | Array containing the x,y,z coordinates of the center | + | ~n~ | input | Number of computed Gaussians | + | ~A(n)~ | input | Exponents of the Gaussians | + | ~VGL(ldv,5)~ | output | Value, gradients and Laplacian of the Gaussians | + | ~ldv~ | input | Leading dimension of array ~VGL~ | + |--------------+--------+------------------------------------------------------| + + Requirements + + - ~context~ is not 0 + - ~n~ > 0 + - ~ldv~ >= 5 + - ~A(i)~ > 0 for all ~i~ + - ~X~ is allocated with at least $3 \times 8$ bytes + - ~R~ is allocated with at least $3 \times 8$ bytes + - ~A~ is allocated with at least $n \times 8$ bytes + - ~VGL~ is allocated with at least $n \times 5 \times 8$ bytes + + #+begin_src c :tangle (eval h_func) +qmckl_exit_code +qmckl_ao_gaussian_vgl(const qmckl_context context, + const double *X, + const double *R, + const int64_t *n, + const int64_t *A, + const double *VGL, + const int64_t ldv); + #+end_src + + #+begin_src f90 :tangle (eval f) +integer function qmckl_ao_gaussian_vgl_f(context, X, R, n, A, VGL, ldv) result(info) + use qmckl + implicit none + integer*8 , intent(in) :: context + real*8 , intent(in) :: X(3), R(3) + integer*8 , intent(in) :: n + real*8 , intent(in) :: A(n) + real*8 , intent(out) :: VGL(ldv,5) + integer*8 , intent(in) :: ldv + + integer*8 :: i,j + real*8 :: Y(3), r2, t, u, v + + info = QMCKL_SUCCESS + + if (context == QMCKL_NULL_CONTEXT) then + info = QMCKL_INVALID_CONTEXT + return + endif + + if (n <= 0) then + info = QMCKL_INVALID_ARG_4 + return + endif + + if (ldv < n) then + info = QMCKL_INVALID_ARG_7 + return + endif + + + do i=1,3 + Y(i) = X(i) - R(i) + end do + r2 = Y(1)*Y(1) + Y(2)*Y(2) + Y(3)*Y(3) + + do i=1,n + VGL(i,1) = dexp(-A(i) * r2) + end do + + do i=1,n + VGL(i,5) = A(i) * VGL(i,1) + end do + + t = -2.d0 * ( X(1) - R(1) ) + u = -2.d0 * ( X(2) - R(2) ) + v = -2.d0 * ( X(3) - R(3) ) + + do i=1,n + VGL(i,2) = t * VGL(i,5) + VGL(i,3) = u * VGL(i,5) + VGL(i,4) = v * VGL(i,5) + end do + + t = 4.d0 * r2 + do i=1,n + VGL(i,5) = (t * A(i) - 6.d0) * VGL(i,5) + end do + +end function qmckl_ao_gaussian_vgl_f + #+end_src + + #+begin_src f90 :tangle (eval f) :exports none +integer(c_int32_t) function qmckl_ao_gaussian_vgl(context, X, R, n, A, VGL, ldv) & + bind(C) result(info) + use, intrinsic :: iso_c_binding + implicit none + integer (c_int64_t) , intent(in) , value :: context + real (c_double) , intent(in) :: X(3), R(3) + integer (c_int64_t) , intent(in) , value :: n + real (c_double) , intent(in) :: A(n) + real (c_double) , intent(out) :: VGL(ldv,5) + integer (c_int64_t) , intent(in) , value :: ldv + + integer, external :: qmckl_ao_gaussian_vgl_f + info = qmckl_ao_gaussian_vgl_f(context, X, R, n, A, VGL, ldv) +end function qmckl_ao_gaussian_vgl + #+end_src + + #+begin_src f90 :tangle (eval fh_func) :exports none + interface + integer(c_int32_t) function qmckl_ao_gaussian_vgl(context, X, R, n, A, VGL, ldv) & + bind(C) + use, intrinsic :: iso_c_binding + integer (c_int64_t) , intent(in) , value :: context + integer (c_int64_t) , intent(in) , value :: ldv + integer (c_int64_t) , intent(in) , value :: n + real (c_double) , intent(in) :: X(3), R(3), A(n) + real (c_double) , intent(out) :: VGL(ldv,5) + end function qmckl_ao_gaussian_vgl + end interface + #+end_src + + # Test + #+begin_src f90 :tangle (eval f_test) +integer(c_int32_t) function test_qmckl_ao_gaussian_vgl(context) bind(C) + use qmckl + implicit none + + integer(c_int64_t), intent(in), value :: context + + integer*8 :: n, ldv, j, i + double precision :: X(3), R(3), Y(3), r2 + double precision, allocatable :: VGL(:,:), A(:) + double precision :: epsilon + + epsilon = qmckl_get_numprec_epsilon(context) + + X = (/ 1.1 , 2.2 , 3.3 /) + R = (/ 0.1 , 1.2 , -2.3 /) + Y(:) = X(:) - R(:) + r2 = Y(1)**2 + Y(2)**2 + Y(3)**2 + + n = 10; + ldv = 100; + + allocate (A(n), VGL(ldv,5)) + do i=1,n + A(i) = 0.0013 * dble(ishft(1,i)) + end do + + + test_qmckl_ao_gaussian_vgl = & + qmckl_ao_gaussian_vgl(context, X, R, n, A, VGL, ldv) + if (test_qmckl_ao_gaussian_vgl /= 0) return + + test_qmckl_ao_gaussian_vgl = -1 + + do i=1,n + test_qmckl_ao_gaussian_vgl = -11 + if (dabs(1.d0 - VGL(i,1) / (& + dexp(-A(i) * r2) & + )) > epsilon ) return + + test_qmckl_ao_gaussian_vgl = -12 + if (dabs(1.d0 - VGL(i,2) / (& + -2.d0 * A(i) * Y(1) * dexp(-A(i) * r2) & + )) > epsilon ) return + + test_qmckl_ao_gaussian_vgl = -13 + if (dabs(1.d0 - VGL(i,3) / (& + -2.d0 * A(i) * Y(2) * dexp(-A(i) * r2) & + )) > epsilon ) return + + test_qmckl_ao_gaussian_vgl = -14 + if (dabs(1.d0 - VGL(i,4) / (& + -2.d0 * A(i) * Y(3) * dexp(-A(i) * r2) & + )) > epsilon ) return + + test_qmckl_ao_gaussian_vgl = -15 + if (dabs(1.d0 - VGL(i,5) / (& + A(i) * (4.d0*r2*A(i) - 6.d0) * dexp(-A(i) * r2) & + )) > epsilon ) return + end do + + test_qmckl_ao_gaussian_vgl = 0 + + deallocate(VGL) +end function test_qmckl_ao_gaussian_vgl + #+end_src + + #+begin_src c :tangle (eval c_test) :exports none + int test_qmckl_ao_gaussian_vgl(qmckl_context context); + assert(0 == test_qmckl_ao_gaussian_vgl(context)); + #+end_src + +** TODO General functions for Slater basis functions +** TODO General functions for Radial functions on a grid +** Computation of primitives +*** Get + + #+begin_src c :comments org :tangle (eval h_func) :noweb yes +qmckl_exit_code qmckl_get_ao_basis_primitive_vgl(qmckl_context context, double* const primitive_vgl); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_get_ao_basis_primitive_vgl(qmckl_context context, double* const primitive_vgl) { + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + + qmckl_exit_code rc; + + rc = qmckl_provide_ao_basis_primitive_vgl(context); + if (rc != QMCKL_SUCCESS) return rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + size_t sze = ctx->ao_basis.prim_num * 5 * ctx->electron.num * ctx->electron.walk_num; + memcpy(primitive_vgl, ctx->ao_basis.primitive_vgl, sze * sizeof(double)); + + return QMCKL_SUCCESS; +} + #+end_src + +*** Provide + + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none +qmckl_exit_code qmckl_provide_ao_basis_primitive_vgl(qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_provide_ao_basis_primitive_vgl(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); + + if (!ctx->ao_basis.provided) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_ao_basis_primitive_vgl", + NULL); + } + + /* Compute if necessary */ + if (ctx->electron.coord_new_date > ctx->ao_basis.primitive_vgl_date) { + + /* Allocate array */ + if (ctx->ao_basis.primitive_vgl == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->ao_basis.prim_num * 5 * ctx->electron.num * + ctx->electron.walk_num * sizeof(double); + double* primitive_vgl = (double*) qmckl_malloc(context, mem_info); + + if (primitive_vgl == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_ao_basis_primitive_vgl", + NULL); + } + ctx->ao_basis.primitive_vgl = primitive_vgl; + } + + qmckl_exit_code rc; + if (ctx->ao_basis.type == 'G') { + rc = qmckl_compute_ao_basis_primitive_gaussian_vgl(context, + ctx->ao_basis.prim_num, + ctx->electron.num, + ctx->nucleus.num, + ctx->electron.walk_num, + ctx->ao_basis.nucleus_prim_index, + ctx->electron.coord_new, + ctx->nucleus.coord, + ctx->ao_basis.exponent, + ctx->ao_basis.primitive_vgl); + } else { + return qmckl_failwith( context, + QMCKL_FAILURE, + "compute_ao_basis_primitive_vgl", + "Not yet implemented"); + } + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->ao_basis.primitive_vgl_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + #+end_src + +*** Compute + :PROPERTIES: + :Name: qmckl_compute_ao_basis_primitive_gaussian_vgl + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_ao_basis_primitive_gaussian_vgl_args + | qmckl_context | context | in | Global state | + | int64_t | prim_num | in | Number of primitives | + | int64_t | elec_num | in | Number of electrons | + | int64_t | nucl_num | in | Number of nuclei | + | int64_t | walk_num | in | Number of walkers | + | int64_t | nucleus_prim_index[nucl_num] | in | Index of the 1st primitive of each nucleus | + | double | elec_coord[walk_num][3][elec_num] | in | Electron coordinates | + | double | nucl_coord[3][elec_num] | in | Nuclear coordinates | + | double | expo[prim_num] | in | Exponents of the primitives | + | double | primitive_vgl[prim_num][5][walk_num][elec_num] | out | Value, gradients and Laplacian of the primitives | + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer function qmckl_compute_ao_basis_primitive_gaussian_vgl_f(context, & + prim_num, elec_num, nucl_num, walk_num, & + nucleus_prim_index, elec_coord, nucl_coord, expo, primitive_vgl) & + result(info) + use qmckl + implicit none + integer(qmckl_context), intent(in) :: context + integer*8 , intent(in) :: prim_num + integer*8 , intent(in) :: nucl_num + integer*8 , intent(in) :: elec_num + integer*8 , intent(in) :: walk_num + integer*8 , intent(in) :: nucleus_prim_index(nucl_num+1) + double precision , intent(in) :: elec_coord(elec_num,3,walk_num) + double precision , intent(in) :: nucl_coord(nucl_num,3) + double precision , intent(in) :: expo(prim_num) + double precision , intent(out) :: primitive_vgl(elec_num,walk_num,5,prim_num) + + integer*8 :: inucl, iprim, iwalk, ielec + double precision :: x, y, z, two_a, ar2, r2, v + double precision :: r2_cut(elec_num,walk_num) + + info = QMCKL_SUCCESS + + if (context == QMCKL_NULL_CONTEXT) then + info = QMCKL_INVALID_CONTEXT + return + endif + + if (prim_num <= 0) then + info = QMCKL_INVALID_ARG_2 + return + endif + + if (elec_num <= 0) then + info = QMCKL_INVALID_ARG_4 + return + endif + + if (nucl_num <= 0) then + info = QMCKL_INVALID_ARG_5 + return + endif + + if (walk_num <= 0) then + info = QMCKL_INVALID_ARG_6 + return + endif + + do inucl=1,nucl_num + ! C is zero-based, so shift bounds by one + do iprim = nucleus_prim_index(inucl)+1, nucleus_prim_index(inucl+1) + do iwalk = 1, walk_num + do ielec = 1, elec_num + x = elec_coord(ielec,1,iwalk) - nucl_coord(inucl,1) + y = elec_coord(ielec,2,iwalk) - nucl_coord(inucl,2) + z = elec_coord(ielec,3,iwalk) - nucl_coord(inucl,3) + + r2 = x*x + y*y + z*z + ar2 = expo(iprim)*r2 + v = dexp(-ar2) + two_a = -2.d0 * expo(iprim) * v + + primitive_vgl(ielec, iwalk, 1, iprim) = v + primitive_vgl(ielec, iwalk, 2, iprim) = two_a * x + primitive_vgl(ielec, iwalk, 3, iprim) = two_a * y + primitive_vgl(ielec, iwalk, 4, iprim) = two_a * z + primitive_vgl(ielec, iwalk, 5, iprim) = two_a * (3.d0 - 2.d0*ar2) + + end do + end do + end do + end do + +end function qmckl_compute_ao_basis_primitive_gaussian_vgl_f + #+end_src + + #+begin_src c :tangle (eval h_private_func) :comments org :exports none +qmckl_exit_code qmckl_compute_ao_basis_primitive_gaussian_vgl( + const qmckl_context context, + const int64_t prim_num, + const int64_t elec_num, + const int64_t nucl_num, + const int64_t walk_num, + const int64_t* nucleus_prim_index, + const double* elec_coord, + const double* nucl_coord, + const double* expo, + double* const primitive_vgl); + #+end_src + + #+CALL: generate_c_interface(table=qmckl_ao_basis_primitive_gaussian_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_ao_basis_primitive_gaussian_vgl")) + + #+RESULTS: + #+begin_src f90 :tangle (eval f) :comments org :exports none + integer(c_int32_t) function qmckl_compute_ao_basis_primitive_gaussian_vgl & + (context, prim_num, elec_num, nucl_num, walk_num, nucleus_prim_index, elec_coord, nucl_coord, expo, primitive_vgl) & + 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 :: prim_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 :: walk_num + integer (c_int64_t) , intent(in) :: nucleus_prim_index(nucl_num) + real (c_double ) , intent(in) :: elec_coord(elec_num,3,walk_num) + real (c_double ) , intent(in) :: nucl_coord(elec_num,3) + real (c_double ) , intent(in) :: expo(prim_num) + real (c_double ) , intent(out) :: primitive_vgl(elec_num,walk_num,5,prim_num) + + integer(c_int32_t), external :: qmckl_compute_ao_basis_primitive_gaussian_vgl_f + info = qmckl_compute_ao_basis_primitive_gaussian_vgl_f & + (context, prim_num, elec_num, nucl_num, walk_num, nucleus_prim_index, elec_coord, nucl_coord, expo, primitive_vgl) + + end function qmckl_compute_ao_basis_primitive_gaussian_vgl + #+end_src + + #+begin_src python :results output :exports none +import numpy as np + +def f(a,x,y): + return np.exp( -a*(np.linalg.norm(x-y))**2 ) + +def df(a,x,y,n): + h0 = 1.e-6 + if n == 1: h = np.array([h0,0.,0.]) + elif n == 2: h = np.array([0.,h0,0.]) + elif n == 3: h = np.array([0.,0.,h0]) + return ( f(a,x+h,y) - f(a,x-h,y) ) / (2.*h0) + +def d2f(a,x,y,n): + h0 = 1.e-6 + if n == 1: h = np.array([h0,0.,0.]) + elif n == 2: h = np.array([0.,h0,0.]) + elif n == 3: h = np.array([0.,0.,h0]) + return ( f(a,x+h,y) - 2.*f(a,x,y) + f(a,x-h,y) ) / h0**2 + +def lf(a,x,y): + return d2f(a,x,y,1) + d2f(a,x,y,2) + d2f(a,x,y,3) + +elec_26_w1 = np.array( [ 1.49050402641, 2.90106987953, -1.05920815468 ] ) +elec_15_w2 = np.array( [ -2.20180344582,-1.9113150239, 2.2193744778600002 ] ) +nucl_1 = np.array( [ 1.096243353458458e+00, 8.907054016973815e-01, 7.777092280258892e-01 ] ) +nucl_2 = np.array( [ 1.168459237342663e+00, 1.125660720053393e+00, 2.833370314829343e+00 ] ) + +#double prim_vgl[prim_num][5][walk_num][elec_num]; +a = 0.9059; x = elec_26_w1 ; y = nucl_1 +print ( "[7][0][0][26] : %e"% f(a,x,y)) +print ( "[7][1][0][26] : %e"% df(a,x,y,1)) +print ( "[7][2][0][26] : %e"% df(a,x,y,2)) +print ( "[7][3][0][26] : %e"% df(a,x,y,3)) +print ( "[7][4][0][26] : %e"% lf(a,x,y)) + +a = 0.32578; x = elec_15_w2 ; y = nucl_2 +print ( "[39][0][1][15] : %e"% f(a,x,y)) +print ( "[39][1][1][15] : %e"% df(a,x,y,1)) +print ( "[39][2][1][15] : %e"% df(a,x,y,2)) +print ( "[39][3][1][15] : %e"% df(a,x,y,3)) +print ( "[39][4][1][15] : %e"% lf(a,x,y)) + + #+end_src + + #+RESULTS: + #+begin_example + [7][0][0][26] : 1.050157e-03 + [7][1][0][26] : -7.501497e-04 + [7][2][0][26] : -3.825069e-03 + [7][3][0][26] : 3.495056e-03 + [7][4][0][26] : 2.040013e-02 + [39][0][1][15] : 1.083038e-03 + [39][1][1][15] : 2.378275e-03 + [39][2][1][15] : 2.143086e-03 + [39][3][1][15] : 4.332750e-04 + [39][4][1][15] : 7.514605e-03 + #+end_example + +*** Test + + #+begin_src c :tangle (eval c_test) +#define walk_num chbrclf_walk_num +#define elec_num chbrclf_elec_num +#define prim_num chbrclf_prim_num + +int64_t elec_up_num = chbrclf_elec_up_num; +int64_t elec_dn_num = chbrclf_elec_dn_num; +double* elec_coord = &(chbrclf_elec_coord[0][0][0]); + +rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num); +assert (rc == QMCKL_SUCCESS); + +rc = qmckl_set_electron_walk_num (context, walk_num); +assert (rc == QMCKL_SUCCESS); + +assert(qmckl_electron_provided(context)); + +rc = qmckl_set_electron_coord (context, 'N', elec_coord); +assert(rc == QMCKL_SUCCESS); + + + +double prim_vgl[prim_num][5][walk_num][elec_num]; + +rc = qmckl_get_ao_basis_primitive_vgl(context, &(prim_vgl[0][0][0][0])); +assert (rc == QMCKL_SUCCESS); + +assert( fabs(prim_vgl[7][0][0][26] - ( 1.0501570432064878E-003)) < 1.e-14 ); +assert( fabs(prim_vgl[7][1][0][26] - (-7.5014974095310560E-004)) < 1.e-14 ); +assert( fabs(prim_vgl[7][2][0][26] - (-3.8250692897610380E-003)) < 1.e-14 ); +assert( fabs(prim_vgl[7][3][0][26] - ( 3.4950559194080275E-003)) < 1.e-14 ); +assert( fabs(prim_vgl[7][4][0][26] - ( 2.0392163767356572E-002)) < 1.e-14 ); + +assert( fabs(prim_vgl[39][0][1][15] - ( 1.0825844173157661E-003)) < 1.e-14 ); +assert( fabs(prim_vgl[39][1][1][15] - ( 2.3774237611651531E-003)) < 1.e-14 ); +assert( fabs(prim_vgl[39][2][1][15] - ( 2.1423191526963063E-003)) < 1.e-14 ); +assert( fabs(prim_vgl[39][3][1][15] - ( 4.3312003523048492E-004)) < 1.e-14 ); +assert( fabs(prim_vgl[39][4][1][15] - ( 7.5174404780004771E-003)) < 1.e-14 ); + + + + + #+end_src * Polynomial part -** Powers of $x-X_i$ +** General functions for Powers of $x-X_i$ :PROPERTIES: :Name: qmckl_ao_power :CRetType: qmckl_exit_code @@ -1238,7 +1872,7 @@ end function test_qmckl_ao_power assert(0 == test_qmckl_ao_power(context)); #+end_src -** Value, Gradient and Laplacian of a polynomial +** General functions for Value, Gradient and Laplacian of a polynomial :PROPERTIES: :Name: qmckl_ao_polynomial_vgl :CRetType: qmckl_exit_code @@ -1614,220 +2248,6 @@ end function test_qmckl_ao_polynomial_vgl assert(0 == test_qmckl_ao_polynomial_vgl(context)); #+end_src -* Radial part -** Gaussian basis functions - - ~qmckl_ao_gaussian_vgl~ computes the values, gradients and - Laplacians at a given point of ~n~ Gaussian functions centered at - the same point: - - \[ v_i = \exp(-a_i |X-R|^2) \] - \[ \nabla_x v_i = -2 a_i (X_x - R_x) v_i \] - \[ \nabla_y v_i = -2 a_i (X_y - R_y) v_i \] - \[ \nabla_z v_i = -2 a_i (X_z - R_z) v_i \] - \[ \Delta v_i = a_i (4 |X-R|^2 a_i - 6) v_i \] - - |--------------+--------+------------------------------------------------------| - | ~context~ | input | Global state | - | ~X(3)~ | input | Array containing the coordinates of the points | - | ~R(3)~ | input | Array containing the x,y,z coordinates of the center | - | ~n~ | input | Number of computed Gaussians | - | ~A(n)~ | input | Exponents of the Gaussians | - | ~VGL(ldv,5)~ | output | Value, gradients and Laplacian of the Gaussians | - | ~ldv~ | input | Leading dimension of array ~VGL~ | - |--------------+--------+------------------------------------------------------| - - Requirements - - - ~context~ is not 0 - - ~n~ > 0 - - ~ldv~ >= 5 - - ~A(i)~ > 0 for all ~i~ - - ~X~ is allocated with at least $3 \times 8$ bytes - - ~R~ is allocated with at least $3 \times 8$ bytes - - ~A~ is allocated with at least $n \times 8$ bytes - - ~VGL~ is allocated with at least $n \times 5 \times 8$ bytes - - #+begin_src c :tangle (eval h_func) -qmckl_exit_code -qmckl_ao_gaussian_vgl(const qmckl_context context, - const double *X, - const double *R, - const int64_t *n, - const int64_t *A, - const double *VGL, - const int64_t ldv); - #+end_src - - #+begin_src f90 :tangle (eval f) -integer function qmckl_ao_gaussian_vgl_f(context, X, R, n, A, VGL, ldv) result(info) - use qmckl - implicit none - integer*8 , intent(in) :: context - real*8 , intent(in) :: X(3), R(3) - integer*8 , intent(in) :: n - real*8 , intent(in) :: A(n) - real*8 , intent(out) :: VGL(ldv,5) - integer*8 , intent(in) :: ldv - - integer*8 :: i,j - real*8 :: Y(3), r2, t, u, v - - info = QMCKL_SUCCESS - - if (context == QMCKL_NULL_CONTEXT) then - info = QMCKL_INVALID_CONTEXT - return - endif - - if (n <= 0) then - info = QMCKL_INVALID_ARG_4 - return - endif - - if (ldv < n) then - info = QMCKL_INVALID_ARG_7 - return - endif - - - do i=1,3 - Y(i) = X(i) - R(i) - end do - r2 = Y(1)*Y(1) + Y(2)*Y(2) + Y(3)*Y(3) - - do i=1,n - VGL(i,1) = dexp(-A(i) * r2) - end do - - do i=1,n - VGL(i,5) = A(i) * VGL(i,1) - end do - - t = -2.d0 * ( X(1) - R(1) ) - u = -2.d0 * ( X(2) - R(2) ) - v = -2.d0 * ( X(3) - R(3) ) - - do i=1,n - VGL(i,2) = t * VGL(i,5) - VGL(i,3) = u * VGL(i,5) - VGL(i,4) = v * VGL(i,5) - end do - - t = 4.d0 * r2 - do i=1,n - VGL(i,5) = (t * A(i) - 6.d0) * VGL(i,5) - end do - -end function qmckl_ao_gaussian_vgl_f - #+end_src - - #+begin_src f90 :tangle (eval f) :exports none -integer(c_int32_t) function qmckl_ao_gaussian_vgl(context, X, R, n, A, VGL, ldv) & - bind(C) result(info) - use, intrinsic :: iso_c_binding - implicit none - integer (c_int64_t) , intent(in) , value :: context - real (c_double) , intent(in) :: X(3), R(3) - integer (c_int64_t) , intent(in) , value :: n - real (c_double) , intent(in) :: A(n) - real (c_double) , intent(out) :: VGL(ldv,5) - integer (c_int64_t) , intent(in) , value :: ldv - - integer, external :: qmckl_ao_gaussian_vgl_f - info = qmckl_ao_gaussian_vgl_f(context, X, R, n, A, VGL, ldv) -end function qmckl_ao_gaussian_vgl - #+end_src - - #+begin_src f90 :tangle (eval fh_func) :exports none - interface - integer(c_int32_t) function qmckl_ao_gaussian_vgl(context, X, R, n, A, VGL, ldv) & - bind(C) - use, intrinsic :: iso_c_binding - integer (c_int64_t) , intent(in) , value :: context - integer (c_int64_t) , intent(in) , value :: ldv - integer (c_int64_t) , intent(in) , value :: n - real (c_double) , intent(in) :: X(3), R(3), A(n) - real (c_double) , intent(out) :: VGL(ldv,5) - end function qmckl_ao_gaussian_vgl - end interface - #+end_src - - # Test - #+begin_src f90 :tangle (eval f_test) -integer(c_int32_t) function test_qmckl_ao_gaussian_vgl(context) bind(C) - use qmckl - implicit none - - integer(c_int64_t), intent(in), value :: context - - integer*8 :: n, ldv, j, i - double precision :: X(3), R(3), Y(3), r2 - double precision, allocatable :: VGL(:,:), A(:) - double precision :: epsilon - - epsilon = qmckl_get_numprec_epsilon(context) - - X = (/ 1.1 , 2.2 , 3.3 /) - R = (/ 0.1 , 1.2 , -2.3 /) - Y(:) = X(:) - R(:) - r2 = Y(1)**2 + Y(2)**2 + Y(3)**2 - - n = 10; - ldv = 100; - - allocate (A(n), VGL(ldv,5)) - do i=1,n - A(i) = 0.0013 * dble(ishft(1,i)) - end do - - - test_qmckl_ao_gaussian_vgl = & - qmckl_ao_gaussian_vgl(context, X, R, n, A, VGL, ldv) - if (test_qmckl_ao_gaussian_vgl /= 0) return - - test_qmckl_ao_gaussian_vgl = -1 - - do i=1,n - test_qmckl_ao_gaussian_vgl = -11 - if (dabs(1.d0 - VGL(i,1) / (& - dexp(-A(i) * r2) & - )) > epsilon ) return - - test_qmckl_ao_gaussian_vgl = -12 - if (dabs(1.d0 - VGL(i,2) / (& - -2.d0 * A(i) * Y(1) * dexp(-A(i) * r2) & - )) > epsilon ) return - - test_qmckl_ao_gaussian_vgl = -13 - if (dabs(1.d0 - VGL(i,3) / (& - -2.d0 * A(i) * Y(2) * dexp(-A(i) * r2) & - )) > epsilon ) return - - test_qmckl_ao_gaussian_vgl = -14 - if (dabs(1.d0 - VGL(i,4) / (& - -2.d0 * A(i) * Y(3) * dexp(-A(i) * r2) & - )) > epsilon ) return - - test_qmckl_ao_gaussian_vgl = -15 - if (dabs(1.d0 - VGL(i,5) / (& - A(i) * (4.d0*r2*A(i) - 6.d0) * dexp(-A(i) * r2) & - )) > epsilon ) return - end do - - test_qmckl_ao_gaussian_vgl = 0 - - deallocate(VGL) -end function test_qmckl_ao_gaussian_vgl - #+end_src - - #+begin_src c :tangle (eval c_test) :exports none - int test_qmckl_ao_gaussian_vgl(qmckl_context context); - assert(0 == test_qmckl_ao_gaussian_vgl(context)); - #+end_src - -** TODO Slater basis functions -** TODO Radial functions on a grid * Combining radial and polynomial parts * End of files :noexport: diff --git a/org/qmckl_context.org b/org/qmckl_context.org index 6c7a3ac..11eb235 100644 --- a/org/qmckl_context.org +++ b/org/qmckl_context.org @@ -31,6 +31,9 @@ int main() { #include "qmckl_nucleus_private_type.h" #include "qmckl_electron_private_type.h" #include "qmckl_ao_private_type.h" +#include "qmckl_nucleus_private_func.h" +#include "qmckl_electron_private_func.h" +#include "qmckl_ao_private_func.h" #+end_src #+begin_src c :tangle (eval c) @@ -213,17 +216,26 @@ qmckl_context qmckl_context_create() { } /* Initialize data */ - ctx->tag = VALID_TAG; + { + ctx->tag = VALID_TAG; + + const qmckl_context context = (const qmckl_context) ctx; + assert ( qmckl_context_check(context) != QMCKL_NULL_CONTEXT ); + + qmckl_exit_code rc; + + ctx->numprec.precision = QMCKL_DEFAULT_PRECISION; + ctx->numprec.range = QMCKL_DEFAULT_RANGE; - const qmckl_context context = (const qmckl_context) ctx; - assert ( qmckl_context_check(context) != QMCKL_NULL_CONTEXT ); - - ctx->numprec.precision = QMCKL_DEFAULT_PRECISION; - ctx->numprec.range = QMCKL_DEFAULT_RANGE; - - ctx->ao_basis.uninitialized = (1 << 12) - 1; - ctx->nucleus.uninitialized = (1 << 4) - 1; - ctx->electron.uninitialized = (1 << 3) - 1; + rc = qmckl_init_electron(context); + assert (rc == QMCKL_SUCCESS); + + rc = qmckl_init_nucleus(context); + assert (rc == QMCKL_SUCCESS); + + rc = qmckl_init_ao_basis(context); + assert (rc == QMCKL_SUCCESS); + } /* Allocate qmckl_memory_struct */ { diff --git a/org/qmckl_electron.org b/org/qmckl_electron.org index 11f025f..9cdd100 100644 --- a/org/qmckl_electron.org +++ b/org/qmckl_electron.org @@ -63,25 +63,25 @@ 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 distances | - | ~ee_distance_rescaled_date~ | uint64_t | Last modification date of the electron-electron distances | - | ~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 | + | ~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 distances | + | ~ee_distance_rescaled_date~ | ~uint64_t~ | Last modification date of the electron-electron distances | + | ~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 | ** Data structure @@ -107,16 +107,41 @@ typedef struct qmckl_electron_struct { int32_t uninitialized; bool provided; } qmckl_electron_struct; + #+end_src The ~uninitialized~ integer contains one bit set to one for each - initialization function which has not bee called. It becomes equal + 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. - When all the data relative to electrons have been set, the - following function returns ~true~. + #+begin_src c :comments org :tangle (eval h_private_func) +qmckl_exit_code qmckl_init_electron(qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) +qmckl_exit_code qmckl_init_electron(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->electron.uninitialized = (1 << 2) - 1; + + /* Default values */ + ctx->electron.rescale_factor_kappa_ee = 1.0; + ctx->electron.rescale_factor_kappa_en = 1.0; + + return QMCKL_SUCCESS; +} + #+end_src + + #+begin_src c :comments org :tangle (eval h_func) bool qmckl_electron_provided (const qmckl_context context); #+end_src @@ -287,13 +312,13 @@ qmckl_get_electron_walk_num (const qmckl_context context, int64_t* const walk_nu *** Scaling factors Kappa #+begin_src c :comments org :tangle (eval h_func) :exports none -qmckl_exit_code qmckl_get_kappa_ee (const qmckl_context context, double* const rescale_factor_kappa_ee); -qmckl_exit_code qmckl_get_kappa_en (const qmckl_context context, double* const rescale_factor_kappa_en); +qmckl_exit_code qmckl_get_electron_rescale_factor_ee (const qmckl_context context, double* const rescale_factor_kappa_ee); +qmckl_exit_code qmckl_get_electron_rescale_factor_en (const qmckl_context context, double* const rescale_factor_kappa_en); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code -qmckl_get_kappa_ee (const qmckl_context context, double* const rescale_factor_kappa_ee) { +qmckl_get_electron_rescale_factor_ee (const qmckl_context context, double* const rescale_factor_kappa_ee) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_INVALID_CONTEXT; @@ -302,27 +327,22 @@ qmckl_get_kappa_ee (const qmckl_context context, double* const rescale_factor_ka if (rescale_factor_kappa_ee == NULL) { return qmckl_failwith( context, QMCKL_INVALID_ARG_2, - "qmckl_get_kappa_ee", + "qmckl_get_electron_rescale_factor_ee", "rescale_factor_kappa_ee is a null pointer"); } qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - int32_t mask = 1 << 2; + assert (ctx->electron.rescale_factor_kappa_ee > 0.0); - if ( (ctx->electron.uninitialized & mask) != 0) { - return QMCKL_NOT_PROVIDED; - } - - // TODO: assert (ctx->electron.rescale_factor_kappa_ee > (double) 0); - ,*rescale_factor_kappa_ee = ctx->electron.rescale_factor_kappa_ee; + *rescale_factor_kappa_ee = ctx->electron.rescale_factor_kappa_ee; return QMCKL_SUCCESS; } qmckl_exit_code -qmckl_get_kappa_en (const qmckl_context context, double* const rescale_factor_kappa_en) { +qmckl_get_electron_rescale_factor_en (const qmckl_context context, double* const rescale_factor_kappa_en) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_INVALID_CONTEXT; } @@ -330,21 +350,15 @@ qmckl_get_kappa_en (const qmckl_context context, double* const rescale_factor_ka if (rescale_factor_kappa_en == NULL) { return qmckl_failwith( context, QMCKL_INVALID_ARG_2, - "qmckl_get_kappa_en", + "qmckl_get_electron_rescale_factor_en", "rescale_factor_kappa_en is a null pointer"); } qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - int32_t mask = 1 << 2; - - if ( (ctx->electron.uninitialized & mask) != 0) { - return QMCKL_NOT_PROVIDED; - } - - // TODO: assert (ctx->electron.rescale_factor_kappa_en > (double) 0); - ,*rescale_factor_kappa_en = ctx->electron.rescale_factor_kappa_en; + assert (ctx->electron.rescale_factor_kappa_en > 0.0); + *rescale_factor_kappa_en = ctx->electron.rescale_factor_kappa_en; return QMCKL_SUCCESS; } #+end_src @@ -435,10 +449,12 @@ qmckl_get_electron_coord (const qmckl_context context, const char transp, double both allocated. #+begin_src c :comments org :tangle (eval h_func) -qmckl_exit_code qmckl_set_electron_num (qmckl_context context, const int64_t up_num, const int64_t down_num); -qmckl_exit_code qmckl_set_kappa (qmckl_context context, const double rescale_factor_kappa_ee, const double rescale_factor_kappa_en); -qmckl_exit_code qmckl_set_electron_walk_num (qmckl_context context, const int64_t walk_num); -qmckl_exit_code qmckl_set_electron_coord (qmckl_context context, const char transp, const double* coord); +qmckl_exit_code qmckl_set_electron_num (qmckl_context context, const int64_t up_num, const int64_t down_num); +qmckl_exit_code qmckl_set_electron_walk_num (qmckl_context context, const int64_t walk_num); +qmckl_exit_code qmckl_set_electron_coord (qmckl_context context, const char transp, const double* coord); + +qmckl_exit_code qmckl_set_electron_rescale_factor_ee (qmckl_context context, const double kappa_ee); +qmckl_exit_code qmckl_set_electron_rescale_factor_en (qmckl_context context, const double kappa_en); #+end_src #+NAME:pre2 @@ -551,32 +567,37 @@ qmckl_set_electron_walk_num(qmckl_context context, const int64_t walk_num) { #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code -qmckl_set_kappa(qmckl_context context, - const double rescale_factor_kappa_ee, +qmckl_set_electron_rescale_factor_ee(qmckl_context context, + const double rescale_factor_kappa_ee) { + <> + + if (rescale_factor_kappa_ee <= 0.0) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_set_electron_rescale_factor_ee", + "rescale_factor_kappa_ee <= 0.0"); + } + + ctx->electron.rescale_factor_kappa_ee = rescale_factor_kappa_ee; + + return QMCKL_SUCCESS; +} + +qmckl_exit_code +qmckl_set_electron_rescale_factor_en(qmckl_context context, const double rescale_factor_kappa_en) { <> - // TODO: Check for 0 values - //if (rescale_factor_kappa_ee != 0) { - // return qmckl_failwith( context, - // QMCKL_INVALID_ARG_2, - // "qmckl_set_kappa", - // "rescale_factor_kappa_ee == 0"); - //} + if (rescale_factor_kappa_en <= 0.0) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_set_electron_rescale_factor_en", + "rescale_factor_kappa_en <= 0.0"); + } - //if (rescale_factor_kappa_en <= 0) { - // return qmckl_failwith( context, - // QMCKL_INVALID_ARG_3, - // "qmckl_set_kappa", - // "rescale_factor_kappa_en == 0"); - //} - - int32_t mask = 1 << 2; - - ctx->electron.rescale_factor_kappa_ee = rescale_factor_kappa_ee; ctx->electron.rescale_factor_kappa_en = rescale_factor_kappa_en; - <> + return QMCKL_SUCCESS; } #+end_src @@ -680,14 +701,13 @@ int64_t walk_num = chbrclf_walk_num; int64_t elec_num = chbrclf_elec_num; int64_t elec_up_num = chbrclf_elec_up_num; int64_t elec_dn_num = chbrclf_elec_dn_num; -double rescale_factor_kappa_ee = 1.0; // TODO Get rescale_factor_kappa_ee from chbrclf -double rescale_factor_kappa_en = 1.0; // TODO Get rescale_factor_kappa_en from chbrclf +double rescale_factor_kappa_ee = 2.0; +double rescale_factor_kappa_en = 3.0; double* elec_coord = &(chbrclf_elec_coord[0][0][0]); int64_t nucl_num = chbrclf_nucl_num; double* charge = chbrclf_charge; double* nucl_coord = &(chbrclf_nucl_coord[0][0]); -double nucl_rescale_factor_kappa = 1.0; // TODO Change get rescale_factor_kappa from chbrclf example /* --- */ @@ -722,23 +742,27 @@ rc = qmckl_get_electron_num (context, &n); assert(rc == QMCKL_SUCCESS); assert(n == elec_num); -double k_ee; -double k_en; -rc = qmckl_get_kappa_ee (context, &k_ee); -assert(rc == QMCKL_NOT_PROVIDED); - -rc = qmckl_get_kappa_en (context, &k_en); -assert(rc == QMCKL_NOT_PROVIDED); - -rc = qmckl_set_kappa (context, rescale_factor_kappa_ee, rescale_factor_kappa_en); +double k_ee = 0.; +double k_en = 0.; +rc = qmckl_get_electron_rescale_factor_ee (context, &k_ee); assert(rc == QMCKL_SUCCESS); -assert(!qmckl_electron_provided(context)); +assert(k_ee == 1.0); -rc = qmckl_get_kappa_ee (context, &k_ee); +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_kappa_en (context, &k_en); +rc = qmckl_get_electron_rescale_factor_en (context, &k_en); assert(rc == QMCKL_SUCCESS); assert(k_en == rescale_factor_kappa_en); @@ -1294,7 +1318,10 @@ qmckl_exit_code qmckl_provide_en_distance(qmckl_context context) assert (ctx != NULL); if (!(ctx->nucleus.provided)) { - return QMCKL_NOT_PROVIDED; + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_provide_en_distance", + NULL); } /* Compute if necessary */ @@ -1475,9 +1502,6 @@ assert(qmckl_electron_provided(context)); rc = qmckl_set_nucleus_num (context, nucl_num); assert(rc == QMCKL_SUCCESS); -rc = qmckl_set_nucleus_kappa (context, nucl_rescale_factor_kappa); -assert(rc == QMCKL_SUCCESS); - rc = qmckl_set_nucleus_charge (context, charge); assert (rc == QMCKL_SUCCESS); @@ -1753,9 +1777,6 @@ assert(qmckl_electron_provided(context)); rc = qmckl_set_nucleus_num (context, nucl_num); assert(rc == QMCKL_SUCCESS); -rc = qmckl_set_nucleus_kappa (context, nucl_rescale_factor_kappa); -assert(rc == QMCKL_SUCCESS); - rc = qmckl_set_nucleus_charge (context, charge); assert (rc == QMCKL_SUCCESS); diff --git a/org/qmckl_nucleus.org b/org/qmckl_nucleus.org index 004560f..99237af 100644 --- a/org/qmckl_nucleus.org +++ b/org/qmckl_nucleus.org @@ -92,20 +92,49 @@ typedef struct qmckl_nucleus_struct { int32_t uninitialized; bool provided; } qmckl_nucleus_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_private_func) +qmckl_exit_code qmckl_init_nucleus(qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) +qmckl_exit_code qmckl_init_nucleus(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->nucleus.uninitialized = (1 << 3) - 1; + + /* Default values */ + ctx->nucleus.rescale_factor_kappa = 1.0; + + return QMCKL_SUCCESS; +} + #+end_src + + + ** Access functions #+begin_src c :comments org :tangle (eval h_func) :exports none -qmckl_exit_code qmckl_get_nucleus_num (const qmckl_context context, int64_t* const num); -qmckl_exit_code qmckl_get_nucleus_charge (const qmckl_context context, double* const charge); -qmckl_exit_code qmckl_get_nucleus_kappa (const qmckl_context context, double* const rescale_factor_kappa); -qmckl_exit_code qmckl_get_nucleus_coord (const qmckl_context context, const char transp, double* const coord); +qmckl_exit_code qmckl_get_nucleus_num (const qmckl_context context, int64_t* const num); +qmckl_exit_code qmckl_get_nucleus_charge (const qmckl_context context, double* const charge); +qmckl_exit_code qmckl_get_nucleus_coord (const qmckl_context context, const char transp, double* const coord); + +qmckl_exit_code qmckl_get_nucleus_rescale_factor (const qmckl_context context, double* const rescale_factor_kappa); #+end_src #+NAME:post @@ -136,7 +165,7 @@ qmckl_get_nucleus_num (const qmckl_context context, int64_t* const num) { int32_t mask = 1 << 0; if ( (ctx->nucleus.uninitialized & mask) != 0) { - *num = (int64_t) 0; + ,*num = (int64_t) 0; return qmckl_failwith( context, QMCKL_NOT_PROVIDED, "qmckl_get_nucleus_num", @@ -187,23 +216,25 @@ qmckl_get_nucleus_charge (const qmckl_context context, double* const charge) { qmckl_exit_code -qmckl_get_nucleus_kappa (const qmckl_context context, double* const rescale_factor_kappa) { +qmckl_get_nucleus_rescale_factor (const qmckl_context context, + double* const rescale_factor_kappa) +{ if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_INVALID_CONTEXT; } + if (rescale_factor_kappa == NULL) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_get_nucleus_rescale_factor", + "rescale_factor_kappa is a null pointer"); + } + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - int32_t mask = 1 << 2; - - if ( (ctx->nucleus.uninitialized & mask) != 0) { - return qmckl_failwith( context, - QMCKL_NOT_PROVIDED, - "qmckl_get_nucleus_kappa", - "nucleus data is not provided"); - } + assert (ctx->nucleus.rescale_factor_kappa > 0.0); (*rescale_factor_kappa) = ctx->nucleus.rescale_factor_kappa; @@ -235,7 +266,7 @@ qmckl_get_nucleus_coord (const qmckl_context context, const char transp, double* qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - int32_t mask = 1 << 3; + int32_t mask = 1 << 2; if ( (ctx->nucleus.uninitialized & mask) != 0) { return qmckl_failwith( context, @@ -295,8 +326,9 @@ bool qmckl_nucleus_provided(const qmckl_context context) { #+begin_src c :comments org :tangle (eval h_func) qmckl_exit_code qmckl_set_nucleus_num (qmckl_context context, const int64_t num); qmckl_exit_code qmckl_set_nucleus_charge (qmckl_context context, const double* charge); -qmckl_exit_code qmckl_set_nucleus_kappa (qmckl_context context, const double rescale_factor_kappa); qmckl_exit_code qmckl_set_nucleus_coord (qmckl_context context, const char transp, const double* coord); + +qmckl_exit_code qmckl_set_nucleus_rescale_factor (qmckl_context context, const double rescale_factor_kappa); #+end_src #+NAME:pre2 @@ -331,7 +363,7 @@ qmckl_set_nucleus_num(qmckl_context context, const int64_t num) { "num <= 0"); } - int32_t mask = 1; + int32_t mask = 1 << 0; ctx->nucleus.num = num; @@ -388,22 +420,19 @@ qmckl_set_nucleus_charge(qmckl_context context, const double* charge) { #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code -qmckl_set_nucleus_kappa(qmckl_context context, const double rescale_factor_kappa) { +qmckl_set_nucleus_rescale_factor(qmckl_context context, const double rescale_factor_kappa) { <> - //TODO: Check for small values of kappa if (rescale_factor_kappa <= 0.0) { return qmckl_failwith( context, QMCKL_INVALID_ARG_2, "qmckl_set_nucleus_kappa", - "rescale_factor_kappa cannot be 0"); + "rescale_factor_kappa cannot be <= 0."); } - int32_t mask = 1 << 2; - ctx->nucleus.rescale_factor_kappa = rescale_factor_kappa; - <> + return QMCKL_SUCCESS; } #+end_src @@ -418,7 +447,7 @@ qmckl_set_nucleus_coord(qmckl_context context, const char transp, const double* int64_t nucl_num = (int64_t) 0; qmckl_exit_code rc; - int32_t mask = 1 << 3; + int32_t mask = 1 << 2; rc = qmckl_get_nucleus_num(context, &nucl_num); if (rc != QMCKL_SUCCESS) return rc; @@ -459,7 +488,7 @@ qmckl_set_nucleus_coord(qmckl_context context, const char transp, const double* const int64_t nucl_num = chbrclf_nucl_num; const double* nucl_charge = chbrclf_charge; const double* nucl_coord = &(chbrclf_nucl_coord[0][0]); -const double nucl_rescale_factor_kappa = 1.0; // TODO Change get rescale_factor_kappa from chbrclf example +const double nucl_rescale_factor_kappa = 2.0; /* --- */ @@ -481,15 +510,15 @@ assert(rc == QMCKL_SUCCESS); assert(n == nucl_num); double k; -rc = qmckl_get_nucleus_kappa (context, &k); -assert(rc == QMCKL_NOT_PROVIDED); - - -rc = qmckl_set_nucleus_kappa (context, nucl_rescale_factor_kappa); +rc = qmckl_get_nucleus_rescale_factor (context, &k); assert(rc == QMCKL_SUCCESS); -assert(!qmckl_nucleus_provided(context)); +assert(k == 1.0); -rc = qmckl_get_nucleus_kappa (context, &k); + +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); @@ -532,7 +561,7 @@ for (size_t i=0 ; i Date: Mon, 21 Jun 2021 11:18:52 +0200 Subject: [PATCH 8/9] Fix broken documentation link --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 1fa2d3c..095b50c 100644 --- a/README.md +++ b/README.md @@ -9,7 +9,7 @@ simple language and provide a standard API and tests to enable the development of high-performance QMCkl implementations taking advantage of modern hardware. -See the [source code](https://github.com/TREX-CoE/qmckl/tree/main/src) +See the [source code](https://github.com/TREX-CoE/qmckl/blob/master/org/qmckl.org) to read the documentation. From ddeb5f82f408ecdfdf3de1c6a9c8792b7cea7476 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 21 Jun 2021 11:20:34 +0200 Subject: [PATCH 9/9] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 095b50c..97606c1 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ # QMCkl: Quantum Monte Carlo Kernel Library -![Build Status](https://github.com/TREX-CoE/qmckl/workflows/test-build/badge.svg?branch=main) +![Build Status](https://github.com/TREX-CoE/qmckl/workflows/test-build/badge.svg?branch=master) The domain of quantum chemistry needs a library in which the main kernels of Quantum Monte Carlo (QMC) methods are implemented. In the