diff --git a/.github/workflows/test-build.yml b/.github/workflows/test-build.yml index 543abf8..4bfb193 100644 --- a/.github/workflows/test-build.yml +++ b/.github/workflows/test-build.yml @@ -2,9 +2,7 @@ name: test-build on: push: - branches: [ master ] pull_request: - branches: [ master ] jobs: x86_ubuntu: diff --git a/configure.ac b/configure.ac index 31322e1..ca736be 100644 --- a/configure.ac +++ b/configure.ac @@ -93,6 +93,7 @@ AC_PROG_F77 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 +AM_PROG_CC_C_O AC_PROG_FC AC_PROG_FC_C_O AC_FC_PP_DEFINE @@ -137,10 +138,10 @@ case "$with_chameleon" in [PKG_CFLAGS="$PKG_CFLAGS $LIBCHAMELEON_CFLAGS" PKG_LIBS="$PKG_LIBS $LIBCHAMELEON_LIBS"] ,[ - + ## something went wrong. ## try to find the package without pkg-config - + ## check that the library is actually new enough. ## by testing for a 1.0.0+ function which we use AC_CHECK_LIB(chameleon,CHAMELEON_finalize,[LIBCHAMELEON_LIBS="-lchameleon"]) @@ -212,6 +213,9 @@ esac case $CC in + *gcc*) + CFLAGS="$CFLAGS -fPIC" + ;; *nvc*) CFLAGS="$CFLAGS -fPIC" ;; @@ -224,6 +228,131 @@ AS_IF([test "$HAVE_HPC" = "yes"], [ AC_DEFINE([HAVE_HPC], [1], [If defined, activate HPC routines]) ]) +# Enable Verificarlo tests +AC_ARG_ENABLE([vfc_ci], +[ --enable-vfc_ci Build the library with vfc_ci support], +[case "${enableval}" in + yes) vfc_ci=true && FCFLAGS="-D VFC_CI $FCFLAGS" && CFLAGS="-D VFC_CI $CFLAGS";; + no) vfc_ci=false ;; + *) AC_MSG_ERROR([bad value ${enableval} for --enable_vfc_ci]) ;; +esac],[vfc_ci=false]) +AM_CONDITIONAL([VFC_CI], [test x$vfc_ci = xtrue]) + +if test "$FC" = "verificarlo-f"; then + AC_MSG_NOTICE(verificarlo-f detected) + # Arguments order is important here + FCFLAGS="-Mpreprocess $FCFLAGS" +fi + +## Enable GPU offloading + +# GPU offloading +AC_ARG_ENABLE(gpu, [AS_HELP_STRING([--enable-gpu],[openmp|openacc : Use GPU-offloaded functions])], enable_gpu=$enableval, enable_gpu=no) +AS_IF([test "$enable_gpu" = "yes"], [enable_gpu="openmp"]) + +# OpenMP offloading +HAVE_OPENMP_OFFLOAD="no" +AS_IF([test "$enable_gpu" = "openmp"], [ + AC_DEFINE([HAVE_OPENMP_OFFLOAD], [1], [If defined, activate OpenMP-offloaded routines]) + HAVE_OPENMP_OFFLOAD="yes" + case $CC in + + *gcc*) + CFLAGS="$CFLAGS -fopenmp" + ;; + *nvc*) + CFLAGS="$CFLAGS -mp=gpu" + ;; + esac + + case $FC in + + *gfortran*) + FCFLAGS="$FCFLAGS -fopenmp" + ;; + *nvfortran*) + FCFLAGS="$FCFLAGS -mp=gpu" + ;; + esac] +) + +# OpenMP offloading +HAVE_OPENACC_OFFLOAD="no" +AS_IF([test "$enable_gpu" = "openacc"], [ + AC_DEFINE([HAVE_OPENACC_OFFLOAD], [1], [If defined, activate OpenACC-offloaded routines]) + HAVE_OPENACC_OFFLOAD="yes" + case $CC in + + *gcc*) + CFLAGS="$CFLAGS -fopenacc" + ;; + *nvc*) + CFLAGS="$CFLAGS -acc=gpu" + ;; + esac + + case $FC in + + *gfortran*) + FCFLAGS="$FCFLAGS -fopenacc" + ;; + *nvfortran*) + FCFLAGS="$FCFLAGS -acc=gpu" + ;; + esac + +]) + +# cuBLAS offloading +AC_ARG_WITH(cublas, [AS_HELP_STRING([--with-cublas],[Use cuBLAS-offloaded functions])], HAVE_CUBLAS_OFFLOAD=$withval, HAVE_CUBLAS_OFFLOAD=no) +AS_IF([test "$HAVE_CUBLAS_OFFLOAD" = "yes"], [ + AC_DEFINE([HAVE_CUBLAS_OFFLOAD], [1], [If defined, activate cuBLAS-offloaded routines]) + HAVE_OPENACC_OFFLOAD="yes" + case $CC in + + *gcc*) + CFLAGS="$CFLAGS -fopenmp" + LDFLAGS="-lcublas" + ;; + *nvc*) + CFLAGS="$CFLAGS -mp=gpu -cudalib=cublas" + ;; + esac + + case $FC in + + *gfortran*) + FCFLAGS="$FCFLAGS -fopenmp" + ;; + *nvfortran*) + FCFLAGS="$FCFLAGS -mp=gpu -cudalib=cublas" + ;; + esac +]) + +AC_ARG_ENABLE(malloc-trace, [AS_HELP_STRING([--enable-malloc-trace],[use debug malloc/free])], ok=$enableval, ok=no) +if test "$ok" = "yes"; then + AC_DEFINE(MALLOC_TRACE,"malloc_trace.dat",[Define to use debugging malloc/free]) + ARGS="${ARGS} malloc-trace" +fi + +AC_ARG_ENABLE(prof, [AS_HELP_STRING([--enable-prof],[compile for profiling])], ok=$enableval, ok=no) +if test "$ok" = "yes"; then + CFLAGS="${CFLAGS} -pg" + AC_DEFINE(ENABLE_PROF,1,[Define when using the profiler tool]) + ARGS="${ARGS} prof" +fi + +AC_ARG_WITH(efence, [AS_HELP_STRING([--with-efence],[use ElectricFence library])], ok=$withval, ok=no) +if test "$ok" = "yes"; then + AC_CHECK_LIB([efence], [malloc]) + ARGS="${ARGS} efence" +fi + + + +## + AC_ARG_ENABLE(debug, [AS_HELP_STRING([--enable-debug],[compile for debugging])], ok=$enableval, ok=no) if test "$ok" = "yes"; then if test "$GCC" = "yes"; then @@ -247,26 +376,6 @@ if test "$ok" = "yes"; then ARGS="${ARGS} debug" fi -AC_ARG_ENABLE(malloc-trace, [AS_HELP_STRING([--enable-malloc-trace],[use debug malloc/free])], ok=$enableval, ok=no) -if test "$ok" = "yes"; then - AC_DEFINE(MALLOC_TRACE,"malloc_trace.dat",[Define to use debugging malloc/free]) - ARGS="${ARGS} malloc-trace" -fi - -AC_ARG_ENABLE(prof, [AS_HELP_STRING([--enable-prof],[compile for profiling])], ok=$enableval, ok=no) -if test "$ok" = "yes"; then - CFLAGS="${CFLAGS} -pg" - AC_DEFINE(ENABLE_PROF,1,[Define when using the profiler tool]) - ARGS="${ARGS} prof" -fi - -AC_ARG_WITH(efence, [AS_HELP_STRING([--with-efence],[use ElectricFence library])], ok=$withval, ok=no) -if test "$ok" = "yes"; then - AC_CHECK_LIB(efence, malloc) - ARGS="${ARGS} efence" -fi - - # Checks for header files. ## qmckl @@ -321,21 +430,6 @@ if test "x${QMCKL_DEVEL}" != "x"; then fi -# Enable Verificarlo tests -AC_ARG_ENABLE([vfc_ci], -[ --enable-vfc_ci Build the library with vfc_ci support], -[case "${enableval}" in - yes) vfc_ci=true && FCFLAGS="-D VFC_CI $FCFLAGS" && CFLAGS="-D VFC_CI $CFLAGS";; - no) vfc_ci=false ;; - *) AC_MSG_ERROR([bad value ${enableval} for --enable_vfc_ci]) ;; -esac],[vfc_ci=false]) -AM_CONDITIONAL([VFC_CI], [test x$vfc_ci = xtrue]) - -if test "$FC" = "verificarlo-f"; then - AC_MSG_NOTICE(verificarlo-f detected) - # Arguments order is important here - FCFLAGS="-Mpreprocess $FCFLAGS" -fi #PKG-CONFIG #mkl-dynamic-lp64-seq @@ -371,6 +465,9 @@ LDFLAGS:........: ${LDFLAGS} LIBS............: ${LIBS} USE CHAMELEON...: ${with_chameleon} HPC version.....: ${HAVE_HPC} +OpenMP offload..: ${HAVE_OPENMP_OFFLOAD} +OpenACC offload.: ${HAVE_OPENACC_OFFLOAD} +cuBLAS offload..: ${HAVE_CUBLAS_OFFLOAD} Package features: ${ARGS} diff --git a/org/qmckl_ao.org b/org/qmckl_ao.org index 54c5319..4120690 100644 --- a/org/qmckl_ao.org +++ b/org/qmckl_ao.org @@ -299,6 +299,8 @@ typedef struct qmckl_ao_basis_struct { uint64_t shell_vgl_date; double * restrict ao_vgl; uint64_t ao_vgl_date; + double * restrict ao_value; + uint64_t ao_value_date; int32_t uninitialized; bool provided; @@ -331,7 +333,7 @@ qmckl_exit_code qmckl_init_ao_basis(qmckl_context context) { NULL); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); ctx->ao_basis.uninitialized = (1 << 14) - 1; @@ -359,7 +361,7 @@ if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { NULL); } -qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; +qmckl_context_struct* const ctx = (qmckl_context_struct*) context; #+end_src #+NAME:post2 @@ -1376,7 +1378,7 @@ qmckl_get_ao_basis_type (const qmckl_context context, NULL); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int32_t mask = 1; @@ -1416,7 +1418,7 @@ qmckl_get_ao_basis_shell_num (const qmckl_context context, NULL); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int32_t mask = 1 << 1; @@ -1454,7 +1456,7 @@ qmckl_get_ao_basis_prim_num (const qmckl_context context, NULL); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int32_t mask = 1 << 2; @@ -1495,7 +1497,7 @@ qmckl_get_ao_basis_nucleus_shell_num (const qmckl_context context, NULL); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int32_t mask = 1 << 3; @@ -1550,7 +1552,7 @@ qmckl_get_ao_basis_nucleus_index (const qmckl_context context, NULL); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int32_t mask = 1 << 4; @@ -1606,7 +1608,7 @@ qmckl_get_ao_basis_shell_ang_mom (const qmckl_context context, NULL); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int32_t mask = 1 << 5; @@ -1662,7 +1664,7 @@ qmckl_get_ao_basis_shell_prim_num (const qmckl_context context, NULL); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int32_t mask = 1 << 6; @@ -1718,7 +1720,7 @@ qmckl_get_ao_basis_shell_prim_index (const qmckl_context context, NULL); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int32_t mask = 1 << 7; @@ -1772,7 +1774,7 @@ qmckl_get_ao_basis_shell_factor (const qmckl_context context, NULL); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int32_t mask = 1 << 8; @@ -1828,7 +1830,7 @@ qmckl_get_ao_basis_exponent (const qmckl_context context, NULL); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int32_t mask = 1 << 9; @@ -1882,7 +1884,7 @@ qmckl_get_ao_basis_coefficient (const qmckl_context context, NULL); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int32_t mask = 1 << 10; @@ -1937,7 +1939,7 @@ qmckl_get_ao_basis_prim_factor (const qmckl_context context, NULL); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int32_t mask = 1 << 11; @@ -1990,7 +1992,7 @@ qmckl_get_ao_basis_ao_num (const qmckl_context context, NULL); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int32_t mask = 1 << 12; @@ -2032,7 +2034,7 @@ qmckl_get_ao_basis_ao_factor (const qmckl_context context, NULL); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int32_t mask = 1 << 13; @@ -2079,7 +2081,7 @@ bool qmckl_ao_basis_provided(const qmckl_context context) { return false; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); return ctx->ao_basis.provided; @@ -2490,8 +2492,10 @@ free(ao_factor_test); | ~primitive_vgl_date~ | ~uint64_t~ | Last modification date of Value, gradients, Laplacian of the primitives at current positions | | ~shell_vgl~ | ~double[point_num][5][shell_num]~ | Value, gradients, Laplacian of the primitives at current positions | | ~shell_vgl_date~ | ~uint64_t~ | Last modification date of Value, gradients, Laplacian of the AOs at current positions | - | ~ao_vgl~ | ~double[point_num][5][ao_num]~ | Value, gradients, Laplacian of the primitives at current positions | + | ~ao_vgl~ | ~double[point_num][5][ao_num]~ | Value, gradients, Laplacian of the AOs at current positions | | ~ao_vgl_date~ | ~uint64_t~ | Last modification date of Value, gradients, Laplacian of the AOs at current positions | + | ~ao_value~ | ~double[point_num][ao_num]~ | Values of the the AOs at current positions | + | ~ao_value_date~ | ~uint64_t~ | Last modification date of the values of the AOs at current positions | |----------------------+-----------------------------------+----------------------------------------------------------------------------------------------| *** After initialization @@ -2518,7 +2522,7 @@ qmckl_exit_code qmckl_finalize_basis(qmckl_context context) { NULL); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int64_t nucl_num = 0; @@ -2634,9 +2638,10 @@ qmckl_exit_code qmckl_finalize_basis(qmckl_context context) { } } - rc = QMCKL_SUCCESS; #ifdef HAVE_HPC rc = qmckl_finalize_basis_hpc(context); +#else + rc = QMCKL_SUCCESS; #endif return rc; @@ -2696,7 +2701,7 @@ int compare_basis( const void * l, const void * r ) #ifdef HAVE_HPC qmckl_exit_code qmckl_finalize_basis_hpc (qmckl_context context) { - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = ctx->nucleus.num * sizeof(int32_t); @@ -2867,7 +2872,7 @@ qmckl_get_ao_basis_primitive_vgl (qmckl_context context, rc = qmckl_provide_ao_basis_primitive_vgl(context); if (rc != QMCKL_SUCCESS) return rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int64_t sze = ctx->ao_basis.prim_num * 5 * ctx->point.num; @@ -2928,7 +2933,7 @@ qmckl_get_ao_basis_shell_vgl (qmckl_context context, rc = qmckl_provide_ao_basis_shell_vgl(context); if (rc != QMCKL_SUCCESS) return rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int64_t sze = ctx->ao_basis.shell_num * 5 * ctx->point.num; @@ -2991,7 +2996,7 @@ qmckl_get_ao_basis_ao_vgl (qmckl_context context, rc = qmckl_provide_ao_vgl(context); if (rc != QMCKL_SUCCESS) return rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int64_t sze = ctx->ao_basis.ao_num * 5 * ctx->point.num; @@ -3021,7 +3026,7 @@ qmckl_get_ao_basis_ao_vgl (qmckl_context context, end interface #+end_src - Uses the give array to compute the VGL. + Uses the given array to compute the VGL. #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code @@ -3029,7 +3034,7 @@ qmckl_get_ao_basis_ao_vgl_inplace (qmckl_context context, double* const ao_vgl, const int64_t size_max); #+end_src - + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_ao_basis_ao_vgl_inplace (qmckl_context context, @@ -3046,7 +3051,7 @@ qmckl_get_ao_basis_ao_vgl_inplace (qmckl_context context, qmckl_exit_code rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int64_t sze = ctx->ao_basis.ao_num * 5 * ctx->point.num; @@ -3087,6 +3092,133 @@ qmckl_get_ao_basis_ao_vgl_inplace (qmckl_context context, end interface #+end_src + + + #+begin_src c :comments org :tangle (eval h_func) :noweb yes +qmckl_exit_code +qmckl_get_ao_basis_ao_value (qmckl_context context, + double* const ao_value, + const int64_t size_max); + #+end_src + + Returns the array of values of the atomic orbitals evaluated at + the current coordinates. See section [[Combining radial and polynomial parts]]. + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code +qmckl_get_ao_basis_ao_value (qmckl_context context, + double* const ao_value, + const int64_t size_max) +{ + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return qmckl_failwith( context, + QMCKL_INVALID_CONTEXT, + "qmckl_get_ao_basis_ao_value", + NULL); + } + + qmckl_exit_code rc; + + rc = qmckl_provide_ao_value(context); + if (rc != QMCKL_SUCCESS) return rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; + assert (ctx != NULL); + + int64_t sze = ctx->ao_basis.ao_num * ctx->point.num; + if (size_max < sze) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_get_ao_basis_ao_value", + "input array too small"); + } + memcpy(ao_value, ctx->ao_basis.ao_value, (size_t) sze * sizeof(double)); + + return QMCKL_SUCCESS; +} + #+end_src + + #+begin_src f90 :tangle (eval fh_func) :comments org :exports none + interface + integer(c_int32_t) function qmckl_get_ao_basis_ao_value (context, & + ao_value, size_max) bind(C) + use, intrinsic :: iso_c_binding + import + implicit none + integer (c_int64_t) , intent(in) , value :: context + double precision, intent(out) :: ao_value(*) + integer (c_int64_t) , intent(in) , value :: size_max + end function qmckl_get_ao_basis_ao_value + end interface + #+end_src + + Uses the given array to compute the value. + + #+begin_src c :comments org :tangle (eval h_func) :noweb yes +qmckl_exit_code +qmckl_get_ao_basis_ao_value_inplace (qmckl_context context, + double* const ao_value, + const int64_t size_max); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code +qmckl_get_ao_basis_ao_value_inplace (qmckl_context context, + double* const ao_value, + const int64_t size_max) +{ + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return qmckl_failwith( context, + QMCKL_INVALID_CONTEXT, + "qmckl_get_ao_basis_ao_value", + NULL); + } + + qmckl_exit_code rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; + assert (ctx != NULL); + + int64_t sze = ctx->ao_basis.ao_num * ctx->point.num; + if (size_max < sze) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_get_ao_basis_ao_value", + "input array too small"); + } + + rc = qmckl_context_touch(context); + if (rc != QMCKL_SUCCESS) return rc; + + double* old_array = ctx->ao_basis.ao_value; + + ctx->ao_basis.ao_value = ao_value; + + rc = qmckl_provide_ao_value(context); + if (rc != QMCKL_SUCCESS) return rc; + + ctx->ao_basis.ao_value = old_array; + + return QMCKL_SUCCESS; +} + #+end_src + + #+begin_src f90 :tangle (eval fh_func) :comments org :exports none + interface + integer(c_int32_t) function qmckl_get_ao_basis_ao_value_inplace (context, & + ao_value, size_max) bind(C) + use, intrinsic :: iso_c_binding + import + implicit none + integer (c_int64_t) , intent(in) , value :: context + double precision, intent(out) :: ao_value(*) + integer (c_int64_t) , intent(in) , value :: size_max + end function qmckl_get_ao_basis_ao_value_inplace + end interface + #+end_src + * Radial part ** General functions for Gaussian basis functions @@ -3487,7 +3619,7 @@ qmckl_exit_code qmckl_provide_ao_basis_primitive_vgl(qmckl_context context) NULL); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); if (!ctx->ao_basis.provided) { @@ -3888,7 +4020,7 @@ qmckl_exit_code qmckl_provide_ao_basis_shell_vgl(qmckl_context context) NULL); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); if (!ctx->ao_basis.provided) { @@ -4796,13 +4928,13 @@ qmckl_ao_polynomial_transp_vgl_hpc (const qmckl_context context, const double* restrict X, const double* restrict R, const int32_t lmax, - int64_t* restrict n, + int64_t* n, int32_t* restrict const L, const int64_t ldl, double* restrict const VGL, const int64_t ldv ) { - const qmckl_context_struct* ctx = (qmckl_context_struct* const) context; + const qmckl_context_struct* ctx = (qmckl_context_struct*) context; assert (ctx != NULL && X != NULL && R != NULL && n != NULL && L != NULL && VGL != NULL); if (lmax < 0) return QMCKL_INVALID_ARG_4; if (ldl < 3) return QMCKL_INVALID_ARG_7; @@ -5284,33 +5416,780 @@ for (int32_t ldl=3 ; ldl<=5 ; ++ldl) { #+end_src * Combining radial and polynomial parts + +** Values only + :PROPERTIES: + :Name: qmckl_compute_ao_value + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: +*** Unoptimized version + #+NAME: qmckl_ao_value_args_doc + | Variable | Type | In/Out | Description | + |-----------------------+-----------------------------------+--------+----------------------------------------------| + | ~context~ | ~qmckl_context~ | in | Global state | + | ~ao_num~ | ~int64_t~ | in | Number of AOs | + | ~shell_num~ | ~int64_t~ | in | Number of shells | + | ~point_num~ | ~int64_t~ | in | Number of points | + | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | + | ~coord~ | ~double[3][point_num]~ | in | Coordinates | + | ~nucl_coord~ | ~double[3][nucl_num]~ | in | Nuclear coordinates | + | ~nucleus_index~ | ~int64_t[nucl_num]~ | in | Index of the 1st shell of each nucleus | + | ~nucleus_shell_num~ | ~int64_t[nucl_num]~ | in | Number of shells per nucleus | + | ~nucleus_range~ | ~double[nucl_num]~ | in | Range beyond which all is zero | + | ~nucleus_max_ang_mom~ | ~int32_t[nucl_num]~ | in | Maximum angular momentum per nucleus | + | ~shell_ang_mom~ | ~int32_t[shell_num]~ | in | Angular momentum of each shell | + | ~ao_factor~ | ~double[ao_num]~ | in | Normalization factor of the AOs | + | ~shell_vgl~ | ~double[point_num][5][shell_num]~ | in | Value, gradients and Laplacian of the shells | + | ~ao_value~ | ~double[point_num][ao_num]~ | out | Values of the AOs | + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer function qmckl_compute_ao_value_doc_f(context, & + ao_num, shell_num, point_num, nucl_num, & + coord, nucl_coord, nucleus_index, nucleus_shell_num, & + nucleus_range, nucleus_max_ang_mom, shell_ang_mom, & + ao_factor, shell_vgl, ao_value) & + result(info) + use qmckl + implicit none + integer(qmckl_context), intent(in) :: context + integer*8 , intent(in) :: ao_num + integer*8 , intent(in) :: shell_num + integer*8 , intent(in) :: point_num + integer*8 , intent(in) :: nucl_num + double precision , intent(in) :: coord(point_num,3) + double precision , intent(in) :: nucl_coord(nucl_num,3) + integer*8 , intent(in) :: nucleus_index(nucl_num) + integer*8 , intent(in) :: nucleus_shell_num(nucl_num) + double precision , intent(in) :: nucleus_range(nucl_num) + integer , intent(in) :: nucleus_max_ang_mom(nucl_num) + integer , intent(in) :: shell_ang_mom(shell_num) + double precision , intent(in) :: ao_factor(ao_num) + double precision , intent(in) :: shell_vgl(shell_num,5,point_num) + double precision , intent(out) :: ao_value(ao_num,point_num) + + double precision :: e_coord(3), n_coord(3) + integer*8 :: n_poly + integer :: l, il, k + integer*8 :: ipoint, inucl, ishell + integer*8 :: ishell_start, ishell_end + integer :: lstart(0:20) + double precision :: x, y, z, r2 + double precision :: cutoff + integer, external :: qmckl_ao_polynomial_vgl_doc_f + + double precision, allocatable :: poly_vgl(:,:) + integer , allocatable :: powers(:,:), ao_index(:) + + allocate(poly_vgl(5,ao_num), powers(3,ao_num), ao_index(ao_num)) + + ! Pre-computed data + do l=0,20 + lstart(l) = l*(l+1)*(l+2)/6 +1 + end do + + k=1 + do inucl=1,nucl_num + ishell_start = nucleus_index(inucl) + 1 + ishell_end = nucleus_index(inucl) + nucleus_shell_num(inucl) + do ishell = ishell_start, ishell_end + l = shell_ang_mom(ishell) + ao_index(ishell) = k + k = k + lstart(l+1) - lstart(l) + end do + end do + info = QMCKL_SUCCESS + + ! Don't compute polynomials when the radial part is zero. + cutoff = -dlog(1.d-12) + + do ipoint = 1, point_num + e_coord(1) = coord(ipoint,1) + e_coord(2) = coord(ipoint,2) + e_coord(3) = coord(ipoint,3) + do inucl=1,nucl_num + n_coord(1) = nucl_coord(inucl,1) + n_coord(2) = nucl_coord(inucl,2) + n_coord(3) = nucl_coord(inucl,3) + + ! Test if the point is in the range of the nucleus + x = e_coord(1) - n_coord(1) + y = e_coord(2) - n_coord(2) + z = e_coord(3) - n_coord(3) + + r2 = x*x + y*y + z*z + + if (r2 > cutoff*nucleus_range(inucl)) then + cycle + end if + + ! Compute polynomials + info = qmckl_ao_polynomial_vgl_doc_f(context, e_coord, n_coord, & + nucleus_max_ang_mom(inucl), n_poly, powers, 3_8, & + poly_vgl, 5_8) + + ! Loop over shells + ishell_start = nucleus_index(inucl) + 1 + ishell_end = nucleus_index(inucl) + nucleus_shell_num(inucl) + do ishell = ishell_start, ishell_end + k = ao_index(ishell) + l = shell_ang_mom(ishell) + do il = lstart(l), lstart(l+1)-1 + ! Value + ao_value(k,ipoint) = & + poly_vgl(1,il) * shell_vgl(ishell,1,ipoint) * ao_factor(k) + k = k+1 + end do + end do + end do + end do + + deallocate(poly_vgl, powers) +end function qmckl_compute_ao_value_doc_f + #+end_src + +*** HPC version + #+NAME: qmckl_ao_value_args_hpc_gaussian + | Variable | Type | In/Out | Description | + |-----------------------+-----------------------------+--------+----------------------------------------------| + | ~context~ | ~qmckl_context~ | in | Global state | + | ~ao_num~ | ~int64_t~ | in | Number of AOs | + | ~shell_num~ | ~int64_t~ | in | Number of shells | + | ~prim_num~ | ~int64_t~ | in | Number of primitives | + | ~point_num~ | ~int64_t~ | in | Number of points | + | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | + | ~coord~ | ~double[3][point_num]~ | in | Coordinates | + | ~nucl_coord~ | ~double[3][nucl_num]~ | in | Nuclear coordinates | + | ~nucleus_index~ | ~int64_t[nucl_num]~ | in | Index of the 1st shell of each nucleus | + | ~nucleus_shell_num~ | ~int64_t[nucl_num]~ | in | Number of shells per nucleus | + | ~nucleus_range~ | ~double[nucl_num]~ | in | Range beyond which all is zero | + | ~nucleus_max_ang_mom~ | ~int32_t[nucl_num]~ | in | Maximum angular momentum per nucleus | + | ~shell_ang_mom~ | ~int32_t[shell_num]~ | in | Angular momentum of each shell | + | ~shell_prim_index~ | ~int64_t[shell_num]~ | in | Index of the 1st primitive of each shell | + | ~shell_prim_num~ | ~int64_t[shell_num]~ | in | Number of primitives per shell | + | ~ao_factor~ | ~double[ao_num]~ | in | Normalization factor of the AOs | + | ~ao_expo~ | ~double[prim_num]~ | in | Value, gradients and Laplacian of the shells | + | ~coef_normalized~ | ~double[prim_num]~ | in | Value, gradients and Laplacian of the shells | + | ~ao_value~ | ~double[point_num][ao_num]~ | out | Values of the AOs | + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +#ifdef HAVE_HPC +qmckl_exit_code +qmckl_compute_ao_value_hpc_gaussian (const qmckl_context context, + const int64_t ao_num, + const int64_t shell_num, + const int32_t* restrict prim_num_per_nucleus, + const int64_t point_num, + const int64_t nucl_num, + const double* restrict coord, + const double* restrict nucl_coord, + const int64_t* restrict nucleus_index, + const int64_t* restrict nucleus_shell_num, + const double* nucleus_range, + const int32_t* restrict nucleus_max_ang_mom, + const int32_t* restrict shell_ang_mom, + const double* restrict ao_factor, + const qmckl_matrix expo_per_nucleus, + const qmckl_tensor coef_per_nucleus, + double* restrict const ao_value ) +{ + int32_t lstart[32]; + for (int32_t l=0 ; l<32 ; ++l) { + lstart[l] = l*(l+1)*(l+2)/6; + } + + int64_t ao_index[shell_num+1]; + int64_t size_max = 0; + int64_t prim_max = 0; + int64_t shell_max = 0; + { + int64_t k=0; + for (int inucl=0 ; inucl < nucl_num ; ++inucl) { + prim_max = prim_num_per_nucleus[inucl] > prim_max ? + prim_num_per_nucleus[inucl] : prim_max; + shell_max = nucleus_shell_num[inucl] > shell_max ? + nucleus_shell_num[inucl] : shell_max; + const int64_t ishell_start = nucleus_index[inucl]; + const int64_t ishell_end = nucleus_index[inucl] + nucleus_shell_num[inucl]; + for (int64_t ishell = ishell_start ; ishell < ishell_end ; ++ishell) { + const int l = shell_ang_mom[ishell]; + ao_index[ishell] = k; + k += lstart[l+1] - lstart[l]; + size_max = size_max < lstart[l+1] ? lstart[l+1] : size_max; + } + } + ao_index[shell_num] = ao_num+1; + } + + /* Don't compute polynomials when the radial part is zero. */ + double cutoff = -log(1.e-12); + +#ifdef HAVE_OPENMP +#pragma omp parallel +#endif + { + qmckl_exit_code rc; + double ar2[prim_max]; + int32_t powers[3*size_max]; + double poly_vgl[5*size_max]; + + double exp_mat[prim_max]; + double ce_mat[shell_max]; + + double coef_mat[nucl_num][shell_max][prim_max]; + for (int i=0 ; i cutoff * nucleus_range[inucl]) { + continue; + } + + int64_t n_poly; + switch (nucleus_max_ang_mom[inucl]) { + case 0: + break; + + case 1: + poly_vgl[0] = 0.; + poly_vgl[1] = x; + poly_vgl[2] = y; + poly_vgl[3] = z; + break; + + case 2: + poly_vgl[0] = 0.; + poly_vgl[1] = x; + poly_vgl[2] = y; + poly_vgl[3] = z; + poly_vgl[4] = x*x; + poly_vgl[5] = x*y; + poly_vgl[6] = x*z; + poly_vgl[7] = y*y; + poly_vgl[8] = y*z; + poly_vgl[9] = z*z; + break; + + default: + rc = qmckl_ao_polynomial_transp_vgl_hpc(context, e_coord, n_coord, + nucleus_max_ang_mom[inucl], + &n_poly, powers, (int64_t) 3, + poly_vgl, size_max); + assert (rc == QMCKL_SUCCESS); + break; + } + + /* Compute all exponents */ + + int64_t nidx = 0; + for (int64_t iprim = 0 ; iprim < prim_num_per_nucleus[inucl] ; ++iprim) { + const double v = qmckl_mat(expo_per_nucleus, iprim, inucl) * r2; + if (v <= cutoff) { + ar2[iprim] = v; + ++nidx; + } else { + break; + } + } + + for (int64_t iprim = 0 ; iprim < nidx ; ++iprim) { + exp_mat[iprim] = exp(-ar2[iprim]); + } + + for (int i=0 ; i 0) { + const double* restrict f = ao_factor + k; + const int64_t idx = lstart[l]; + + poly_vgl_1 = &(poly_vgl[idx]); + + switch (n) { + case(1): + ao_value_1[0] = s1 * f[0]; + break; + case (3): +#ifdef HAVE_OPENMP +#pragma omp simd +#endif + for (int il=0 ; il<3 ; ++il) { + ao_value_1[il] = poly_vgl_1[il] * s1 * f[il]; + } + break; + case(6): +#ifdef HAVE_OPENMP +#pragma omp simd +#endif + for (int il=0 ; il<6 ; ++il) { + ao_value_1[il] = poly_vgl_1[il] * s1 * f[il]; + } + break; + default: +#ifdef HAVE_OPENMP +#pragma omp simd simdlen(8) +#endif + for (int il=0 ; ilao_basis.provided) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_ao_value", + NULL); + } + + /* Compute if necessary */ + if (ctx->point.date > ctx->ao_basis.ao_value_date) { + + qmckl_exit_code rc; + + /* Provide required data */ +#ifndef HAVE_HPC + rc = qmckl_provide_ao_basis_shell_vgl(context); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, "qmckl_provide_ao_basis_shell_vgl", NULL); + } +#endif + + /* Allocate array */ + if (ctx->ao_basis.ao_value == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->ao_basis.ao_num * 5 * ctx->point.num * sizeof(double); + double* ao_value = (double*) qmckl_malloc(context, mem_info); + + if (ao_value == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_ao_basis_ao_value", + NULL); + } + ctx->ao_basis.ao_value = ao_value; + } + + if (ctx->ao_basis.ao_vgl_date == ctx->point.date) { + + // ao_vgl has been computed at this step: Just copy the data. + + double * v = &(ctx->ao_basis.ao_value[0]); + double * vgl = &(ctx->ao_basis.ao_vgl[0]); + for (int i=0 ; ipoint.num ; ++i) { + for (int k=0 ; kao_basis.ao_num ; ++k) { + v[k] = vgl[k]; + } + v += ctx->ao_basis.ao_num; + vgl += ctx->ao_basis.ao_num * 5; + } + + } else { + +#ifdef HAVE_HPC + if (ctx->ao_basis.type == 'G') { + rc = qmckl_compute_ao_value_hpc_gaussian(context, + ctx->ao_basis.ao_num, + ctx->ao_basis.shell_num, + ctx->ao_basis.prim_num_per_nucleus, + ctx->point.num, + ctx->nucleus.num, + ctx->point.coord.data, + ctx->nucleus.coord.data, + ctx->ao_basis.nucleus_index, + ctx->ao_basis.nucleus_shell_num, + ctx->ao_basis.nucleus_range, + ctx->ao_basis.nucleus_max_ang_mom, + ctx->ao_basis.shell_ang_mom, + ctx->ao_basis.ao_factor, + ctx->ao_basis.expo_per_nucleus, + ctx->ao_basis.coef_per_nucleus, + ctx->ao_basis.ao_value); + /* + } else if (ctx->ao_basis.type == 'S') { + rc = qmck_compute_ao_value_hpc_slater(context, + ctx->ao_basis.ao_num, + ctx->ao_basis.shell_num, + ctx->ao_basis.prim_num, + ctx->point.num, + ctx->nucleus.num, + ctx->point.coord.data, + ctx->nucleus.coord.data, + ctx->ao_basis.nucleus_index, + ctx->ao_basis.nucleus_shell_num, + ctx->ao_basis.nucleus_range, + ctx->ao_basis.nucleus_max_ang_mom, + ctx->ao_basis.shell_ang_mom, + ctx->ao_basis.shell_prim_index, + ctx->ao_basis.shell_prim_num, + ctx->ao_basis.ao_factor, + ctx->ao_basis.exponent, + ctx->ao_basis.coefficient_normalized, + ctx->ao_basis.ao_value); + ,*/ + } else { + rc = qmckl_compute_ao_value_doc(context, + ctx->ao_basis.ao_num, + ctx->ao_basis.shell_num, + ctx->point.num, + ctx->nucleus.num, + ctx->point.coord.data, + ctx->nucleus.coord.data, + ctx->ao_basis.nucleus_index, + ctx->ao_basis.nucleus_shell_num, + ctx->ao_basis.nucleus_range, + ctx->ao_basis.nucleus_max_ang_mom, + ctx->ao_basis.shell_ang_mom, + ctx->ao_basis.ao_factor, + ctx->ao_basis.shell_vgl, + ctx->ao_basis.ao_value); + } +#else + rc = qmckl_compute_ao_value_doc(context, + ctx->ao_basis.ao_num, + ctx->ao_basis.shell_num, + ctx->point.num, + ctx->nucleus.num, + ctx->point.coord.data, + ctx->nucleus.coord.data, + ctx->ao_basis.nucleus_index, + ctx->ao_basis.nucleus_shell_num, + ctx->ao_basis.nucleus_range, + ctx->ao_basis.nucleus_max_ang_mom, + ctx->ao_basis.shell_ang_mom, + ctx->ao_basis.ao_factor, + ctx->ao_basis.shell_vgl, + ctx->ao_basis.ao_value); +#endif + if (rc != QMCKL_SUCCESS) { + return rc; + } + + } + + ctx->ao_basis.ao_value_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + #+end_src + +**** Test :noexport: + + #+begin_src python :results output :exports none +import numpy as np +from math import sqrt + +h0 = 1.e-4 +def f(a,x,y): + return np.sum( [c * np.exp( -b*(np.linalg.norm(x-y))**2) for b,c in a] ) + +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( [ -2.302574592081335e+00, -3.542027060505035e-01, -5.334129934317614e-02] ) + +#double ao_value[prim_num][5][elec_num]; +x = elec_26_w1 ; y = nucl_1 +a = [( 4.0382999999999998e+02, 1.4732000000000000e-03 * 5.9876577632594533e+04), + ( 1.2117000000000000e+02, 1.2672500000000000e-02 * 7.2836806319891484e+03), + ( 4.6344999999999999e+01, 5.8045100000000002e-02 * 1.3549226646722386e+03), + ( 1.9721000000000000e+01, 1.7051030000000000e-01 * 3.0376315094739988e+02), + ( 8.8623999999999992e+00, 3.1859579999999998e-01 * 7.4924579607137730e+01), + ( 3.9962000000000000e+00, 3.8450230000000002e-01 * 1.8590543353806009e+01), + ( 1.7636000000000001e+00, 2.7377370000000001e-01 * 4.4423176930919421e+00), + ( 7.0618999999999998e-01, 7.4396699999999996e-02 * 8.9541051939952665e-01)] + +norm = sqrt(3.) +# x^2 * g(r) +print ( "[26][0][219] : %25.15e"%(fx(a,x,y)) ) +print ( "[26][1][219] : %25.15e"%(df(a,x,y,1)) ) +print ( "[26][2][219] : %25.15e"%(df(a,x,y,2)) ) +print ( "[26][3][219] : %25.15e"%(df(a,x,y,3)) ) +print ( "[26][4][219] : %25.15e"%(lf(a,x,y)) ) + +print ( "[26][0][220] : %25.15e"%(norm*f(a,x,y) * (x[0] - y[0]) * (x[1] - y[1]) )) +print ( "[26][1][220] : %25.15e"%(norm*df(a,x,y,1)* (x[0] - y[0]) * (x[1] - y[1]) + norm*f(a,x,y) * (x[1] - y[1])) ) + +print ( "[26][0][221] : %25.15e"%(norm*f(a,x,y) * (x[0] - y[0]) * (x[2] - y[2])) ) +print ( "[26][1][221] : %25.15e"%(norm*df(a,x,y,1)* (x[0] - y[0]) * (x[2] - y[2]) + norm*f(a,x,y) * (x[2] - y[2])) ) + +print ( "[26][0][222] : %25.15e"%(f(a,x,y) * (x[1] - y[1]) * (x[1] - y[1])) ) +print ( "[26][1][222] : %25.15e"%(df(a,x,y,1)* (x[1] - y[1]) * (x[1] - y[1])) ) + +print ( "[26][0][223] : %25.15e"%(norm*f(a,x,y) * (x[1] - y[1]) * (x[2] - y[2])) ) +print ( "[26][1][223] : %25.15e"%(norm*df(a,x,y,1)* (x[1] - y[1]) * (x[2] - y[2])) ) + +print ( "[26][0][224] : %25.15e"%(f(a,x,y) * (x[2] - y[2]) * (x[2] - y[2])) ) +print ( "[26][1][224] : %25.15e"%(df(a,x,y,1)* (x[2] - y[2]) * (x[2] - y[2])) ) + + #+end_src + + #+RESULTS: + + #+begin_src c :tangle (eval c_test) :exports none +{ +#define walk_num 1 // chbrclf_walk_num +#define elec_num chbrclf_elec_num +#define shell_num chbrclf_shell_num +#define ao_num chbrclf_ao_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, walk_num*elec_num*3); +assert(rc == QMCKL_SUCCESS); + + +double ao_value[elec_num][ao_num]; + +rc = qmckl_get_ao_basis_ao_value(context, &(ao_value[0][0]), + (int64_t) elec_num*ao_num); +assert (rc == QMCKL_SUCCESS); + +printf("\n"); +printf(" ao_value ao_value[26][219] %25.15e\n", ao_value[26][219]); +printf(" ao_value ao_value[26][220] %25.15e\n", ao_value[26][220]); +printf(" ao_value ao_value[26][221] %25.15e\n", ao_value[26][221]); +printf(" ao_value ao_value[26][222] %25.15e\n", ao_value[26][222]); +printf(" ao_value ao_value[26][223] %25.15e\n", ao_value[26][223]); +printf(" ao_value ao_value[26][224] %25.15e\n", ao_value[26][224]); +printf("\n"); + +assert( fabs(ao_value[26][219] - ( 1.020298798341620e-08)) < 1.e-14 ); +assert( fabs(ao_value[26][220] - ( 1.516643537739178e-08)) < 1.e-14 ); +assert( fabs(ao_value[26][221] - ( -4.686370882518819e-09)) < 1.e-14 ); +assert( fabs(ao_value[26][222] - ( 7.514816980753531e-09)) < 1.e-14 ); +assert( fabs(ao_value[26][223] - ( -4.021908374204471e-09)) < 1.e-14 ); +assert( fabs(ao_value[26][224] - ( 7.175045873560788e-10)) < 1.e-14 ); + +} + + #+end_src + +** Value, gradients, Laplacian :PROPERTIES: :Name: qmckl_compute_ao_vgl :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: +*** Unoptimized version + #+NAME: qmckl_ao_vgl_args_doc + | Variable | Type | In/Out | Description | + |-----------------------+-----------------------------------+--------+----------------------------------------------| + | ~context~ | ~qmckl_context~ | in | Global state | + | ~ao_num~ | ~int64_t~ | in | Number of AOs | + | ~shell_num~ | ~int64_t~ | in | Number of shells | + | ~point_num~ | ~int64_t~ | in | Number of points | + | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | + | ~coord~ | ~double[3][point_num]~ | in | Coordinates | + | ~nucl_coord~ | ~double[3][nucl_num]~ | in | Nuclear coordinates | + | ~nucleus_index~ | ~int64_t[nucl_num]~ | in | Index of the 1st shell of each nucleus | + | ~nucleus_shell_num~ | ~int64_t[nucl_num]~ | in | Number of shells per nucleus | + | ~nucleus_range~ | ~double[nucl_num]~ | in | Range beyond which all is zero | + | ~nucleus_max_ang_mom~ | ~int32_t[nucl_num]~ | in | Maximum angular momentum per nucleus | + | ~shell_ang_mom~ | ~int32_t[shell_num]~ | in | Angular momentum of each shell | + | ~ao_factor~ | ~double[ao_num]~ | in | Normalization factor of the AOs | + | ~shell_vgl~ | ~double[point_num][5][shell_num]~ | in | Value, gradients and Laplacian of the shells | + | ~ao_vgl~ | ~double[point_num][5][ao_num]~ | out | Value, gradients and Laplacian of the AOs | -** Unoptimized version - #+NAME: qmckl_ao_vgl_args_doc - | Variable | Type | In/Out | Description | - |-----------------------+-----------------------------------+--------+----------------------------------------------| - | ~context~ | ~qmckl_context~ | in | Global state | - | ~ao_num~ | ~int64_t~ | in | Number of AOs | - | ~shell_num~ | ~int64_t~ | in | Number of shells | - | ~point_num~ | ~int64_t~ | in | Number of points | - | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | - | ~coord~ | ~double[3][point_num]~ | in | Coordinates | - | ~nucl_coord~ | ~double[3][nucl_num]~ | in | Nuclear coordinates | - | ~nucleus_index~ | ~int64_t[nucl_num]~ | in | Index of the 1st shell of each nucleus | - | ~nucleus_shell_num~ | ~int64_t[nucl_num]~ | in | Number of shells per nucleus | - | ~nucleus_range~ | ~double[nucl_num]~ | in | Range beyond which all is zero | - | ~nucleus_max_ang_mom~ | ~int32_t[nucl_num]~ | in | Maximum angular momentum per nucleus | - | ~shell_ang_mom~ | ~int32_t[shell_num]~ | in | Angular momentum of each shell | - | ~ao_factor~ | ~double[ao_num]~ | in | Normalization factor of the AOs | - | ~shell_vgl~ | ~double[point_num][5][shell_num]~ | in | Value, gradients and Laplacian of the shells | - | ~ao_vgl~ | ~double[point_num][5][ao_num]~ | out | Value, gradients and Laplacian of the AOs | - - #+begin_src f90 :comments org :tangle (eval f) :noweb yes + #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_ao_vgl_doc_f(context, & ao_num, shell_num, point_num, nucl_num, & coord, nucl_coord, nucleus_index, nucleus_shell_num, & @@ -5442,34 +6321,34 @@ integer function qmckl_compute_ao_vgl_doc_f(context, & deallocate(poly_vgl, powers) end function qmckl_compute_ao_vgl_doc_f - #+end_src + #+end_src -** HPC version - #+NAME: qmckl_ao_vgl_args_hpc_gaussian - | Variable | Type | In/Out | Description | - |-----------------------+--------------------------------+--------+----------------------------------------------| - | ~context~ | ~qmckl_context~ | in | Global state | - | ~ao_num~ | ~int64_t~ | in | Number of AOs | - | ~shell_num~ | ~int64_t~ | in | Number of shells | - | ~prim_num~ | ~int64_t~ | in | Number of primitives | - | ~point_num~ | ~int64_t~ | in | Number of points | - | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | - | ~coord~ | ~double[3][point_num]~ | in | Coordinates | - | ~nucl_coord~ | ~double[3][nucl_num]~ | in | Nuclear coordinates | - | ~nucleus_index~ | ~int64_t[nucl_num]~ | in | Index of the 1st shell of each nucleus | - | ~nucleus_shell_num~ | ~int64_t[nucl_num]~ | in | Number of shells per nucleus | - | ~nucleus_range~ | ~double[nucl_num]~ | in | Range beyond which all is zero | - | ~nucleus_max_ang_mom~ | ~int32_t[nucl_num]~ | in | Maximum angular momentum per nucleus | - | ~shell_ang_mom~ | ~int32_t[shell_num]~ | in | Angular momentum of each shell | - | ~shell_prim_index~ | ~int64_t[shell_num]~ | in | Index of the 1st primitive of each shell | - | ~shell_prim_num~ | ~int64_t[shell_num]~ | in | Number of primitives per shell | - | ~ao_factor~ | ~double[ao_num]~ | in | Normalization factor of the AOs | - | ~ao_expo~ | ~double[prim_num]~ | in | Value, gradients and Laplacian of the shells | - | ~coef_normalized~ | ~double[prim_num]~ | in | Value, gradients and Laplacian of the shells | - | ~ao_vgl~ | ~double[point_num][5][ao_num]~ | out | Value, gradients and Laplacian of the AOs | +*** HPC version + #+NAME: qmckl_ao_vgl_args_hpc_gaussian + | Variable | Type | In/Out | Description | + |-----------------------+--------------------------------+--------+----------------------------------------------| + | ~context~ | ~qmckl_context~ | in | Global state | + | ~ao_num~ | ~int64_t~ | in | Number of AOs | + | ~shell_num~ | ~int64_t~ | in | Number of shells | + | ~prim_num~ | ~int64_t~ | in | Number of primitives | + | ~point_num~ | ~int64_t~ | in | Number of points | + | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | + | ~coord~ | ~double[3][point_num]~ | in | Coordinates | + | ~nucl_coord~ | ~double[3][nucl_num]~ | in | Nuclear coordinates | + | ~nucleus_index~ | ~int64_t[nucl_num]~ | in | Index of the 1st shell of each nucleus | + | ~nucleus_shell_num~ | ~int64_t[nucl_num]~ | in | Number of shells per nucleus | + | ~nucleus_range~ | ~double[nucl_num]~ | in | Range beyond which all is zero | + | ~nucleus_max_ang_mom~ | ~int32_t[nucl_num]~ | in | Maximum angular momentum per nucleus | + | ~shell_ang_mom~ | ~int32_t[shell_num]~ | in | Angular momentum of each shell | + | ~shell_prim_index~ | ~int64_t[shell_num]~ | in | Index of the 1st primitive of each shell | + | ~shell_prim_num~ | ~int64_t[shell_num]~ | in | Number of primitives per shell | + | ~ao_factor~ | ~double[ao_num]~ | in | Normalization factor of the AOs | + | ~ao_expo~ | ~double[prim_num]~ | in | Value, gradients and Laplacian of the shells | + | ~coef_normalized~ | ~double[prim_num]~ | in | Value, gradients and Laplacian of the shells | + | ~ao_vgl~ | ~double[point_num][5][ao_num]~ | out | Value, gradients and Laplacian of the AOs | - #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none #ifdef HAVE_HPC qmckl_exit_code qmckl_compute_ao_vgl_hpc_gaussian ( @@ -5528,7 +6407,7 @@ qmckl_compute_ao_vgl_hpc_gaussian ( { qmckl_exit_code rc; double ar2[prim_max]; - int32_t powers[prim_max]; + int32_t powers[3*size_max]; double poly_vgl_l1[4][4] = {{1.0, 0.0, 0.0, 0.0}, {0.0, 1.0, 0.0, 0.0}, {0.0, 0.0, 1.0, 0.0}, @@ -5613,7 +6492,6 @@ qmckl_compute_ao_vgl_hpc_gaussian ( nucleus_max_ang_mom[inucl], &n_poly, powers, (int64_t) 3, &(poly_vgl[0][0]), size_max); - assert (rc == QMCKL_SUCCESS); break; } @@ -5784,14 +6662,13 @@ qmckl_compute_ao_vgl_hpc_gaussian ( return QMCKL_SUCCESS; } #endif - #+end_src + #+end_src -** Interfaces -# #+CALL: generate_c_header(table=qmckl_ao_vgl_args_doc,rettyp=get_value("CRetType"),fname="qmckl_compute_ao_vgl")) -# (Commented because the header needs to go into h_private_func) +*** Interfaces + # #+CALL: generate_c_header(table=qmckl_ao_vgl_args_doc,rettyp=get_value("CRetType"),fname="qmckl_compute_ao_vgl")) + # (Commented because the header needs to go into h_private_func) - #+RESULTS: - #+begin_src c :tangle (eval h_private_func) :comments org + #+begin_src c :tangle (eval h_private_func) :comments org qmckl_exit_code qmckl_compute_ao_vgl_doc ( const qmckl_context context, const int64_t ao_num, @@ -5808,8 +6685,9 @@ qmckl_compute_ao_vgl_hpc_gaussian ( const double* ao_factor, const double* shell_vgl, double* const ao_vgl ); - #+end_src - #+begin_src c :tangle (eval h_private_func) :comments org + #+end_src + + #+begin_src c :tangle (eval h_private_func) :comments org #ifdef HAVE_HPC qmckl_exit_code qmckl_compute_ao_vgl_hpc_gaussian ( const qmckl_context context, @@ -5830,12 +6708,12 @@ qmckl_compute_ao_vgl_hpc_gaussian ( const qmckl_tensor coef_per_nucleus, double* const ao_vgl ); #endif - #+end_src + #+end_src - #+CALL: generate_c_interface(table=qmckl_ao_vgl_args_doc,rettyp=get_value("CRetType"),fname="qmckl_compute_ao_vgl_doc")) + #+CALL: generate_c_interface(table=qmckl_ao_vgl_args_doc,rettyp=get_value("CRetType"),fname="qmckl_compute_ao_vgl_doc")) - #+RESULTS: - #+begin_src f90 :tangle (eval f) :comments org :exports none + #+RESULTS: + #+begin_src f90 :tangle (eval f) :comments org :exports none integer(c_int32_t) function qmckl_compute_ao_vgl_doc & (context, & ao_num, & @@ -5892,16 +6770,15 @@ qmckl_compute_ao_vgl_hpc_gaussian ( ao_vgl) end function qmckl_compute_ao_vgl_doc - #+end_src + #+end_src +**** Provide :noexport: -*** Provide :noexport: - - #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context); - #+end_src + #+end_src - #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context) { @@ -5912,7 +6789,7 @@ qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context) NULL); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); if (!ctx->ao_basis.provided) { @@ -6034,11 +6911,11 @@ qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context) return QMCKL_SUCCESS; } - #+end_src + #+end_src -*** Test :noexport: +**** Test :noexport: - #+begin_src python :results output :exports none + #+begin_src python :results output :exports none import numpy as np from math import sqrt @@ -6107,11 +6984,11 @@ print ( "[26][1][223] : %25.15e"%(norm*df(a,x,y,1)* (x[1] - y[1]) * (x[2] - y[2] print ( "[26][0][224] : %25.15e"%(f(a,x,y) * (x[2] - y[2]) * (x[2] - y[2])) ) print ( "[26][1][224] : %25.15e"%(df(a,x,y,1)* (x[2] - y[2]) * (x[2] - y[2])) ) - #+end_src + #+end_src - #+RESULTS: + #+RESULTS: - #+begin_src c :tangle (eval c_test) :exports none + #+begin_src c :tangle (eval c_test) :exports none { #define walk_num 1 // chbrclf_walk_num #define elec_num chbrclf_elec_num @@ -6206,7 +7083,7 @@ assert( fabs(ao_vgl[26][4][224] - ( 3.153244195820293e-08)) < 1.e-14 ); } - #+end_src + #+end_src * End of files :noexport: diff --git a/org/qmckl_blas.org b/org/qmckl_blas.org index 9cd7e18..1cf76e4 100644 --- a/org/qmckl_blas.org +++ b/org/qmckl_blas.org @@ -84,8 +84,8 @@ are not intended to be passed to external codes. #+begin_src c :comments org :tangle (eval h_private_type) :exports none typedef struct qmckl_vector { - int64_t size; double* restrict data; + int64_t size; } qmckl_vector; #+end_src @@ -160,8 +160,8 @@ qmckl_vector_free( qmckl_context context, #+begin_src c :comments org :tangle (eval h_private_type) :exports none typedef struct qmckl_matrix { - int64_t size[2]; double* restrict data; + int64_t size[2]; } qmckl_matrix; #+end_src @@ -245,9 +245,9 @@ qmckl_matrix_free( qmckl_context context, #define QMCKL_TENSOR_ORDER_MAX 16 typedef struct qmckl_tensor { + double* restrict data; int64_t order; int64_t size[QMCKL_TENSOR_ORDER_MAX]; - double* restrict data; } qmckl_tensor; #+end_src diff --git a/org/qmckl_context.org b/org/qmckl_context.org index ac389f0..9054ba1 100644 --- a/org/qmckl_context.org +++ b/org/qmckl_context.org @@ -169,7 +169,7 @@ qmckl_context qmckl_context_check(const qmckl_context context) { if (context == QMCKL_NULL_CONTEXT) return QMCKL_NULL_CONTEXT; - const qmckl_context_struct* const ctx = (const qmckl_context_struct*) context; + const qmckl_context_struct* const ctx = (qmckl_context_struct*) context; /* Try to access memory */ if (ctx->tag != VALID_TAG) { @@ -267,7 +267,7 @@ qmckl_context qmckl_context_create() { { ctx->tag = VALID_TAG; - const qmckl_context context = (const qmckl_context) ctx; + const qmckl_context context = (qmckl_context) ctx; assert ( qmckl_context_check(context) != QMCKL_NULL_CONTEXT ); qmckl_exit_code rc; diff --git a/org/qmckl_determinant.org b/org/qmckl_determinant.org index 0412db6..ec4c9d9 100644 --- a/org/qmckl_determinant.org +++ b/org/qmckl_determinant.org @@ -182,7 +182,7 @@ qmckl_exit_code qmckl_init_determinant(qmckl_context context) { return false; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); ctx->det.uninitialized = (1 << 6) - 1; @@ -216,7 +216,7 @@ bool qmckl_determinant_provided(const qmckl_context context) { return false; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); return ctx->det.provided; @@ -238,7 +238,7 @@ char qmckl_get_determinant_type (const qmckl_context context) { return (char) 0; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int32_t mask = 1; @@ -256,7 +256,7 @@ int64_t qmckl_get_determinant_walk_num (const qmckl_context context) { return (int64_t) 0; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int32_t mask = 1 << 1; @@ -274,7 +274,7 @@ int64_t qmckl_get_determinant_det_num_alpha (const qmckl_context context) { return (int64_t) 0; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int32_t mask = 1 << 2; @@ -292,7 +292,7 @@ int64_t qmckl_get_determinant_det_num_beta (const qmckl_context context) { return (int64_t) 0; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int32_t mask = 1 << 3; @@ -310,7 +310,7 @@ int64_t* qmckl_get_determinant_mo_index_alpha (const qmckl_context context) { return NULL; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int32_t mask = 1 << 4; @@ -328,7 +328,7 @@ int64_t* qmckl_get_determinant_mo_index_beta (const qmckl_context context) { return NULL; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int32_t mask = 1 << 5; @@ -363,7 +363,7 @@ if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } -qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; +qmckl_context_struct* const ctx = (qmckl_context_struct*) context; #+end_src #+NAME:post2 @@ -525,7 +525,7 @@ qmckl_exit_code qmckl_finalize_determinant(qmckl_context context) { NULL); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); qmckl_exit_code rc; @@ -596,7 +596,7 @@ qmckl_exit_code qmckl_get_det_vgl_alpha(qmckl_context context, double * const de rc = qmckl_provide_det_vgl_alpha(context); if (rc != QMCKL_SUCCESS) return rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); size_t sze = 5 * ctx->det.det_num_alpha * ctx->det.walk_num * @@ -623,7 +623,7 @@ qmckl_exit_code qmckl_get_det_vgl_beta(qmckl_context context, double * const det rc = qmckl_provide_det_vgl_beta(context); if (rc != QMCKL_SUCCESS) return rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); size_t sze = 5 * ctx->det.det_num_beta * ctx->det.walk_num * @@ -649,7 +649,7 @@ qmckl_exit_code qmckl_provide_det_vgl_alpha(qmckl_context context) { return QMCKL_NULL_CONTEXT; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); if(!(ctx->nucleus.provided)) { @@ -748,7 +748,7 @@ qmckl_exit_code qmckl_provide_det_vgl_beta(qmckl_context context) { return QMCKL_NULL_CONTEXT; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); if(!(ctx->nucleus.provided)) { @@ -1134,36 +1134,28 @@ end function qmckl_compute_det_vgl_beta_f #+begin_src c :tangle (eval c_test) :exports none -#define walk_num chbrclf_walk_num -#define elec_num chbrclf_elec_num -#define shell_num chbrclf_shell_num -#define ao_num chbrclf_ao_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]); -const int64_t nucl_num = chbrclf_nucl_num; const double* nucl_charge = chbrclf_charge; const double* nucl_coord = &(chbrclf_nucl_coord[0][0]); -rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num); +rc = qmckl_set_electron_num (context, chbrclf_elec_up_num, chbrclf_elec_dn_num); assert (rc == QMCKL_SUCCESS); -rc = qmckl_set_electron_walk_num (context, walk_num); +rc = qmckl_set_electron_walk_num (context, chbrclf_walk_num); assert (rc == QMCKL_SUCCESS); assert(qmckl_electron_provided(context)); -rc = qmckl_set_electron_coord (context, 'N', elec_coord, walk_num*elec_num*3); +rc = qmckl_set_electron_coord (context, 'N', elec_coord, chbrclf_walk_num*chbrclf_elec_num*3); assert(rc == QMCKL_SUCCESS); -rc = qmckl_set_nucleus_num (context, nucl_num); +rc = qmckl_set_nucleus_num (context, chbrclf_nucl_num); assert(rc == QMCKL_SUCCESS); -rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), nucl_num*3); +rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), chbrclf_nucl_num*3); assert(rc == QMCKL_SUCCESS); -rc = qmckl_set_nucleus_charge(context, nucl_charge, nucl_num); +rc = qmckl_set_nucleus_charge(context, nucl_charge, chbrclf_nucl_num); assert(rc == QMCKL_SUCCESS); assert(qmckl_nucleus_provided(context)); @@ -1195,27 +1187,27 @@ rc = qmckl_set_ao_basis_prim_num (context, chbrclf_prim_num); assert(rc == QMCKL_SUCCESS); assert(!qmckl_ao_basis_provided(context)); -rc = qmckl_set_ao_basis_nucleus_index (context, nucleus_index, nucl_num); +rc = qmckl_set_ao_basis_nucleus_index (context, nucleus_index, chbrclf_nucl_num); assert(rc == QMCKL_SUCCESS); assert(!qmckl_ao_basis_provided(context)); -rc = qmckl_set_ao_basis_nucleus_shell_num (context, nucleus_shell_num, nucl_num); +rc = qmckl_set_ao_basis_nucleus_shell_num (context, nucleus_shell_num, chbrclf_nucl_num); assert(rc == QMCKL_SUCCESS); assert(!qmckl_ao_basis_provided(context)); -rc = qmckl_set_ao_basis_shell_ang_mom (context, shell_ang_mom, shell_num); +rc = qmckl_set_ao_basis_shell_ang_mom (context, shell_ang_mom, chbrclf_shell_num); assert(rc == QMCKL_SUCCESS); assert(!qmckl_ao_basis_provided(context)); -rc = qmckl_set_ao_basis_shell_factor (context, shell_factor, shell_num); +rc = qmckl_set_ao_basis_shell_factor (context, shell_factor, chbrclf_shell_num); assert(rc == QMCKL_SUCCESS); assert(!qmckl_ao_basis_provided(context)); -rc = qmckl_set_ao_basis_shell_prim_num (context, shell_prim_num, shell_num); +rc = qmckl_set_ao_basis_shell_prim_num (context, shell_prim_num, chbrclf_shell_num); assert(rc == QMCKL_SUCCESS); assert(!qmckl_ao_basis_provided(context)); -rc = qmckl_set_ao_basis_shell_prim_index (context, shell_prim_index, shell_num); +rc = qmckl_set_ao_basis_shell_prim_index (context, shell_prim_index, chbrclf_shell_num); assert(rc == QMCKL_SUCCESS); assert(!qmckl_ao_basis_provided(context)); @@ -1239,14 +1231,13 @@ assert(rc == QMCKL_SUCCESS); assert(qmckl_ao_basis_provided(context)); -double ao_vgl[walk_num*elec_num][5][chbrclf_ao_num]; +double ao_vgl[chbrclf_walk_num*chbrclf_elec_num][5][chbrclf_ao_num]; -rc = qmckl_get_ao_basis_ao_vgl(context, &(ao_vgl[0][0][0]), (int64_t) 5*walk_num*elec_num*chbrclf_ao_num); +rc = qmckl_get_ao_basis_ao_vgl(context, &(ao_vgl[0][0][0]), (int64_t) 5*chbrclf_walk_num*chbrclf_elec_num*chbrclf_ao_num); assert (rc == QMCKL_SUCCESS); /* Set up MO data */ -const int64_t mo_num = chbrclf_mo_num; -rc = qmckl_set_mo_basis_mo_num(context, mo_num); +rc = qmckl_set_mo_basis_mo_num(context, chbrclf_mo_num); assert (rc == QMCKL_SUCCESS); const double * mo_coefficient = &(chbrclf_mo_coef[0]); @@ -1256,31 +1247,31 @@ assert (rc == QMCKL_SUCCESS); assert(qmckl_mo_basis_provided(context)); -double mo_vgl[walk_num*elec_num][5][chbrclf_mo_num]; -rc = qmckl_get_mo_basis_mo_vgl(context, &(mo_vgl[0][0][0]), 5*walk_num*elec_num*chbrclf_mo_num); +double mo_vgl[chbrclf_walk_num*chbrclf_elec_num][5][chbrclf_mo_num]; +rc = qmckl_get_mo_basis_mo_vgl(context, &(mo_vgl[0][0][0]), 5*chbrclf_walk_num*chbrclf_elec_num*chbrclf_mo_num); assert (rc == QMCKL_SUCCESS); /* Set up determinant data */ -const int64_t det_num_alpha = 1; -const int64_t det_num_beta = 1; -int64_t mo_index_alpha[det_num_alpha][walk_num][elec_up_num]; -int64_t mo_index_beta[det_num_alpha][walk_num][elec_dn_num]; +#define det_num_alpha 1 +#define det_num_beta 1 +int64_t mo_index_alpha[det_num_alpha][chbrclf_walk_num][chbrclf_elec_up_num]; +int64_t mo_index_beta[det_num_alpha][chbrclf_walk_num][chbrclf_elec_dn_num]; int i, j, k; for(k = 0; k < det_num_alpha; ++k) - for(i = 0; i < walk_num; ++i) - for(j = 0; j < elec_up_num; ++j) + for(i = 0; i < chbrclf_walk_num; ++i) + for(j = 0; j < chbrclf_elec_up_num; ++j) mo_index_alpha[k][i][j] = j + 1; for(k = 0; k < det_num_beta; ++k) - for(i = 0; i < walk_num; ++i) - for(j = 0; j < elec_up_num; ++j) + for(i = 0; i < chbrclf_walk_num; ++i) + for(j = 0; j < chbrclf_elec_up_num; ++j) mo_index_beta[k][i][j] = j + 1; rc = qmckl_set_determinant_type (context, typ); assert(rc == QMCKL_SUCCESS); -rc = qmckl_set_determinant_walk_num (context, walk_num); +rc = qmckl_set_determinant_walk_num (context, chbrclf_walk_num); assert (rc == QMCKL_SUCCESS); rc = qmckl_set_determinant_det_num_alpha (context, det_num_alpha); @@ -1297,8 +1288,8 @@ assert (rc == QMCKL_SUCCESS); // Get slater-determinant -double det_vgl_alpha[det_num_alpha][walk_num][5][elec_up_num][elec_up_num]; -double det_vgl_beta[det_num_beta][walk_num][5][elec_dn_num][elec_dn_num]; +double det_vgl_alpha[det_num_alpha][chbrclf_walk_num][5][chbrclf_elec_up_num][chbrclf_elec_up_num]; +double det_vgl_beta[det_num_beta][chbrclf_walk_num][5][chbrclf_elec_dn_num][chbrclf_elec_dn_num]; rc = qmckl_get_det_vgl_alpha(context, &(det_vgl_alpha[0][0][0][0][0])); assert (rc == QMCKL_SUCCESS); @@ -1347,7 +1338,7 @@ qmckl_exit_code qmckl_get_det_inv_matrix_alpha(qmckl_context context, double * c rc = qmckl_provide_det_inv_matrix_alpha(context); if (rc != QMCKL_SUCCESS) return rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); size_t sze = ctx->det.det_num_alpha * ctx->det.walk_num * ctx->electron.up_num * ctx->electron.up_num; @@ -1376,7 +1367,7 @@ qmckl_exit_code qmckl_get_det_inv_matrix_beta(qmckl_context context, double * co rc = qmckl_provide_det_inv_matrix_beta(context); if (rc != QMCKL_SUCCESS) return rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); size_t sze = ctx->det.det_num_alpha * ctx->det.walk_num * ctx->electron.down_num * ctx->electron.down_num; @@ -1405,7 +1396,7 @@ qmckl_exit_code qmckl_get_det_adj_matrix_alpha(qmckl_context context, double * c rc = qmckl_provide_det_inv_matrix_alpha(context); if (rc != QMCKL_SUCCESS) return rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); size_t sze = ctx->det.det_num_alpha * ctx->det.walk_num * ctx->electron.up_num * ctx->electron.up_num; @@ -1434,7 +1425,7 @@ qmckl_exit_code qmckl_get_det_adj_matrix_beta(qmckl_context context, double * co rc = qmckl_provide_det_inv_matrix_beta(context); if (rc != QMCKL_SUCCESS) return rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); size_t sze = ctx->det.det_num_alpha * ctx->det.walk_num * ctx->electron.down_num * ctx->electron.down_num; @@ -1463,7 +1454,7 @@ qmckl_exit_code qmckl_get_det_alpha(qmckl_context context, double * const det_va rc = qmckl_provide_det_inv_matrix_alpha(context); if (rc != QMCKL_SUCCESS) return rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); size_t sze = ctx->det.det_num_alpha * ctx->det.walk_num; @@ -1492,7 +1483,7 @@ qmckl_exit_code qmckl_get_det_beta(qmckl_context context, double * const det_val rc = qmckl_provide_det_inv_matrix_beta(context); if (rc != QMCKL_SUCCESS) return rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); size_t sze = ctx->det.det_num_alpha * ctx->det.walk_num; @@ -1517,7 +1508,7 @@ qmckl_exit_code qmckl_provide_det_inv_matrix_alpha(qmckl_context context) { return QMCKL_NULL_CONTEXT; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); if(!(ctx->nucleus.provided)) { @@ -1640,7 +1631,7 @@ qmckl_exit_code qmckl_provide_det_inv_matrix_beta(qmckl_context context) { return QMCKL_NULL_CONTEXT; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); if(!(ctx->nucleus.provided)) { @@ -2047,8 +2038,8 @@ end function qmckl_compute_det_inv_matrix_beta_f #+begin_src c :tangle (eval c_test) :exports none // Get adjoint of the slater-determinant -double det_inv_matrix_alpha[det_num_alpha][walk_num][elec_up_num][elec_up_num]; -double det_inv_matrix_beta[det_num_beta][walk_num][elec_dn_num][elec_dn_num]; +double det_inv_matrix_alpha[det_num_alpha][chbrclf_walk_num][chbrclf_elec_up_num][chbrclf_elec_up_num]; +double det_inv_matrix_beta[det_num_beta][chbrclf_walk_num][chbrclf_elec_dn_num][chbrclf_elec_dn_num]; rc = qmckl_get_det_inv_matrix_alpha(context, &(det_inv_matrix_alpha[0][0][0][0])); assert (rc == QMCKL_SUCCESS); diff --git a/org/qmckl_electron.org b/org/qmckl_electron.org index 33ac366..c0c0254 100644 --- a/org/qmckl_electron.org +++ b/org/qmckl_electron.org @@ -97,8 +97,8 @@ int main() { | ~ee_distance_rescaled_deriv_e_date~ | ~uint64_t~ | Last modification date of the electron-electron distance derivatives | | ~ee_pot~ | ~double[walk_num]~ | Electron-electron rescaled distances derivatives | | ~ee_pot_date~ | ~uint64_t~ | Last modification date of the electron-electron distance derivatives | - | ~en_pot~ | double[walk_num] | Electron-nucleus potential energy | - | ~en_pot_date~ | int64_t | Date when the electron-nucleus potential energy was computed | + | ~en_pot~ | ~double[walk_num]~ | Electron-nucleus potential energy | + | ~en_pot_date~ | ~int64_t~ | Date when the electron-nucleus potential energy was computed | | ~en_distance_rescaled~ | ~double[walk_num][nucl_num][num]~ | Electron-nucleus distances | | ~en_distance_rescaled_date~ | ~uint64_t~ | Last modification date of the electron-electron distances | | ~en_distance_rescaled_deriv_e~ | ~double[walk_num][4][nucl_num][num]~ | Electron-electron rescaled distances derivatives | @@ -157,7 +157,7 @@ qmckl_exit_code qmckl_init_electron(qmckl_context context) { return false; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); ctx->electron.uninitialized = (1 << 2) - 1; @@ -182,7 +182,7 @@ bool qmckl_electron_provided(const qmckl_context context) { return false; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); return ctx->electron.provided; @@ -228,7 +228,7 @@ qmckl_get_electron_num (const qmckl_context context, int64_t* const num) { "num is a null pointer"); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int32_t mask = 1 << 0; @@ -256,7 +256,7 @@ qmckl_get_electron_up_num (const qmckl_context context, int64_t* const up_num) { "up_num is a null pointer"); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int32_t mask = 1 << 0; @@ -284,7 +284,7 @@ qmckl_get_electron_down_num (const qmckl_context context, int64_t* const down_nu "down_num is a null pointer"); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int32_t mask = 1 << 0; @@ -323,7 +323,7 @@ qmckl_get_electron_walk_num (const qmckl_context context, int64_t* const walk_nu "walk_num is a null pointer"); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int32_t mask = 1 << 1; @@ -360,7 +360,7 @@ qmckl_get_electron_rescale_factor_ee (const qmckl_context context, double* const "rescale_factor_kappa_ee is a null pointer"); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); assert (ctx->electron.rescale_factor_kappa_ee > 0.0); @@ -383,7 +383,7 @@ qmckl_get_electron_rescale_factor_en (const qmckl_context context, double* const "rescale_factor_kappa_en is a null pointer"); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); assert (ctx->electron.rescale_factor_kappa_en > 0.0); @@ -448,7 +448,7 @@ qmckl_get_electron_coord (const qmckl_context context, return QMCKL_INVALID_CONTEXT; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); if (!ctx->electron.provided) { @@ -489,7 +489,7 @@ if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } -qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; +qmckl_context_struct* const ctx = (qmckl_context_struct*) context; #+end_src #+NAME:post2 @@ -897,7 +897,7 @@ qmckl_exit_code qmckl_get_electron_ee_distance(qmckl_context context, double* co rc = qmckl_provide_ee_distance(context); if (rc != QMCKL_SUCCESS) return rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); size_t sze = ctx->electron.num * ctx->electron.num * ctx->electron.walk_num; @@ -921,7 +921,7 @@ qmckl_exit_code qmckl_provide_ee_distance(qmckl_context context) return QMCKL_NULL_CONTEXT; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); @@ -1138,7 +1138,7 @@ qmckl_exit_code qmckl_get_electron_ee_distance_rescaled(qmckl_context context, d rc = qmckl_provide_ee_distance_rescaled(context); if (rc != QMCKL_SUCCESS) return rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); size_t sze = ctx->electron.num * ctx->electron.num * ctx->electron.walk_num; @@ -1162,7 +1162,7 @@ qmckl_exit_code qmckl_provide_ee_distance_rescaled(qmckl_context context) return QMCKL_NULL_CONTEXT; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); @@ -1218,7 +1218,7 @@ qmckl_exit_code qmckl_provide_ee_distance_rescaled(qmckl_context context) | ~elec_num~ | ~int64_t~ | in | Number of electrons | | ~rescale_factor_kappa_ee~ | ~double~ | in | Factor to rescale ee distances | | ~walk_num~ | ~int64_t~ | in | Number of walkers | - | ~coord~ | ~double[walk_num][3][elec_num]~ | in | Electron coordinates | + | ~coord~ | ~double[3][walk_num][elec_num]~ | in | Electron coordinates | | ~ee_distance~ | ~double[walk_num][elec_num][elec_num]~ | out | Electron-electron rescaled distances | #+begin_src f90 :comments org :tangle (eval f) :noweb yes @@ -1231,7 +1231,7 @@ integer function qmckl_compute_ee_distance_rescaled_f(context, elec_num, rescale integer*8 , intent(in) :: elec_num double precision , intent(in) :: rescale_factor_kappa_ee integer*8 , intent(in) :: walk_num - double precision , intent(in) :: coord(elec_num,3,walk_num) + double precision , intent(in) :: coord(elec_num,walk_num,3) double precision , intent(out) :: ee_distance_rescaled(elec_num,elec_num,walk_num) integer*8 :: k @@ -1357,7 +1357,7 @@ assert(fabs(ee_distance_rescaled[elec_num*elec_num+1]-0.9985724058042633) < 1.e- #+end_src -** Electron-electron rescaled distance gradients and laplacian with respect to electron coords +** Electron-electron rescaled distance gradients and Laplacian with respect to electron coords The rescaled distances which is given as $R = (1 - \exp{-\kappa r})/\kappa$ needs to be perturbed with respect to the electorn coordinates. @@ -1384,7 +1384,7 @@ qmckl_exit_code qmckl_get_electron_ee_distance_rescaled_deriv_e(qmckl_context co rc = qmckl_provide_ee_distance_rescaled_deriv_e(context); if (rc != QMCKL_SUCCESS) return rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); size_t sze = 4 * ctx->electron.num * ctx->electron.num * ctx->electron.walk_num; @@ -1408,7 +1408,7 @@ qmckl_exit_code qmckl_provide_ee_distance_rescaled_deriv_e(qmckl_context context return QMCKL_NULL_CONTEXT; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); @@ -1464,7 +1464,7 @@ qmckl_exit_code qmckl_provide_ee_distance_rescaled_deriv_e(qmckl_context context | ~elec_num~ | ~int64_t~ | in | Number of electrons | | ~rescale_factor_kappa_ee~ | ~double~ | in | Factor to rescale ee distances | | ~walk_num~ | ~int64_t~ | in | Number of walkers | - | ~coord~ | ~double[walk_num][3][elec_num]~ | in | Electron coordinates | + | ~coord~ | ~double[3][walk_num][elec_num]~ | in | Electron coordinates | | ~ee_distance_deriv_e~ | ~double[walk_num][4][elec_num][elec_num]~ | out | Electron-electron rescaled distance derivatives | #+begin_src f90 :comments org :tangle (eval f) :noweb yes @@ -1477,7 +1477,7 @@ integer function qmckl_compute_ee_distance_rescaled_deriv_e_f(context, elec_num, integer*8 , intent(in) :: elec_num double precision , intent(in) :: rescale_factor_kappa_ee integer*8 , intent(in) :: walk_num - double precision , intent(in) :: coord(elec_num,3,walk_num) + double precision , intent(in) :: coord(elec_num,walk_num,3) double precision , intent(out) :: ee_distance_rescaled_deriv_e(4,elec_num,elec_num,walk_num) integer*8 :: k @@ -1501,8 +1501,8 @@ integer function qmckl_compute_ee_distance_rescaled_deriv_e_f(context, elec_num, do k=1,walk_num info = qmckl_distance_rescaled_deriv_e(context, 'T', 'T', elec_num, elec_num, & - coord(1,1,k), elec_num, & - coord(1,1,k), elec_num, & + coord(1,k,1), elec_num*walk_num, & + coord(1,k,1), elec_num*walk_num, & ee_distance_rescaled_deriv_e(1,1,1,k), elec_num, rescale_factor_kappa_ee) if (info /= QMCKL_SUCCESS) then exit @@ -1613,7 +1613,7 @@ qmckl_exit_code qmckl_get_electron_ee_potential(qmckl_context context, double* c rc = qmckl_provide_ee_potential(context); if (rc != QMCKL_SUCCESS) return rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); size_t sze = ctx->electron.walk_num * sizeof(double); @@ -1637,7 +1637,7 @@ qmckl_exit_code qmckl_provide_ee_potential(qmckl_context context) return QMCKL_NULL_CONTEXT; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); if (!ctx->electron.provided) return QMCKL_NOT_PROVIDED; @@ -1818,7 +1818,7 @@ qmckl_exit_code qmckl_get_electron_en_distance(qmckl_context context, double* di rc = qmckl_provide_en_distance(context); if (rc != QMCKL_SUCCESS) return rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); size_t sze = ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num; @@ -1842,7 +1842,7 @@ qmckl_exit_code qmckl_provide_en_distance(qmckl_context context) return QMCKL_NULL_CONTEXT; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); if (!(ctx->nucleus.provided)) { @@ -1905,7 +1905,7 @@ qmckl_exit_code qmckl_provide_en_distance(qmckl_context context) | ~elec_num~ | ~int64_t~ | in | Number of electrons | | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | | ~walk_num~ | ~int64_t~ | in | Number of walkers | - | ~elec_coord~ | ~double[walk_num][3][elec_num]~ | in | Electron coordinates | + | ~elec_coord~ | ~double[3][walk_num][elec_num]~ | in | Electron coordinates | | ~nucl_coord~ | ~double[3][elec_num]~ | in | Nuclear coordinates | | ~en_distance~ | ~double[walk_num][nucl_num][elec_num]~ | out | Electron-nucleus distances | @@ -2097,7 +2097,7 @@ qmckl_exit_code qmckl_get_electron_en_distance_rescaled(qmckl_context context, d rc = qmckl_provide_en_distance_rescaled(context); if (rc != QMCKL_SUCCESS) return rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); size_t sze = ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num; @@ -2121,7 +2121,7 @@ qmckl_exit_code qmckl_provide_en_distance_rescaled(qmckl_context context) return QMCKL_NULL_CONTEXT; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); if (!(ctx->nucleus.provided)) { @@ -2183,7 +2183,7 @@ qmckl_exit_code qmckl_provide_en_distance_rescaled(qmckl_context context) | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | | ~rescale_factor_kappa_en~ | ~double~ | in | The factor for rescaled distances | | ~walk_num~ | ~int64_t~ | in | Number of walkers | - | ~elec_coord~ | ~double[walk_num][3][elec_num]~ | in | Electron coordinates | + | ~elec_coord~ | ~double[3][walk_num][elec_num]~ | in | Electron coordinates | | ~nucl_coord~ | ~double[3][elec_num]~ | in | Nuclear coordinates | | ~en_distance_rescaled~ | ~double[walk_num][nucl_num][elec_num]~ | out | Electron-nucleus distances | @@ -2385,7 +2385,7 @@ qmckl_exit_code qmckl_get_electron_en_distance_rescaled_deriv_e(qmckl_context co rc = qmckl_provide_en_distance_rescaled_deriv_e(context); if (rc != QMCKL_SUCCESS) return rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); size_t sze = 4 * ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num; @@ -2409,7 +2409,7 @@ qmckl_exit_code qmckl_provide_en_distance_rescaled_deriv_e(qmckl_context context return QMCKL_NULL_CONTEXT; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); if (!(ctx->nucleus.provided)) { @@ -2471,9 +2471,9 @@ qmckl_exit_code qmckl_provide_en_distance_rescaled_deriv_e(qmckl_context context | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | | ~rescale_factor_kappa_en~ | ~double~ | in | The factor for rescaled distances | | ~walk_num~ | ~int64_t~ | in | Number of walkers | - | ~elec_coord~ | ~double[walk_num][3][elec_num]~ | in | Electron coordinates | + | ~elec_coord~ | ~double[3][walk_num][elec_num]~ | in | Electron coordinates | | ~nucl_coord~ | ~double[3][elec_num]~ | in | Nuclear coordinates | - | ~en_distance_rescaled_deriv_e~ | ~double[walk_num][4][nucl_num][elec_num]~ | out | Electron-nucleus distance derivatives | + | ~en_distance_rescaled_deriv_e~ | ~double[walk_num][nucl_num][elec_num][4]~ | out | Electron-nucleus distance derivatives | #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_en_distance_rescaled_deriv_e_f(context, elec_num, nucl_num, & @@ -2656,7 +2656,7 @@ qmckl_exit_code qmckl_get_electron_en_potential(qmckl_context context, double* c rc = qmckl_provide_en_potential(context); if (rc != QMCKL_SUCCESS) return rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); size_t sze = ctx->electron.walk_num * sizeof(double); @@ -2680,7 +2680,7 @@ qmckl_exit_code qmckl_provide_en_potential(qmckl_context context) return QMCKL_NULL_CONTEXT; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); if (!ctx->electron.provided) return QMCKL_NOT_PROVIDED; @@ -2843,7 +2843,7 @@ assert (rc == QMCKL_SUCCESS); *** Compute :noexport: - # begin_src f90 :comments org :tangle (eval f) :noweb yes + # begin_src f90 :comments org :tangle (eval f) :noweb yes subroutine draw_init_points implicit none BEGIN_DOC diff --git a/org/qmckl_error.org b/org/qmckl_error.org index 97cc7cc..d23fa8e 100644 --- a/org/qmckl_error.org +++ b/org/qmckl_error.org @@ -239,7 +239,7 @@ for (text, code, message) in table: message = message.replace("'",'"') result += [ f"""case {text}: return {message}; - break;""" ] + """ ] return '\n'.join(result) #+end_src @@ -247,89 +247,91 @@ return '\n'.join(result) #+RESULTS: cases #+begin_example case QMCKL_SUCCESS: - return "Success"; - break; + return "Success"; + case QMCKL_INVALID_ARG_1: - return "Invalid argument 1"; - break; + return "Invalid argument 1"; + case QMCKL_INVALID_ARG_2: - return "Invalid argument 2"; - break; + return "Invalid argument 2"; + case QMCKL_INVALID_ARG_3: - return "Invalid argument 3"; - break; + return "Invalid argument 3"; + case QMCKL_INVALID_ARG_4: - return "Invalid argument 4"; - break; + return "Invalid argument 4"; + case QMCKL_INVALID_ARG_5: - return "Invalid argument 5"; - break; + return "Invalid argument 5"; + case QMCKL_INVALID_ARG_6: - return "Invalid argument 6"; - break; + return "Invalid argument 6"; + case QMCKL_INVALID_ARG_7: - return "Invalid argument 7"; - break; + return "Invalid argument 7"; + case QMCKL_INVALID_ARG_8: - return "Invalid argument 8"; - break; + return "Invalid argument 8"; + case QMCKL_INVALID_ARG_9: - return "Invalid argument 9"; - break; + return "Invalid argument 9"; + case QMCKL_INVALID_ARG_10: - return "Invalid argument 10"; - break; + return "Invalid argument 10"; + case QMCKL_INVALID_ARG_11: - return "Invalid argument 11"; - break; + return "Invalid argument 11"; + case QMCKL_INVALID_ARG_12: - return "Invalid argument 12"; - break; + return "Invalid argument 12"; + case QMCKL_INVALID_ARG_13: - return "Invalid argument 13"; - break; + return "Invalid argument 13"; + case QMCKL_INVALID_ARG_14: - return "Invalid argument 14"; - break; + return "Invalid argument 14"; + case QMCKL_INVALID_ARG_15: - return "Invalid argument 15"; - break; + return "Invalid argument 15"; + case QMCKL_INVALID_ARG_16: - return "Invalid argument 16"; - break; + return "Invalid argument 16"; + case QMCKL_INVALID_ARG_17: - return "Invalid argument 17"; - break; + return "Invalid argument 17"; + case QMCKL_INVALID_ARG_18: - return "Invalid argument 18"; - break; + return "Invalid argument 18"; + case QMCKL_INVALID_ARG_19: - return "Invalid argument 19"; - break; + return "Invalid argument 19"; + case QMCKL_INVALID_ARG_20: - return "Invalid argument 20"; - break; + return "Invalid argument 20"; + case QMCKL_FAILURE: - return "Failure"; - break; + return "Failure"; + case QMCKL_ERRNO: - return strerror(errno); - break; + return strerror(errno); + case QMCKL_INVALID_CONTEXT: - return "Invalid context"; - break; + return "Invalid context"; + case QMCKL_ALLOCATION_FAILED: - return "Allocation failed"; - break; + return "Allocation failed"; + case QMCKL_DEALLOCATION_FAILED: - return "De-allocation failed"; - break; + return "De-allocation failed"; + case QMCKL_NOT_PROVIDED: - return "Not provided"; - break; + return "Not provided"; + + case QMCKL_OUT_OF_BOUNDS: + return "Index out of bounds"; + case QMCKL_INVALID_EXIT_CODE: - return "Invalid exit code"; - break; + return "Invalid exit code"; #+end_example # Source @@ -414,7 +416,7 @@ qmckl_set_error(qmckl_context context, qmckl_lock(context); { - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); /* Impossible because the context is valid. */ ctx->error.exit_code = exit_code; @@ -460,7 +462,7 @@ qmckl_get_error(qmckl_context context, qmckl_lock(context); { - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); /* Impossible because the context is valid. */ /* Turn off annoying GCC warning */ diff --git a/org/qmckl_jastrow.org b/org/qmckl_jastrow.org index 61062af..3007181 100644 --- a/org/qmckl_jastrow.org +++ b/org/qmckl_jastrow.org @@ -11,7 +11,7 @@ \[ J(\mathbf{r},\mathbf{R}) = J_{\text{eN}}(\mathbf{r},\mathbf{R}) + J_{\text{ee}}(\mathbf{r}) + J_{\text{eeN}}(\mathbf{r},\mathbf{R}) \] - + In the following, we us the notations $r_{ij} = |\mathbf{r}_i - \mathbf{r}_j|$ and $R_{i\alpha} = |\mathbf{r}_i - \mathbf{R}_\alpha|$. @@ -58,7 +58,6 @@ The terms $J_{\text{ee}}^\infty$ and $J_{\text{eN}}^\infty$ are shifts to ensure that $J_{\text{ee}}$ and $J_{\text{eN}}$ have an asymptotic value of zero. - * Headers :noexport: #+begin_src elisp :noexport :results none (org-babel-lob-ingest "../tools/lib.org") @@ -108,6 +107,7 @@ int main() { #include #include + #include #include "qmckl.h" @@ -116,6 +116,12 @@ int main() { #include "qmckl_memory_private_func.h" #include "qmckl_jastrow_private_func.h" #include "qmckl_jastrow_private_type.h" + +#ifdef HAVE_CUBLAS_OFFLOAD +#include "cublas_v2.h" +#endif + + #+end_src * Context @@ -154,25 +160,25 @@ int main() { computed data: - | Variable | Type | In/Out | Description | - |----------------------------+---------------------------------------------------------------------+-------------------------------------------------+---------------------------------| - | ~dim_cord_vect~ | ~int64_t~ | Number of unique C coefficients | | - | ~dim_cord_vect_date~ | ~uint64_t~ | Number of unique C coefficients | | - | ~asymp_jasb~ | ~double[2]~ | Asymptotic component | | - | ~asymp_jasb_date~ | ~uint64_t~ | Asymptotic component | | - | ~cord_vect_full~ | ~double[dim_cord_vect][nucl_num]~ | vector of non-zero coefficients | | - | ~cord_vect_full_date~ | ~uint64_t~ | Keep track of changes here | | - | ~lkpm_combined_index~ | ~int64_t[4][dim_cord_vect]~ | Transform l,k,p, and m into consecutive indices | | - | ~lkpm_combined_index_date~ | ~uint64_t~ | Transform l,k,p, and m into consecutive indices | | - | ~tmp_c~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]~ | out | vector of non-zero coefficients | - | ~dtmp_c~ | ~double[walk_num][elec_num][4][nucl_num][0:cord_num][0:cord_num-1]~ | vector of non-zero coefficients | | + | Variable | Type | In/Out | + |----------------------------+-----------------------------------------------------------------+-------------------------------------------------| + | ~dim_cord_vect~ | ~int64_t~ | Number of unique C coefficients | + | ~dim_cord_vect_date~ | ~uint64_t~ | Number of unique C coefficients | + | ~asymp_jasb~ | ~double[2]~ | Asymptotic component | + | ~asymp_jasb_date~ | ~uint64_t~ | Asymptotic component | + | ~cord_vect_full~ | ~double[dim_cord_vect][nucl_num]~ | vector of non-zero coefficients | + | ~cord_vect_full_date~ | ~uint64_t~ | Keep track of changes here | + | ~lkpm_combined_index~ | ~int64_t[4][dim_cord_vect]~ | Transform l,k,p, and m into consecutive indices | + | ~lkpm_combined_index_date~ | ~uint64_t~ | Transform l,k,p, and m into consecutive indices | + | ~tmp_c~ | ~double[walk_num][cord_num][cord_num+1][nucl_num][elec_num]~ | vector of non-zero coefficients | + | ~dtmp_c~ | ~double[walk_num][elec_num][4][nucl_num][cord_num+1][cord_num]~ | vector of non-zero coefficients | - | ~een_rescaled_n~ | ~double[walk_num][elec_num][nucl_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord | | - | ~een_rescaled_n_date~ | ~uint64_t~ | Keep track of the date of creation | | - | ~een_rescaled_e_deriv_e~ | ~double[walk_num][elec_num][4][elec_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons | | - | ~een_rescaled_e_deriv_e_date~ | ~uint64_t~ | Keep track of the date of creation | | - | ~een_rescaled_n_deriv_e~ | ~double[walk_num][elec_num][4][nucl_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons | | - | ~een_rescaled_n_deriv_e_date~ | ~uint64_t~ | Keep track of the date of creation | | + | ~een_rescaled_n~ | ~double[walk_num][cord_num+1][nucl_num][elec_num]~ | The electron-electron rescaled distances raised to the powers defined by cord | + | ~een_rescaled_n_date~ | ~uint64_t~ | Keep track of the date of creation | + | ~een_rescaled_e_deriv_e~ | ~double[walk_num][cord_num+1][elec_num][4][elec_num]~ | The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons | + | ~een_rescaled_e_deriv_e_date~ | ~uint64_t~ | Keep track of the date of creation | + | ~een_rescaled_n_deriv_e~ | ~double[walk_num][cord_num+1][nucl_num][4][elec_num]~ | The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons | + | ~een_rescaled_n_deriv_e_date~ | ~uint64_t~ | Keep track of the date of creation | #+NAME: jastrow_data #+BEGIN_SRC python :results none :exports none @@ -372,6 +378,10 @@ typedef struct qmckl_jastrow_struct{ uint64_t een_rescaled_n_deriv_e_date; bool provided; char * type; + + #ifdef HAVE_HPC + bool gpu_offload; + #endif } qmckl_jastrow_struct; #+end_src @@ -394,7 +404,7 @@ qmckl_exit_code qmckl_init_jastrow(qmckl_context context) { return false; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); ctx->jastrow.uninitialized = (1 << 5) - 1; @@ -436,7 +446,7 @@ bool qmckl_jastrow_provided(const qmckl_context context) { return false; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); return ctx->jastrow.provided; @@ -464,7 +474,7 @@ qmckl_exit_code qmckl_get_jastrow_aord_num (const qmckl_context context, int64_t "aord_num is a null pointer"); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int32_t mask = 1 << 0; @@ -491,7 +501,7 @@ qmckl_exit_code qmckl_get_jastrow_bord_num (const qmckl_context context, int64_t "aord_num is a null pointer"); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int32_t mask = 1 << 0; @@ -518,7 +528,7 @@ qmckl_exit_code qmckl_get_jastrow_cord_num (const qmckl_context context, int64_t "aord_num is a null pointer"); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int32_t mask = 1 << 0; @@ -545,7 +555,7 @@ qmckl_exit_code qmckl_get_jastrow_type_nucl_num (const qmckl_context context, in "type_nucl_num is a null pointer"); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int32_t mask = 1 << 1; @@ -576,7 +586,7 @@ qmckl_get_jastrow_type_nucl_vector (const qmckl_context context, "type_nucl_vector is a null pointer"); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int32_t mask = 1 << 2; @@ -613,7 +623,7 @@ qmckl_get_jastrow_aord_vector (const qmckl_context context, "aord_vector is a null pointer"); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int32_t mask = 1 << 3; @@ -650,7 +660,7 @@ qmckl_get_jastrow_bord_vector (const qmckl_context context, "bord_vector is a null pointer"); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int32_t mask = 1 << 4; @@ -687,7 +697,7 @@ qmckl_get_jastrow_cord_vector (const qmckl_context context, "cord_vector is a null pointer"); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int32_t mask = 1 << 5; @@ -735,7 +745,7 @@ if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } -qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; +qmckl_context_struct* const ctx = (qmckl_context_struct*) context; #+end_src #+NAME:post2 @@ -837,7 +847,7 @@ qmckl_set_jastrow_type_nucl_vector(qmckl_context context, } if (ctx->jastrow.type_nucl_vector != NULL) { - qmckl_exit_code rc = qmckl_free(context, ctx->jastrow.type_nucl_vector); + rc = qmckl_free(context, ctx->jastrow.type_nucl_vector); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_set_type_nucl_vector", @@ -896,7 +906,7 @@ qmckl_set_jastrow_aord_vector(qmckl_context context, } if (ctx->jastrow.aord_vector != NULL) { - qmckl_exit_code rc = qmckl_free(context, ctx->jastrow.aord_vector); + rc = qmckl_free(context, ctx->jastrow.aord_vector); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_set_ord_vector", @@ -959,7 +969,7 @@ qmckl_set_jastrow_bord_vector(qmckl_context context, } if (ctx->jastrow.bord_vector != NULL) { - qmckl_exit_code rc = qmckl_free(context, ctx->jastrow.bord_vector); + rc = qmckl_free(context, ctx->jastrow.bord_vector); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_set_ord_vector", @@ -1029,7 +1039,7 @@ qmckl_set_jastrow_cord_vector(qmckl_context context, } if (ctx->jastrow.cord_vector != NULL) { - qmckl_exit_code rc = qmckl_free(context, ctx->jastrow.cord_vector); + rc = qmckl_free(context, ctx->jastrow.cord_vector); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_set_ord_vector", @@ -1080,7 +1090,7 @@ qmckl_exit_code qmckl_finalize_jastrow(qmckl_context context) { return QMCKL_INVALID_CONTEXT; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); /* ----------------------------------- */ @@ -1109,6 +1119,11 @@ qmckl_exit_code qmckl_finalize_jastrow(qmckl_context context) { NULL); } + /* Decide if the Jastrow if offloaded on GPU or not */ +#if defined(HAVE_HPC) && (defined(HAVE_CUBLAS_OFFLOAD) || defined(HAVE_OPENACC_OFFLOAD) || defined(HAVE_OPENMP_OFFLOAD)) + ctx->jastrow.gpu_offload = true; // ctx->electron.num > 100; +#endif + qmckl_exit_code rc = QMCKL_SUCCESS; return rc; @@ -1330,7 +1345,7 @@ qmckl_get_jastrow_asymp_jasb(qmckl_context context, rc = qmckl_provide_asymp_jasb(context); if (rc != QMCKL_SUCCESS) return rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int64_t sze = 2; @@ -1361,7 +1376,7 @@ qmckl_exit_code qmckl_provide_asymp_jasb(qmckl_context context) return QMCKL_NULL_CONTEXT; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); /* Check if ee kappa is provided */ @@ -1388,12 +1403,11 @@ qmckl_exit_code qmckl_provide_asymp_jasb(qmckl_context context) ctx->jastrow.asymp_jasb = asymp_jasb; } - qmckl_exit_code rc = - qmckl_compute_asymp_jasb(context, - ctx->jastrow.bord_num, - ctx->jastrow.bord_vector, - rescale_factor_kappa_ee, - ctx->jastrow.asymp_jasb); + rc = qmckl_compute_asymp_jasb(context, + ctx->jastrow.bord_num, + ctx->jastrow.bord_vector, + rescale_factor_kappa_ee, + ctx->jastrow.asymp_jasb); if (rc != QMCKL_SUCCESS) { return rc; } @@ -1470,10 +1484,6 @@ qmckl_exit_code qmckl_compute_asymp_jasb ( const double rescale_factor_kappa_ee, double* const asymp_jasb ) { - double kappa_inv, x, asym_one; - - kappa_inv = 1.0 / rescale_factor_kappa_ee; - if (context == QMCKL_NULL_CONTEXT){ return QMCKL_INVALID_CONTEXT; } @@ -1482,14 +1492,15 @@ qmckl_exit_code qmckl_compute_asymp_jasb ( return QMCKL_INVALID_ARG_2; } - asym_one = bord_vector[0] * kappa_inv / (1.0 + bord_vector[1] * kappa_inv); + const double kappa_inv = 1.0 / rescale_factor_kappa_ee; + const double asym_one = bord_vector[0] * kappa_inv / (1.0 + bord_vector[1] * kappa_inv); asymp_jasb[0] = asym_one; asymp_jasb[1] = 0.5 * asym_one; for (int i = 0 ; i <= 1; ++i) { - x = kappa_inv; + double x = kappa_inv; for (int p = 1; p < bord_num; ++p){ - x = x * kappa_inv; + x *= kappa_inv; asymp_jasb[i] = asymp_jasb[i] + bord_vector[p + 1] * x; } } @@ -1498,16 +1509,15 @@ qmckl_exit_code qmckl_compute_asymp_jasb ( } #+end_src - #+CALL: generate_c_header(table=qmckl_asymp_jasb_args,rettyp=get_value("CRetType"),fname=get_value("Name")) +# #+CALL: generate_c_header(table=qmckl_asymp_jasb_args,rettyp=get_value("CRetType"),fname=get_value("Name")) - #+RESULTS: - #+begin_src c :tangle (eval h_func) :comments org + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_asymp_jasb ( - const qmckl_context context, - const int64_t bord_num, - const double* bord_vector, - const double rescale_factor_kappa_ee, - double* const asymp_jasb ); + const qmckl_context context, + const int64_t bord_num, + const double* bord_vector, + const double rescale_factor_kappa_ee, + double* const asymp_jasb ); #+end_src @@ -1615,7 +1625,7 @@ qmckl_get_jastrow_factor_ee(qmckl_context context, rc = qmckl_provide_factor_ee(context); if (rc != QMCKL_SUCCESS) return rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int64_t sze=ctx->electron.walk_num; @@ -1646,7 +1656,7 @@ qmckl_exit_code qmckl_provide_factor_ee(qmckl_context context) return QMCKL_NULL_CONTEXT; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); /* Check if ee rescaled distance is provided */ @@ -1672,16 +1682,15 @@ qmckl_exit_code qmckl_provide_factor_ee(qmckl_context context) ctx->jastrow.factor_ee = factor_ee; } - qmckl_exit_code rc = - qmckl_compute_factor_ee(context, - ctx->electron.walk_num, - ctx->electron.num, - ctx->electron.up_num, - ctx->jastrow.bord_num, - ctx->jastrow.bord_vector, - ctx->electron.ee_distance_rescaled, - ctx->jastrow.asymp_jasb, - ctx->jastrow.factor_ee); + rc = qmckl_compute_factor_ee(context, + ctx->electron.walk_num, + ctx->electron.num, + ctx->electron.up_num, + ctx->jastrow.bord_num, + ctx->jastrow.bord_vector, + ctx->electron.ee_distance_rescaled, + ctx->jastrow.asymp_jasb, + ctx->jastrow.factor_ee); if (rc != QMCKL_SUCCESS) { return rc; } @@ -1727,7 +1736,7 @@ integer function qmckl_compute_factor_ee_f(context, walk_num, elec_num, up_num, double precision , intent(out) :: factor_ee(walk_num) integer*8 :: i, j, p, ipar, nw - double precision :: pow_ser, x, spin_fact, power_ser + double precision :: x, power_ser, spin_fact info = QMCKL_SUCCESS @@ -1766,7 +1775,7 @@ integer function qmckl_compute_factor_ee_f(context, walk_num, elec_num, up_num, power_ser = power_ser + bord_vector(p + 1) * x end do - if(j .LE. up_num .OR. i .GT. up_num) then + if(j <= up_num .OR. i > up_num) then spin_fact = 0.5d0 ipar = 2 endif @@ -1784,66 +1793,89 @@ integer function qmckl_compute_factor_ee_f(context, walk_num, elec_num, up_num, end function qmckl_compute_factor_ee_f #+end_src - #+CALL: generate_c_header(table=qmckl_factor_ee_args,rettyp=get_value("CRetType"),fname=get_value("Name")) +#+begin_src c :comments org :tangle (eval c) :noweb yes +qmckl_exit_code qmckl_compute_factor_ee ( + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t up_num, + const int64_t bord_num, + const double* bord_vector, + const double* ee_distance_rescaled, + const double* asymp_jasb, + double* const factor_ee ) { - #+RESULTS: - #+begin_src c :tangle (eval h_func) :comments org + int ipar; // can we use a smaller integer? + double x, x1, spin_fact, power_ser; + + if (context == QMCKL_NULL_CONTEXT) { + return QMCKL_INVALID_CONTEXT; + } + + if (walk_num <= 0) { + return QMCKL_INVALID_ARG_2; + } + + if (elec_num <= 0) { + return QMCKL_INVALID_ARG_3; + } + + if (bord_num <= 0) { + return QMCKL_INVALID_ARG_4; + } + + for (int nw = 0; nw < walk_num; ++nw) { + factor_ee[nw] = 0.0; // put init array here. + for (int i = 0; i < elec_num; ++i ) { + for (int j = 0; j < i; ++j) { + //x = ee_distance_rescaled[j * (walk_num * elec_num) + i * (walk_num) + nw]; + x = ee_distance_rescaled[j + i * elec_num + nw*(elec_num * elec_num)]; + x1 = x; + power_ser = 0.0; + spin_fact = 1.0; + ipar = 0; // index of asymp_jasb + + for (int p = 1; p < bord_num; ++p) { + x = x * x1; + power_ser = power_ser + bord_vector[p + 1] * x; + } + + if(i < up_num || j >= up_num) { + spin_fact = 0.5; + ipar = 1; + } + + factor_ee[nw] = factor_ee[nw] + spin_fact * bord_vector[0] * \ + x1 / \ + (1.0 + bord_vector[1] * \ + x1) \ + -asymp_jasb[ipar] + power_ser; + + } + } + } + + return QMCKL_SUCCESS; +} +#+end_src + +# #+CALL: generate_c_header(table=qmckl_factor_ee_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_factor_ee ( - const qmckl_context context, - const int64_t walk_num, - const int64_t elec_num, - const int64_t up_num, - const int64_t bord_num, - const double* bord_vector, - const double* ee_distance_rescaled, - const double* asymp_jasb, - double* const factor_ee ); + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t up_num, + const int64_t bord_num, + const double* bord_vector, + const double* ee_distance_rescaled, + const double* asymp_jasb, + double* const factor_ee ); #+end_src - #+CALL: generate_c_interface(table=qmckl_factor_ee_args,rettyp=get_value("CRetType"),fname=get_value("Name")) - #+RESULTS: - #+begin_src f90 :tangle (eval f) :comments org :exports none - integer(c_int32_t) function qmckl_compute_factor_ee & - (context, & - walk_num, & - elec_num, & - up_num, & - bord_num, & - bord_vector, & - ee_distance_rescaled, & - asymp_jasb, & - factor_ee) & - bind(C) result(info) - - use, intrinsic :: iso_c_binding - implicit none - - integer (c_int64_t) , intent(in) , value :: context - integer (c_int64_t) , intent(in) , value :: walk_num - integer (c_int64_t) , intent(in) , value :: elec_num - integer (c_int64_t) , intent(in) , value :: up_num - integer (c_int64_t) , intent(in) , value :: bord_num - real (c_double ) , intent(in) :: bord_vector(bord_num + 1) - real (c_double ) , intent(in) :: ee_distance_rescaled(elec_num,elec_num,walk_num) - real (c_double ) , intent(in) :: asymp_jasb(2) - real (c_double ) , intent(out) :: factor_ee(walk_num) - - integer(c_int32_t), external :: qmckl_compute_factor_ee_f - info = qmckl_compute_factor_ee_f & - (context, & - walk_num, & - elec_num, & - up_num, & - bord_num, & - bord_vector, & - ee_distance_rescaled, & - asymp_jasb, & - factor_ee) - - end function qmckl_compute_factor_ee - #+end_src *** Test #+begin_src python :results output :exports none :noweb yes @@ -1928,7 +1960,7 @@ qmckl_get_jastrow_factor_ee_deriv_e(qmckl_context context, rc = qmckl_provide_factor_ee_deriv_e(context); if (rc != QMCKL_SUCCESS) return rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int64_t sze = ctx->electron.walk_num * 4 * ctx->electron.num; @@ -1960,7 +1992,7 @@ qmckl_exit_code qmckl_provide_factor_ee_deriv_e(qmckl_context context) return QMCKL_NULL_CONTEXT; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); /* Check if ee rescaled distance is provided */ @@ -1990,17 +2022,15 @@ qmckl_exit_code qmckl_provide_factor_ee_deriv_e(qmckl_context context) ctx->jastrow.factor_ee_deriv_e = factor_ee_deriv_e; } - qmckl_exit_code rc = - qmckl_compute_factor_ee_deriv_e(context, - ctx->electron.walk_num, - ctx->electron.num, - ctx->electron.up_num, - ctx->jastrow.bord_num, - ctx->jastrow.bord_vector, - ctx->electron.ee_distance_rescaled, - ctx->electron.ee_distance_rescaled_deriv_e, - ctx->jastrow.asymp_jasb, - ctx->jastrow.factor_ee_deriv_e); + rc = qmckl_compute_factor_ee_deriv_e(context, + ctx->electron.walk_num, + ctx->electron.num, + ctx->electron.up_num, + ctx->jastrow.bord_num, + ctx->jastrow.bord_vector, + ctx->electron.ee_distance_rescaled, + ctx->electron.ee_distance_rescaled_deriv_e, + ctx->jastrow.factor_ee_deriv_e); if (rc != QMCKL_SUCCESS) { return rc; } @@ -2030,13 +2060,13 @@ qmckl_exit_code qmckl_provide_factor_ee_deriv_e(qmckl_context context) | ~bord_vector~ | ~double[bord_num+1]~ | in | List of coefficients | | ~ee_distance_rescaled~ | ~double[walk_num][elec_num][elec_num]~ | in | Electron-electron distances | | ~ee_distance_rescaled_deriv_e~ | ~double[walk_num][4][elec_num][elec_num]~ | in | Electron-electron distances | - | ~asymp_jasb~ | ~double[2]~ | in | Electron-electron distances | | ~factor_ee_deriv_e~ | ~double[walk_num][4][elec_num]~ | out | Electron-electron distances | #+begin_src f90 :comments org :tangle (eval f) :noweb yes -integer function qmckl_compute_factor_ee_deriv_e_f(context, walk_num, elec_num, up_num, bord_num, & - bord_vector, ee_distance_rescaled, ee_distance_rescaled_deriv_e, & - asymp_jasb, factor_ee_deriv_e) & +integer function qmckl_compute_factor_ee_deriv_e_doc_f( & + context, walk_num, elec_num, up_num, bord_num, & + bord_vector, ee_distance_rescaled, ee_distance_rescaled_deriv_e, & + factor_ee_deriv_e) & result(info) use qmckl implicit none @@ -2044,11 +2074,10 @@ integer function qmckl_compute_factor_ee_deriv_e_f(context, walk_num, elec_num, integer*8 , intent(in) :: walk_num, elec_num, bord_num, up_num double precision , intent(in) :: bord_vector(bord_num + 1) double precision , intent(in) :: ee_distance_rescaled(elec_num, elec_num,walk_num) - double precision , intent(in) :: ee_distance_rescaled_deriv_e(4,elec_num, elec_num,walk_num) - double precision , intent(in) :: asymp_jasb(2) + double precision , intent(in) :: ee_distance_rescaled_deriv_e(4,elec_num, elec_num,walk_num) !TODO double precision , intent(out) :: factor_ee_deriv_e(elec_num,4,walk_num) - integer*8 :: i, j, p, ipar, nw, ii + integer*8 :: i, j, p, nw, ii double precision :: x, spin_fact, y double precision :: den, invden, invden2, invden3, xinv double precision :: lap1, lap2, lap3, third @@ -2092,7 +2121,6 @@ integer function qmckl_compute_factor_ee_deriv_e_f(context, walk_num, elec_num, invden2 = invden * invden invden3 = invden2 * invden xinv = 1.0d0 / (x + 1.0d-18) - ipar = 1 dx(1) = ee_distance_rescaled_deriv_e(1, i, j, nw) dx(2) = ee_distance_rescaled_deriv_e(2, i, j, nw) @@ -2134,43 +2162,154 @@ integer function qmckl_compute_factor_ee_deriv_e_f(context, walk_num, elec_num, end do end do -end function qmckl_compute_factor_ee_deriv_e_f +end function qmckl_compute_factor_ee_deriv_e_doc_f #+end_src - #+CALL: generate_c_header(table=qmckl_factor_ee_deriv_e_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + #+begin_src c :comments org :tangle (eval c) :noweb yes +qmckl_exit_code qmckl_compute_factor_ee_deriv_e_hpc( + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t up_num, + const int64_t bord_num, + const double* bord_vector, + const double* ee_distance_rescaled, + const double* ee_distance_rescaled_deriv_e, + double* const factor_ee_deriv_e ) { - #+RESULTS: - #+begin_src c :tangle (eval h_func) :comments org + int64_t ii; + double pow_ser_g[3]; + double dx[4]; + double x, spin_fact, y; + double den, invden, invden2, invden3, xinv; + double lap1, lap2, lap3, third; + + if (context == QMCKL_NULL_CONTEXT) { + return QMCKL_INVALID_CONTEXT; + } + + if (walk_num <= 0) { + return QMCKL_INVALID_ARG_2; + } + + if (elec_num <= 0) { + return QMCKL_INVALID_ARG_3; + } + + if (bord_num <= 0) { + return QMCKL_INVALID_ARG_4; + } + + + for (int nw = 0; nw < walk_num; ++nw) { + for (int ii = 0; ii < 4; ++ii) { + for (int j = 0; j < elec_num; ++j) { + factor_ee_deriv_e[j + ii * elec_num + nw * elec_num * 4] = 0.0; + } + } + } + + third = 1.0 / 3.0; + + for (int nw = 0; nw < walk_num; ++nw) { + for (int i = 0; i < elec_num; ++i) { + for (int j = 0; j < elec_num; ++j) { + x = ee_distance_rescaled[j + i * elec_num + nw * elec_num * elec_num]; + if (fabs(x) < 1.0e-18) continue; + for (int ii = 0; ii < 3; ++ii){ + pow_ser_g[ii] = 0.0; + } + spin_fact = 1.0; + den = 1.0 + bord_vector[1] * x; + invden = 1.0 / den; + invden2 = invden * invden; + invden3 = invden2 * invden; + xinv = 1.0 / (x + 1.0e-18); + + dx[0] = ee_distance_rescaled_deriv_e[0 \ + + j * 4 + i * 4 * elec_num \ + + nw * 4 * elec_num * elec_num]; + dx[1] = ee_distance_rescaled_deriv_e[1 \ + + j * 4 + i * 4 * elec_num \ + + nw * 4 * elec_num * elec_num]; + dx[2] = ee_distance_rescaled_deriv_e[2 \ + + j * 4 + i * 4 * elec_num \ + + nw * 4 * elec_num * elec_num]; + dx[3] = ee_distance_rescaled_deriv_e[3 \ + + j * 4 + i * 4 * elec_num \ + + nw * 4 * elec_num * elec_num]; + + if((i <= (up_num-1) && j <= (up_num-1) ) || (i > (up_num-1) && j > (up_num-1))) { + spin_fact = 0.5; + } + + lap1 = 0.0; + lap2 = 0.0; + lap3 = 0.0; + for (int ii = 0; ii < 3; ++ii) { + x = ee_distance_rescaled[j + i * elec_num + nw * elec_num * elec_num]; + if (fabs(x) < 1.0e-18) continue; + for (int p = 2; p < bord_num+1; ++p) { + y = p * bord_vector[(p-1) + 1] * x; + pow_ser_g[ii] = pow_ser_g[ii] + y * dx[ii]; + lap1 = lap1 + (p - 1) * y * xinv * dx[ii] * dx[ii]; + lap2 = lap2 + y; + x = x * ee_distance_rescaled[j + i * elec_num + nw * elec_num * elec_num]; + } + + lap3 = lap3 - 2.0 * bord_vector[1] * dx[ii] * dx[ii]; + + factor_ee_deriv_e[i + ii * elec_num + nw * elec_num * 4 ] += \ + + spin_fact * bord_vector[0] * dx[ii] * invden2 \ + + pow_ser_g[ii] ; + } + + ii = 3; + lap2 = lap2 * dx[ii] * third; + lap3 = lap3 + den * dx[ii]; + lap3 = lap3 * (spin_fact * bord_vector[0] * invden3); + factor_ee_deriv_e[i + ii*elec_num + nw * elec_num * 4] += lap1 + lap2 + lap3; + + } + } + } + + return QMCKL_SUCCESS; +} + #+end_src + +# #+CALL: generate_c_header(table=qmckl_factor_ee_deriv_e_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + + + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_factor_ee_deriv_e ( - const qmckl_context context, - const int64_t walk_num, - const int64_t elec_num, - const int64_t up_num, - const int64_t bord_num, - const double* bord_vector, - const double* ee_distance_rescaled, - const double* ee_distance_rescaled_deriv_e, - const double* asymp_jasb, - double* const factor_ee_deriv_e ); + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t up_num, + const int64_t bord_num, + const double* bord_vector, + const double* ee_distance_rescaled, + const double* ee_distance_rescaled_deriv_e, + double* const factor_ee_deriv_e ); #+end_src - #+CALL: generate_c_interface(table=qmckl_factor_ee_deriv_e_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + #+CALL: generate_c_interface(table=qmckl_factor_ee_deriv_e_args,rettyp=get_value("CRetType"),fname="qmckl_compute_factor_ee_deriv_e_doc") #+RESULTS: #+begin_src f90 :tangle (eval f) :comments org :exports none - integer(c_int32_t) function qmckl_compute_factor_ee_deriv_e & - (context, & - walk_num, & - elec_num, & - up_num, & - bord_num, & - bord_vector, & - ee_distance_rescaled, & - ee_distance_rescaled_deriv_e, & - asymp_jasb, & - factor_ee_deriv_e) & - bind(C) result(info) +integer(c_int32_t) function qmckl_compute_factor_ee_deriv_e_doc & + (context, & + walk_num, & + elec_num, & + up_num, & + bord_num, & + bord_vector, & + ee_distance_rescaled, & + ee_distance_rescaled_deriv_e, & + factor_ee_deriv_e) & + bind(C) result(info) use, intrinsic :: iso_c_binding implicit none @@ -2180,27 +2319,77 @@ end function qmckl_compute_factor_ee_deriv_e_f integer (c_int64_t) , intent(in) , value :: elec_num integer (c_int64_t) , intent(in) , value :: up_num integer (c_int64_t) , intent(in) , value :: bord_num - real (c_double ) , intent(in) :: bord_vector(bord_num + 1) + real (c_double ) , intent(in) :: bord_vector(bord_num+1) real (c_double ) , intent(in) :: ee_distance_rescaled(elec_num,elec_num,walk_num) real (c_double ) , intent(in) :: ee_distance_rescaled_deriv_e(elec_num,elec_num,4,walk_num) - real (c_double ) , intent(in) :: asymp_jasb(2) real (c_double ) , intent(out) :: factor_ee_deriv_e(elec_num,4,walk_num) - integer(c_int32_t), external :: qmckl_compute_factor_ee_deriv_e_f - info = qmckl_compute_factor_ee_deriv_e_f & - (context, & - walk_num, & - elec_num, & - up_num, & - bord_num, & - bord_vector, & - ee_distance_rescaled, & - ee_distance_rescaled_deriv_e, & - asymp_jasb, & - factor_ee_deriv_e) + integer(c_int32_t), external :: qmckl_compute_factor_ee_deriv_e_doc_f + info = qmckl_compute_factor_ee_deriv_e_doc_f & + (context, & + walk_num, & + elec_num, & + up_num, & + bord_num, & + bord_vector, & + ee_distance_rescaled, & + ee_distance_rescaled_deriv_e, & + factor_ee_deriv_e) - end function qmckl_compute_factor_ee_deriv_e + end function qmckl_compute_factor_ee_deriv_e_doc #+end_src + + #+begin_src c :tangle (eval h_private_func) :comments org + qmckl_exit_code qmckl_compute_factor_ee_deriv_e_hpc ( + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t up_num, + const int64_t bord_num, + const double* bord_vector, + const double* ee_distance_rescaled, + const double* ee_distance_rescaled_deriv_e, + double* const factor_ee_deriv_e ); + + #+end_src + + #+begin_src c :tangle (eval h_private_func) :comments org + qmckl_exit_code qmckl_compute_factor_ee_deriv_e_doc ( + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t up_num, + const int64_t bord_num, + const double* bord_vector, + const double* ee_distance_rescaled, + const double* ee_distance_rescaled_deriv_e, + double* const factor_ee_deriv_e ); + #+end_src + + + #+begin_src c :comments org :tangle (eval c) :noweb yes + qmckl_exit_code qmckl_compute_factor_ee_deriv_e ( + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t up_num, + const int64_t bord_num, + const double* bord_vector, + const double* ee_distance_rescaled, + const double* ee_distance_rescaled_deriv_e, + double* const factor_ee_deriv_e ) { + + #ifdef HAVE_HPC + return qmckl_compute_factor_ee_deriv_e_hpc(context, walk_num, elec_num, up_num, bord_num, bord_vector, ee_distance_rescaled, ee_distance_rescaled_deriv_e, factor_ee_deriv_e ); + #else + return qmckl_compute_factor_ee_deriv_e_doc(context, walk_num, elec_num, up_num, bord_num, bord_vector, ee_distance_rescaled, ee_distance_rescaled_deriv_e, factor_ee_deriv_e ); + #endif +} + #+end_src + + + + *** Test #+begin_src python :results output :exports none :noweb yes @@ -2319,7 +2508,6 @@ assert(fabs(factor_ee_deriv_e[0][0][0]-0.16364894652107934) < 1.e-12); assert(fabs(factor_ee_deriv_e[0][1][0]+0.6927548119830084 ) < 1.e-12); assert(fabs(factor_ee_deriv_e[0][2][0]-0.073267755223968 ) < 1.e-12); assert(fabs(factor_ee_deriv_e[0][3][0]-1.5111672803213185 ) < 1.e-12); - #+end_src ** Electron-nucleus component \(f_{en}\) @@ -2355,7 +2543,7 @@ qmckl_get_jastrow_factor_en(qmckl_context context, rc = qmckl_provide_factor_en(context); if (rc != QMCKL_SUCCESS) return rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int64_t sze=ctx->electron.walk_num; @@ -2386,7 +2574,7 @@ qmckl_exit_code qmckl_provide_factor_en(qmckl_context context) return QMCKL_NULL_CONTEXT; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); /* Check if en rescaled distance is provided */ @@ -2412,21 +2600,20 @@ qmckl_exit_code qmckl_provide_factor_en(qmckl_context context) ctx->jastrow.factor_en = factor_en; } - qmckl_exit_code rc = - qmckl_compute_factor_en(context, - ctx->electron.walk_num, - ctx->electron.num, - ctx->nucleus.num, - ctx->jastrow.type_nucl_num, - ctx->jastrow.type_nucl_vector, - ctx->jastrow.aord_num, - ctx->jastrow.aord_vector, - ctx->electron.en_distance_rescaled, - ctx->jastrow.factor_en); + rc = qmckl_compute_factor_en(context, + ctx->electron.walk_num, + ctx->electron.num, + ctx->nucleus.num, + ctx->jastrow.type_nucl_num, + ctx->jastrow.type_nucl_vector, + ctx->jastrow.aord_num, + ctx->jastrow.aord_vector, + ctx->electron.en_distance_rescaled, + ctx->jastrow.factor_en); if (rc != QMCKL_SUCCESS) { return rc; } - + ctx->jastrow.factor_en_date = ctx->date; } @@ -2456,9 +2643,10 @@ qmckl_exit_code qmckl_provide_factor_en(qmckl_context context) | ~factor_en~ | ~double[walk_num]~ | out | Electron-nucleus jastrow | #+begin_src f90 :comments org :tangle (eval f) :noweb yes -integer function qmckl_compute_factor_en_f(context, walk_num, elec_num, nucl_num, type_nucl_num, & - type_nucl_vector, aord_num, aord_vector, & - en_distance_rescaled, factor_en) & +integer function qmckl_compute_factor_en_f( & + context, walk_num, elec_num, nucl_num, type_nucl_num, & + type_nucl_vector, aord_num, aord_vector, & + en_distance_rescaled, factor_en) & result(info) use qmckl implicit none @@ -2469,8 +2657,8 @@ integer function qmckl_compute_factor_en_f(context, walk_num, elec_num, nucl_num double precision , intent(in) :: en_distance_rescaled(elec_num, nucl_num, walk_num) double precision , intent(out) :: factor_en(walk_num) - integer*8 :: i, a, p, ipar, nw - double precision :: x, spin_fact, power_ser + integer*8 :: i, a, p, nw + double precision :: x, power_ser info = QMCKL_SUCCESS @@ -2525,71 +2713,109 @@ integer function qmckl_compute_factor_en_f(context, walk_num, elec_num, nucl_num end function qmckl_compute_factor_en_f #+end_src - #+CALL: generate_c_header(table=qmckl_factor_en_args,rettyp=get_value("CRetType"),fname=get_value("Name")) - #+RESULTS: - #+begin_src c :tangle (eval h_func) :comments org + #+begin_src c :comments org :tangle (eval c) :noweb yes +qmckl_exit_code qmckl_compute_factor_en ( + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t nucl_num, + const int64_t type_nucl_num, + const int64_t* type_nucl_vector, + const int64_t aord_num, + const double* aord_vector, + const double* en_distance_rescaled, + double* const factor_en ) { + + double x, x1, power_ser; + + + if (context == QMCKL_NULL_CONTEXT) { + return QMCKL_INVALID_CONTEXT; + } + + if (walk_num <= 0) { + return QMCKL_INVALID_ARG_2; + } + + if (elec_num <= 0) { + return QMCKL_INVALID_ARG_3; + } + + if (nucl_num <= 0) { + return QMCKL_INVALID_ARG_4; + } + + if (type_nucl_num <= 0) { + return QMCKL_INVALID_ARG_5; + } + + if (type_nucl_vector == NULL) { + return QMCKL_INVALID_ARG_6; + } + + if (aord_num <= 0) { + return QMCKL_INVALID_ARG_7; + } + + if (aord_vector == NULL) { + return QMCKL_INVALID_ARG_8; + } + + if (en_distance_rescaled == NULL) { + return QMCKL_INVALID_ARG_9; + } + + if (factor_en == NULL) { + return QMCKL_INVALID_ARG_10; + } + + + for (int nw = 0; nw < walk_num; ++nw ) { + // init array + factor_en[nw] = 0.0; + for (int a = 0; a < nucl_num; ++a ) { + for (int i = 0; i < elec_num; ++i ) { + // x = ee_distance_rescaled[j * (walk_num * elec_num) + i * (walk_num) + nw]; + x = en_distance_rescaled[i + a * elec_num + nw * (elec_num * nucl_num)]; + x1 = x; + power_ser = 0.0; + + for (int p = 2; p < aord_num+1; ++p) { + x = x * x1; + power_ser = power_ser + aord_vector[(p+1)-1 + (type_nucl_vector[a]-1) * aord_num] * x; + } + + factor_en[nw] = factor_en[nw] + aord_vector[0 + (type_nucl_vector[a]-1)*aord_num] * x1 / \ + (1.0 + aord_vector[1 + (type_nucl_vector[a]-1) * aord_num] * x1) + \ + power_ser; + + } + } + } + + return QMCKL_SUCCESS; +} + #+end_src + + +# #+CALL: generate_c_header(table=qmckl_factor_en_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_factor_en ( - const qmckl_context context, - const int64_t walk_num, - const int64_t elec_num, - const int64_t nucl_num, - const int64_t type_nucl_num, - const int64_t* type_nucl_vector, - const int64_t aord_num, - const double* aord_vector, - const double* en_distance_rescaled, - double* const factor_en ); + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t nucl_num, + const int64_t type_nucl_num, + const int64_t* type_nucl_vector, + const int64_t aord_num, + const double* aord_vector, + const double* en_distance_rescaled, + double* const factor_en ); #+end_src - #+CALL: generate_c_interface(table=qmckl_factor_en_args,rettyp=get_value("CRetType"),fname=get_value("Name")) - - #+RESULTS: - #+begin_src f90 :tangle (eval f) :comments org :exports none - integer(c_int32_t) function qmckl_compute_factor_en & - (context, & - walk_num, & - elec_num, & - nucl_num, & - type_nucl_num, & - type_nucl_vector, & - aord_num, & - aord_vector, & - en_distance_rescaled, & - factor_en) & - bind(C) result(info) - - use, intrinsic :: iso_c_binding - implicit none - - integer (c_int64_t) , intent(in) , value :: context - integer (c_int64_t) , intent(in) , value :: walk_num - integer (c_int64_t) , intent(in) , value :: elec_num - integer (c_int64_t) , intent(in) , value :: nucl_num - integer (c_int64_t) , intent(in) , value :: type_nucl_num - integer (c_int64_t) , intent(in) :: type_nucl_vector(nucl_num) - integer (c_int64_t) , intent(in) , value :: aord_num - real (c_double ) , intent(in) :: aord_vector(aord_num + 1, type_nucl_num) - real (c_double ) , intent(in) :: en_distance_rescaled(elec_num, nucl_num, walk_num) - real (c_double ) , intent(out) :: factor_en(walk_num) - - integer(c_int32_t), external :: qmckl_compute_factor_en_f - info = qmckl_compute_factor_en_f & - (context, & - walk_num, & - elec_num, & - nucl_num, & - type_nucl_num, & - type_nucl_vector, & - aord_num, & - aord_vector, & - en_distance_rescaled, & - factor_en) - - end function qmckl_compute_factor_en - #+end_src - *** Test #+begin_src python :results output :exports none :noweb yes import numpy as np @@ -2659,7 +2885,7 @@ qmckl_get_jastrow_factor_en_deriv_e(qmckl_context context, rc = qmckl_provide_factor_en_deriv_e(context); if (rc != QMCKL_SUCCESS) return rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int64_t sze = ctx->electron.walk_num * 4 * ctx->electron.num; @@ -2690,7 +2916,7 @@ qmckl_exit_code qmckl_provide_factor_en_deriv_e(qmckl_context context) return QMCKL_NULL_CONTEXT; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); /* Check if en rescaled distance is provided */ @@ -2720,18 +2946,17 @@ qmckl_exit_code qmckl_provide_factor_en_deriv_e(qmckl_context context) ctx->jastrow.factor_en_deriv_e = factor_en_deriv_e; } - qmckl_exit_code rc = - qmckl_compute_factor_en_deriv_e(context, - ctx->electron.walk_num, - ctx->electron.num, - ctx->nucleus.num, - ctx->jastrow.type_nucl_num, - ctx->jastrow.type_nucl_vector, - ctx->jastrow.aord_num, - ctx->jastrow.aord_vector, - ctx->electron.en_distance_rescaled, - ctx->electron.en_distance_rescaled_deriv_e, - ctx->jastrow.factor_en_deriv_e); + rc = qmckl_compute_factor_en_deriv_e(context, + ctx->electron.walk_num, + ctx->electron.num, + ctx->nucleus.num, + ctx->jastrow.type_nucl_num, + ctx->jastrow.type_nucl_vector, + ctx->jastrow.aord_num, + ctx->jastrow.aord_vector, + ctx->electron.en_distance_rescaled, + ctx->electron.en_distance_rescaled_deriv_e, + ctx->jastrow.factor_en_deriv_e); if (rc != QMCKL_SUCCESS) { return rc; } @@ -2766,9 +2991,10 @@ qmckl_exit_code qmckl_provide_factor_en_deriv_e(qmckl_context context) | ~factor_en_deriv_e~ | ~double[walk_num][4][elec_num]~ | out | Electron-nucleus jastrow | #+begin_src f90 :comments org :tangle (eval f) :noweb yes -integer function qmckl_compute_factor_en_deriv_e_f(context, walk_num, elec_num, nucl_num, type_nucl_num, & - type_nucl_vector, aord_num, aord_vector, & - en_distance_rescaled, en_distance_rescaled_deriv_e, factor_en_deriv_e) & +integer function qmckl_compute_factor_en_deriv_e_f( & + context, walk_num, elec_num, nucl_num, type_nucl_num, & + type_nucl_vector, aord_num, aord_vector, & + en_distance_rescaled, en_distance_rescaled_deriv_e, factor_en_deriv_e) & result(info) use qmckl implicit none @@ -2781,7 +3007,7 @@ integer function qmckl_compute_factor_en_deriv_e_f(context, walk_num, elec_num, double precision , intent(out) :: factor_en_deriv_e(elec_num,4,walk_num) integer*8 :: i, a, p, ipar, nw, ii - double precision :: x, spin_fact, den, invden, invden2, invden3, xinv + double precision :: x, den, invden, invden2, invden3, xinv double precision :: y, lap1, lap2, lap3, third double precision, dimension(3) :: power_ser_g double precision, dimension(4) :: dx @@ -2848,7 +3074,7 @@ integer function qmckl_compute_factor_en_deriv_e_f(context, walk_num, elec_num, lap3 = lap3 - 2.0d0 * aord_vector(2, type_nucl_vector(a)) * dx(ii) * dx(ii) factor_en_deriv_e(i, ii, nw) = factor_en_deriv_e(i, ii, nw) + aord_vector(1, type_nucl_vector(a)) & - * dx(ii) * invden2 & + ,* dx(ii) * invden2 & + power_ser_g(ii) end do @@ -2866,22 +3092,21 @@ integer function qmckl_compute_factor_en_deriv_e_f(context, walk_num, elec_num, end function qmckl_compute_factor_en_deriv_e_f #+end_src - #+CALL: generate_c_header(table=qmckl_factor_en_deriv_e_args,rettyp=get_value("CRetType"),fname=get_value("Name")) +# #+CALL: generate_c_header(table=qmckl_factor_en_deriv_e_args,rettyp=get_value("CRetType"),fname=get_value("Name")) - #+RESULTS: - #+begin_src c :tangle (eval h_func) :comments org + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_factor_en_deriv_e ( - const qmckl_context context, - const int64_t walk_num, - const int64_t elec_num, - const int64_t nucl_num, - const int64_t type_nucl_num, - const int64_t* type_nucl_vector, - const int64_t aord_num, - const double* aord_vector, - const double* en_distance_rescaled, - const double* en_distance_rescaled_deriv_e, - double* const factor_en_deriv_e ); + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t nucl_num, + const int64_t type_nucl_num, + const int64_t* type_nucl_vector, + const int64_t aord_num, + const double* aord_vector, + const double* en_distance_rescaled, + const double* en_distance_rescaled_deriv_e, + double* const factor_en_deriv_e ); #+end_src @@ -2890,18 +3115,18 @@ end function qmckl_compute_factor_en_deriv_e_f #+RESULTS: #+begin_src f90 :tangle (eval f) :comments org :exports none integer(c_int32_t) function qmckl_compute_factor_en_deriv_e & - (context, & - walk_num, & - elec_num, & - nucl_num, & - type_nucl_num, & - type_nucl_vector, & - aord_num, & - aord_vector, & - en_distance_rescaled, & - en_distance_rescaled_deriv_e, & - factor_en_deriv_e) & - bind(C) result(info) + (context, & + walk_num, & + elec_num, & + nucl_num, & + type_nucl_num, & + type_nucl_vector, & + aord_num, & + aord_vector, & + en_distance_rescaled, & + en_distance_rescaled_deriv_e, & + factor_en_deriv_e) & + bind(C) result(info) use, intrinsic :: iso_c_binding implicit none @@ -2913,24 +3138,24 @@ end function qmckl_compute_factor_en_deriv_e_f integer (c_int64_t) , intent(in) , value :: type_nucl_num integer (c_int64_t) , intent(in) :: type_nucl_vector(nucl_num) integer (c_int64_t) , intent(in) , value :: aord_num - real (c_double ) , intent(in) :: aord_vector(aord_num + 1, type_nucl_num) + real (c_double ) , intent(in) :: aord_vector(type_nucl_num,aord_num+1) real (c_double ) , intent(in) :: en_distance_rescaled(elec_num,nucl_num,walk_num) real (c_double ) , intent(in) :: en_distance_rescaled_deriv_e(elec_num,nucl_num,4,walk_num) real (c_double ) , intent(out) :: factor_en_deriv_e(elec_num,4,walk_num) integer(c_int32_t), external :: qmckl_compute_factor_en_deriv_e_f info = qmckl_compute_factor_en_deriv_e_f & - (context, & - walk_num, & - elec_num, & - nucl_num, & - type_nucl_num, & - type_nucl_vector, & - aord_num, & - aord_vector, & - en_distance_rescaled, & - en_distance_rescaled_deriv_e, & - factor_en_deriv_e) + (context, & + walk_num, & + elec_num, & + nucl_num, & + type_nucl_num, & + type_nucl_vector, & + aord_num, & + aord_vector, & + en_distance_rescaled, & + en_distance_rescaled_deriv_e, & + factor_en_deriv_e) end function qmckl_compute_factor_en_deriv_e #+end_src @@ -3081,7 +3306,7 @@ qmckl_get_jastrow_een_rescaled_e(qmckl_context context, rc = qmckl_provide_een_rescaled_e(context); if (rc != QMCKL_SUCCESS) return rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int64_t sze = ctx->electron.num * ctx->electron.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1); @@ -3111,7 +3336,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_e(qmckl_context context) return QMCKL_NULL_CONTEXT; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); /* Check if ee distance is provided */ @@ -3138,14 +3363,13 @@ qmckl_exit_code qmckl_provide_een_rescaled_e(qmckl_context context) ctx->jastrow.een_rescaled_e = een_rescaled_e; } - qmckl_exit_code rc = - qmckl_compute_een_rescaled_e(context, - ctx->electron.walk_num, - ctx->electron.num, - ctx->jastrow.cord_num, - ctx->electron.rescale_factor_kappa_ee, - ctx->electron.ee_distance, - ctx->jastrow.een_rescaled_e); + rc = qmckl_compute_een_rescaled_e(context, + ctx->electron.walk_num, + ctx->electron.num, + ctx->jastrow.cord_num, + ctx->electron.rescale_factor_kappa_ee, + ctx->electron.ee_distance, + ctx->jastrow.een_rescaled_e); if (rc != QMCKL_SUCCESS) { return rc; } @@ -3176,7 +3400,8 @@ qmckl_exit_code qmckl_provide_een_rescaled_e(qmckl_context context) | ~een_rescaled_e~ | ~double[walk_num][0:cord_num][elec_num][elec_num]~ | out | Electron-electron rescaled distances | #+begin_src f90 :comments org :tangle (eval f) :noweb yes -integer function qmckl_compute_een_rescaled_e_f(context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, & +integer function qmckl_compute_een_rescaled_e_doc_f( & + context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, & ee_distance, een_rescaled_e) & result(info) use qmckl @@ -3194,7 +3419,6 @@ integer function qmckl_compute_een_rescaled_e_f(context, walk_num, elec_num, cor allocate(een_rescaled_e_ij(elec_num * (elec_num - 1) / 2, cord_num + 1)) - info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then @@ -3223,6 +3447,7 @@ integer function qmckl_compute_een_rescaled_e_f(context, walk_num, elec_num, cor een_rescaled_e_ij = 0.0d0 een_rescaled_e_ij(:, 1) = 1.0d0 + k = 0 do j = 1, elec_num do i = 1, j - 1 @@ -3231,6 +3456,7 @@ integer function qmckl_compute_een_rescaled_e_f(context, walk_num, elec_num, cor end do end do + do l = 2, cord_num do k = 1, elec_num * (elec_num - 1)/2 een_rescaled_e_ij(k, l + 1) = een_rescaled_e_ij(k, l + 1 - 1) * een_rescaled_e_ij(k, 2) @@ -3239,6 +3465,7 @@ integer function qmckl_compute_een_rescaled_e_f(context, walk_num, elec_num, cor ! prepare the actual een table een_rescaled_e(:, :, 0, nw) = 1.0d0 + do l = 1, cord_num k = 0 do j = 1, elec_num @@ -3259,30 +3486,30 @@ integer function qmckl_compute_een_rescaled_e_f(context, walk_num, elec_num, cor end do -end function qmckl_compute_een_rescaled_e_f +end function qmckl_compute_een_rescaled_e_doc_f #+end_src - #+CALL: generate_c_header(table=qmckl_factor_een_rescaled_e_args,rettyp=get_value("CRetType"),fname=get_value("Name")) +# #+CALL: generate_c_header(table=qmckl_factor_een_rescaled_e_args,rettyp=get_value("CRetType"),fname=get_value("Name")) - #+RESULTS: - #+begin_src c :tangle (eval h_func) :comments org + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_een_rescaled_e ( - const qmckl_context context, - const int64_t walk_num, - const int64_t elec_num, - const int64_t cord_num, - const double rescale_factor_kappa_ee, - const double* ee_distance, - double* const een_rescaled_e ); + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t cord_num, + const double rescale_factor_kappa_ee, + const double* ee_distance, + double* const een_rescaled_e ); #+end_src - #+CALL: generate_c_interface(table=qmckl_factor_een_rescaled_e_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + #+CALL: generate_c_interface(table=qmckl_factor_een_rescaled_e_args,rettyp=get_value("CRetType"),fname="qmckl_compute_een_rescaled_e_doc") #+RESULTS: #+begin_src f90 :tangle (eval f) :comments org :exports none - integer(c_int32_t) function qmckl_compute_een_rescaled_e & - (context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, ee_distance, een_rescaled_e) & - bind(C) result(info) + integer(c_int32_t) function qmckl_compute_een_rescaled_e_doc & + (context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, & + ee_distance, een_rescaled_e) & + bind(C) result(info) use, intrinsic :: iso_c_binding implicit none @@ -3295,13 +3522,186 @@ end function qmckl_compute_een_rescaled_e_f real (c_double ) , intent(in) :: ee_distance(elec_num,elec_num,walk_num) real (c_double ) , intent(out) :: een_rescaled_e(elec_num,elec_num,0:cord_num,walk_num) - integer(c_int32_t), external :: qmckl_compute_een_rescaled_e_f - info = qmckl_compute_een_rescaled_e_f & + integer(c_int32_t), external :: qmckl_compute_een_rescaled_e_doc_f + info = qmckl_compute_een_rescaled_e_doc_f & (context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, ee_distance, een_rescaled_e) - end function qmckl_compute_een_rescaled_e + end function qmckl_compute_een_rescaled_e_doc #+end_src + #+begin_src c :comments org :tangle (eval c) :noweb yes +qmckl_exit_code qmckl_compute_een_rescaled_e_hpc ( + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t cord_num, + const double rescale_factor_kappa_ee, + const double* ee_distance, + double* const een_rescaled_e ) { + + double *een_rescaled_e_ij; + double x; + const int64_t elec_pairs = (elec_num * (elec_num - 1)) / 2; + const int64_t len_een_ij = elec_pairs * (cord_num + 1); + int64_t k; + + // number of element for the een_rescaled_e_ij[N_e*(N_e-1)/2][cord+1] + // probably in C is better [cord+1, Ne*(Ne-1)/2] + //elec_pairs = (elec_num * (elec_num - 1)) / 2; + //len_een_ij = elec_pairs * (cord_num + 1); + een_rescaled_e_ij = (double *) malloc (len_een_ij * sizeof(double)); + + if (context == QMCKL_NULL_CONTEXT) { + return QMCKL_INVALID_CONTEXT; + } + + if (walk_num <= 0) { + return QMCKL_INVALID_ARG_2; + } + + if (elec_num <= 0) { + return QMCKL_INVALID_ARG_3; + } + + if (cord_num <= 0) { + return QMCKL_INVALID_ARG_4; + } + + // Prepare table of exponentiated distances raised to appropriate power + // init + + for (int kk = 0; kk < walk_num*(cord_num+1)*elec_num*elec_num; ++kk) { + een_rescaled_e[kk]= 0.0; + } + + /* + for (int nw = 0; nw < walk_num; ++nw) { + for (int l = 0; l < (cord_num + 1); ++l) { + for (int i = 0; i < elec_num; ++i) { + for (int j = 0; j < elec_num; ++j) { + een_rescaled_e[j + i*elec_num + l*elec_num*elec_num + nw*(cord_num+1)*elec_num*elec_num]= 0.0; + } + } + } + } + */ + + for (int nw = 0; nw < walk_num; ++nw) { + + for (int kk = 0; kk < len_een_ij; ++kk) { + // this array initialized at 0 except een_rescaled_e_ij(:, 1) = 1.0d0 + // and the arrangement of indices is [cord_num+1, ne*(ne-1)/2] + een_rescaled_e_ij[kk]= ( kk < (elec_pairs) ? 1.0 : 0.0 ); + } + + k = 0; + for (int i = 0; i < elec_num; ++i) { + for (int j = 0; j < i; ++j) { + // een_rescaled_e_ij(k, 2) = dexp(-rescale_factor_kappa_ee * ee_distance(i, j, nw)); + een_rescaled_e_ij[k + elec_pairs] = exp(-rescale_factor_kappa_ee * \ + ee_distance[j + i*elec_num + nw*(elec_num*elec_num)]); + k = k + 1; + } + } + + + for (int l = 2; l < (cord_num+1); ++l) { + for (int k = 0; k < elec_pairs; ++k) { + // een_rescaled_e_ij(k, l + 1) = een_rescaled_e_ij(k, l + 1 - 1) * een_rescaled_e_ij(k, 2) + een_rescaled_e_ij[k+l*elec_pairs] = een_rescaled_e_ij[k + (l - 1)*elec_pairs] * \ + een_rescaled_e_ij[k + elec_pairs]; + } + } + + + // prepare the actual een table + for (int i = 0; i < elec_num; ++i){ + for (int j = 0; j < elec_num; ++j) { + een_rescaled_e[j + i*elec_num + 0 + nw*(cord_num+1)*elec_num*elec_num] = 1.0; + } + } + + // Up to here it should work. + for ( int l = 1; l < (cord_num+1); ++l) { + k = 0; + for (int i = 0; i < elec_num; ++i) { + for (int j = 0; j < i; ++j) { + x = een_rescaled_e_ij[k + l*elec_pairs]; + een_rescaled_e[j + i*elec_num + l*elec_num*elec_num + nw*elec_num*elec_num*(cord_num+1)] = x; + een_rescaled_e[i + j*elec_num + l*elec_num*elec_num + nw*elec_num*elec_num*(cord_num+1)] = x; + k = k + 1; + } + } + } + + for (int l = 0; l < (cord_num + 1); ++l) { + for (int j = 0; j < elec_num; ++j) { + een_rescaled_e[j + j*elec_num + l*elec_num*elec_num + nw*elec_num*elec_num*(cord_num+1)] = 0.0; + } + } + + } + + free(een_rescaled_e_ij); + + return QMCKL_SUCCESS; +} + #+end_src + +# #+CALL: generate_c_header(table=qmckl_factor_een_rescaled_e_args,rettyp=get_value("CRetType"),fname="qmckl_compute_een_rescaled_e_doc") + + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none + qmckl_exit_code qmckl_compute_een_rescaled_e ( + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t cord_num, + const double rescale_factor_kappa_ee, + const double* ee_distance, + double* const een_rescaled_e ); + #+end_src + + #+begin_src c :tangle (eval h_private_func) :comments org + qmckl_exit_code qmckl_compute_een_rescaled_e_doc ( + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t cord_num, + const double rescale_factor_kappa_ee, + const double* ee_distance, + double* const een_rescaled_e ); + #+end_src + + #+begin_src c :tangle (eval h_private_func) :comments org + qmckl_exit_code qmckl_compute_een_rescaled_e_hpc ( + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t cord_num, + const double rescale_factor_kappa_ee, + const double* ee_distance, + double* const een_rescaled_e ); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes + qmckl_exit_code qmckl_compute_een_rescaled_e ( + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t cord_num, + const double rescale_factor_kappa_ee, + const double* ee_distance, + double* const een_rescaled_e ) { + + #ifdef HAVE_HPC + return qmckl_compute_een_rescaled_e_hpc(context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, ee_distance, een_rescaled_e); + #else + return qmckl_compute_een_rescaled_e_doc(context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, ee_distance, een_rescaled_e); + #endif + } + #+end_src + + *** Test #+begin_src python :results output :exports none :noweb yes @@ -3376,7 +3776,6 @@ assert(fabs(een_rescaled_e[0][1][0][4]-0.01754273169464735) < 1.e-12); assert(fabs(een_rescaled_e[0][2][1][3]-0.02214680362033448) < 1.e-12); assert(fabs(een_rescaled_e[0][2][1][4]-0.0005700154999202759) < 1.e-12); assert(fabs(een_rescaled_e[0][2][1][5]-0.3424402276009091) < 1.e-12); - #+end_src ** Electron-electron rescaled distances for each order and derivatives @@ -3413,7 +3812,7 @@ qmckl_get_jastrow_een_rescaled_e_deriv_e(qmckl_context context, rc = qmckl_provide_een_rescaled_e_deriv_e(context); if (rc != QMCKL_SUCCESS) return rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int64_t sze = ctx->electron.num * 4 * ctx->electron.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1); @@ -3443,7 +3842,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_e_deriv_e(qmckl_context context) return QMCKL_NULL_CONTEXT; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); /* Check if ee distance is provided */ @@ -3470,16 +3869,15 @@ qmckl_exit_code qmckl_provide_een_rescaled_e_deriv_e(qmckl_context context) ctx->jastrow.een_rescaled_e_deriv_e = een_rescaled_e_deriv_e; } - qmckl_exit_code rc = - qmckl_compute_factor_een_rescaled_e_deriv_e(context, - ctx->electron.walk_num, - ctx->electron.num, - ctx->jastrow.cord_num, - ctx->electron.rescale_factor_kappa_ee, - ctx->electron.coord_new.data, - ctx->electron.ee_distance, - ctx->jastrow.een_rescaled_e, - ctx->jastrow.een_rescaled_e_deriv_e); + rc = qmckl_compute_factor_een_rescaled_e_deriv_e(context, + ctx->electron.walk_num, + ctx->electron.num, + ctx->jastrow.cord_num, + ctx->electron.rescale_factor_kappa_ee, + ctx->electron.coord_new.data, + ctx->electron.ee_distance, + ctx->jastrow.een_rescaled_e, + ctx->jastrow.een_rescaled_e_deriv_e); if (rc != QMCKL_SUCCESS) { return rc; } @@ -3493,7 +3891,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_e_deriv_e(qmckl_context context) *** Compute :PROPERTIES: - :Name: qmckl_compute_een_rescaled_e_deriv_e + :Name: qmckl_compute_factor_een_rescaled_e_deriv_e :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: @@ -3512,7 +3910,8 @@ qmckl_exit_code qmckl_provide_een_rescaled_e_deriv_e(qmckl_context context) | ~een_rescaled_e_deriv_e~ | ~double[walk_num][0:cord_num][elec_num][4][elec_num]~ | out | Electron-electron rescaled distances | #+begin_src f90 :comments org :tangle (eval f) :noweb yes -integer function qmckl_compute_factor_een_rescaled_e_deriv_e_f(context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, & +integer function qmckl_compute_factor_een_rescaled_e_deriv_e_f( & + context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, & coord_new, ee_distance, een_rescaled_e, een_rescaled_e_deriv_e) & result(info) use qmckl @@ -3599,21 +3998,20 @@ integer function qmckl_compute_factor_een_rescaled_e_deriv_e_f(context, walk_num end function qmckl_compute_factor_een_rescaled_e_deriv_e_f #+end_src - #+CALL: generate_c_header(table=qmckl_factor_een_rescaled_e_deriv_e_args,rettyp=get_value("CRetType"),fname=get_value("Name")) +# #+CALL: generate_c_header(table=qmckl_factor_een_rescaled_e_deriv_e_args,rettyp=get_value("CRetType"),fname=get_value("Name")) - #+RESULTS: - #+begin_src c :tangle (eval h_func) :comments org - qmckl_exit_code qmckl_compute_factor_een_rescaled_e_deriv_e ( - const qmckl_context context, - const int64_t walk_num, - const int64_t elec_num, - const int64_t cord_num, - const double rescale_factor_kappa_ee, - const double* coord_new, - const double* ee_distance, - const double* een_rescaled_e, - double* const een_rescaled_e_deriv_e ); - #+end_src + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none + qmckl_exit_code qmckl_compute_factor_een_rescaled_e_deriv_e ( + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t cord_num, + const double rescale_factor_kappa_ee, + const double* coord_new, + const double* ee_distance, + const double* een_rescaled_e, + double* const een_rescaled_e_deriv_e ); + #+end_src #+CALL: generate_c_interface(table=qmckl_factor_een_rescaled_e_deriv_e_args,rettyp=get_value("CRetType"),fname=get_value("Name")) @@ -3621,16 +4019,16 @@ end function qmckl_compute_factor_een_rescaled_e_deriv_e_f #+RESULTS: #+begin_src f90 :tangle (eval f) :comments org :exports none integer(c_int32_t) function qmckl_compute_factor_een_rescaled_e_deriv_e & - (context, & - walk_num, & - elec_num, & - cord_num, & - rescale_factor_kappa_ee, & - coord_new, & - ee_distance, & - een_rescaled_e, & - een_rescaled_e_deriv_e) & - bind(C) result(info) + (context, & + walk_num, & + elec_num, & + cord_num, & + rescale_factor_kappa_ee, & + coord_new, & + ee_distance, & + een_rescaled_e, & + een_rescaled_e_deriv_e) & + bind(C) result(info) use, intrinsic :: iso_c_binding implicit none @@ -3647,15 +4045,15 @@ end function qmckl_compute_factor_een_rescaled_e_deriv_e_f integer(c_int32_t), external :: qmckl_compute_factor_een_rescaled_e_deriv_e_f info = qmckl_compute_factor_een_rescaled_e_deriv_e_f & - (context, & - walk_num, & - elec_num, & - cord_num, & - rescale_factor_kappa_ee, & - coord_new, & - ee_distance, & - een_rescaled_e, & - een_rescaled_e_deriv_e) + (context, & + walk_num, & + elec_num, & + cord_num, & + rescale_factor_kappa_ee, & + coord_new, & + ee_distance, & + een_rescaled_e, & + een_rescaled_e_deriv_e) end function qmckl_compute_factor_een_rescaled_e_deriv_e #+end_src @@ -3792,7 +4190,7 @@ qmckl_get_jastrow_een_rescaled_n(qmckl_context context, rc = qmckl_provide_een_rescaled_n(context); if (rc != QMCKL_SUCCESS) return rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int64_t sze = ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1); @@ -3822,7 +4220,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_n(qmckl_context context) return QMCKL_NULL_CONTEXT; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); /* Check if ee distance is provided */ @@ -3849,15 +4247,14 @@ qmckl_exit_code qmckl_provide_een_rescaled_n(qmckl_context context) ctx->jastrow.een_rescaled_n = een_rescaled_n; } - qmckl_exit_code rc = - qmckl_compute_een_rescaled_n(context, - ctx->electron.walk_num, - ctx->electron.num, - ctx->nucleus.num, - ctx->jastrow.cord_num, - ctx->electron.rescale_factor_kappa_en, - ctx->electron.en_distance, - ctx->jastrow.een_rescaled_n); + rc = qmckl_compute_een_rescaled_n(context, + ctx->electron.walk_num, + ctx->electron.num, + ctx->nucleus.num, + ctx->jastrow.cord_num, + ctx->electron.rescale_factor_kappa_en, + ctx->electron.en_distance, + ctx->jastrow.een_rescaled_n); if (rc != QMCKL_SUCCESS) { return rc; } @@ -3889,7 +4286,8 @@ qmckl_exit_code qmckl_provide_een_rescaled_n(qmckl_context context) | ~een_rescaled_n~ | ~double[walk_num][0:cord_num][nucl_num][elec_num]~ | out | Electron-nucleus rescaled distances | #+begin_src f90 :comments org :tangle (eval f) :noweb yes -integer function qmckl_compute_een_rescaled_n_f(context, walk_num, elec_num, nucl_num, cord_num, rescale_factor_kappa_en, & +integer function qmckl_compute_een_rescaled_n_f( & + context, walk_num, elec_num, nucl_num, cord_num, rescale_factor_kappa_en, & en_distance, een_rescaled_n) & result(info) use qmckl @@ -3957,60 +4355,82 @@ integer function qmckl_compute_een_rescaled_n_f(context, walk_num, elec_num, nuc end function qmckl_compute_een_rescaled_n_f #+end_src - #+CALL: generate_c_header(table=qmckl_factor_een_rescaled_n_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + #+begin_src c :comments org :tangle (eval c) :noweb yes +qmckl_exit_code qmckl_compute_een_rescaled_n ( + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t nucl_num, + const int64_t cord_num, + const double rescale_factor_kappa_en, + const double* en_distance, + double* const een_rescaled_n ) { - #+RESULTS: - #+begin_src c :tangle (eval h_func) :comments org - qmckl_exit_code qmckl_compute_een_rescaled_n ( - const qmckl_context context, - const int64_t walk_num, - const int64_t elec_num, - const int64_t nucl_num, - const int64_t cord_num, - const double rescale_factor_kappa_en, - const double* en_distance, - double* const een_rescaled_n ); + + if (context == QMCKL_NULL_CONTEXT) { + return QMCKL_INVALID_CONTEXT; + } + + if (walk_num <= 0) { + return QMCKL_INVALID_ARG_2; + } + + if (elec_num <= 0) { + return QMCKL_INVALID_ARG_3; + } + + if (nucl_num <= 0) { + return QMCKL_INVALID_ARG_4; + } + + if (cord_num <= 0) { + return QMCKL_INVALID_ARG_5; + } + + // Prepare table of exponentiated distances raised to appropriate power + for (int i = 0; i < (walk_num*(cord_num+1)*nucl_num*elec_num); ++i) { + een_rescaled_n[i] = 17.0; + } + + for (int nw = 0; nw < walk_num; ++nw) { + for (int a = 0; a < nucl_num; ++a) { + for (int i = 0; i < elec_num; ++i) { + // prepare the actual een table + //een_rescaled_n(:, :, 0, nw) = 1.0d0 + een_rescaled_n[i + a * elec_num + 0 + nw * elec_num*nucl_num*(cord_num+1)] = 1.0; + //een_rescaled_n(i, a, 1, nw) = dexp(-rescale_factor_kappa_en * en_distance(i, a, nw)) + een_rescaled_n[i + a*elec_num + elec_num*nucl_num + nw*elec_num*nucl_num*(cord_num+1)] = exp(-rescale_factor_kappa_en * \ + en_distance[i + a*elec_num + nw*elec_num*nucl_num]); + } + } + + for (int l = 2; l < (cord_num+1); ++l){ + for (int a = 0; a < nucl_num; ++a) { + for (int i = 0; i < elec_num; ++i) { + een_rescaled_n[i + a*elec_num + l*elec_num*nucl_num + nw*elec_num*nucl_num*(cord_num+1)] = een_rescaled_n[i + a*elec_num + (l-1)*elec_num*nucl_num + nw*elec_num*nucl_num*(cord_num+1)] *\ + een_rescaled_n[i + a*elec_num + elec_num*nucl_num + nw*elec_num*nucl_num*(cord_num+1)]; + } + } + } + + } + + return QMCKL_SUCCESS; +} #+end_src - #+CALL: generate_c_interface(table=qmckl_factor_een_rescaled_n_args,rettyp=get_value("CRetType"),fname=get_value("Name")) +# #+CALL: generate_c_header(table=qmckl_factor_een_rescaled_n_args,rettyp=get_value("CRetType"),fname=get_value("Name")) - #+RESULTS: - #+begin_src f90 :tangle (eval f) :comments org :exports none - integer(c_int32_t) function qmckl_compute_een_rescaled_n & - (context, & - walk_num, & - elec_num, & - nucl_num, & - cord_num, & - rescale_factor_kappa_en, & - en_distance, & - een_rescaled_n) & - bind(C) result(info) - - use, intrinsic :: iso_c_binding - implicit none - - integer (c_int64_t) , intent(in) , value :: context - integer (c_int64_t) , intent(in) , value :: walk_num - integer (c_int64_t) , intent(in) , value :: elec_num - integer (c_int64_t) , intent(in) , value :: nucl_num - integer (c_int64_t) , intent(in) , value :: cord_num - real (c_double ) , intent(in) , value :: rescale_factor_kappa_en - real (c_double ) , intent(in) :: en_distance(nucl_num,elec_num,walk_num) - real (c_double ) , intent(out) :: een_rescaled_n(nucl_num,elec_num,0:cord_num,walk_num) - - integer(c_int32_t), external :: qmckl_compute_een_rescaled_n_f - info = qmckl_compute_een_rescaled_n_f & - (context, & - walk_num, & - elec_num, & - nucl_num, & - cord_num, & - rescale_factor_kappa_en, & - en_distance, & - een_rescaled_n) - - end function qmckl_compute_een_rescaled_n + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none + qmckl_exit_code qmckl_compute_een_rescaled_n ( + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t nucl_num, + const int64_t cord_num, + const double rescale_factor_kappa_en, + const double* en_distance, + double* const een_rescaled_n ); #+end_src *** Test @@ -4071,7 +4491,6 @@ assert(fabs(een_rescaled_n[0][1][0][4]-0.023391817607642338) < 1.e-12); assert(fabs(een_rescaled_n[0][2][1][3]-0.880957224822116) < 1.e-12); assert(fabs(een_rescaled_n[0][2][1][4]-0.027185942659395074) < 1.e-12); assert(fabs(een_rescaled_n[0][2][1][5]-0.01343938025140174) < 1.e-12); - #+end_src ** Electron-nucleus rescaled distances for each order and derivatives @@ -4104,7 +4523,7 @@ qmckl_get_jastrow_een_rescaled_n_deriv_e(qmckl_context context, rc = qmckl_provide_een_rescaled_n_deriv_e(context); if (rc != QMCKL_SUCCESS) return rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int64_t sze = ctx->electron.num * 4 * ctx->nucleus.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1); @@ -4134,7 +4553,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_n_deriv_e(qmckl_context context) return QMCKL_NULL_CONTEXT; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); /* Check if ee distance is provided */ @@ -4165,18 +4584,17 @@ qmckl_exit_code qmckl_provide_een_rescaled_n_deriv_e(qmckl_context context) ctx->jastrow.een_rescaled_n_deriv_e = een_rescaled_n_deriv_e; } - qmckl_exit_code rc = - qmckl_compute_factor_een_rescaled_n_deriv_e(context, - ctx->electron.walk_num, - ctx->electron.num, - ctx->nucleus.num, - ctx->jastrow.cord_num, - ctx->electron.rescale_factor_kappa_en, - ctx->electron.coord_new.data, - ctx->nucleus.coord.data, - ctx->electron.en_distance, - ctx->jastrow.een_rescaled_n, - ctx->jastrow.een_rescaled_n_deriv_e); + rc = qmckl_compute_factor_een_rescaled_n_deriv_e(context, + ctx->electron.walk_num, + ctx->electron.num, + ctx->nucleus.num, + ctx->jastrow.cord_num, + ctx->electron.rescale_factor_kappa_en, + ctx->electron.coord_new.data, + ctx->nucleus.coord.data, + ctx->electron.en_distance, + ctx->jastrow.een_rescaled_n, + ctx->jastrow.een_rescaled_n_deriv_e); if (rc != QMCKL_SUCCESS) { return rc; } @@ -4211,9 +4629,10 @@ qmckl_exit_code qmckl_provide_een_rescaled_n_deriv_e(qmckl_context context) | ~een_rescaled_n_deriv_e~ | ~double[walk_num][0:cord_num][nucl_num][4][elec_num]~ | out | Electron-nucleus rescaled distances | #+begin_src f90 :comments org :tangle (eval f) :noweb yes -integer function qmckl_compute_factor_een_rescaled_n_deriv_e_f(context, walk_num, elec_num, nucl_num, & - cord_num, rescale_factor_kappa_en, & - coord_new, coord, en_distance, een_rescaled_n, een_rescaled_n_deriv_e) & +integer function qmckl_compute_factor_een_rescaled_n_deriv_e_f( & + context, walk_num, elec_num, nucl_num, & + cord_num, rescale_factor_kappa_en, & + coord_new, coord, en_distance, een_rescaled_n, een_rescaled_n_deriv_e) & result(info) use qmckl implicit none @@ -4306,22 +4725,21 @@ integer function qmckl_compute_factor_een_rescaled_n_deriv_e_f(context, walk_num end function qmckl_compute_factor_een_rescaled_n_deriv_e_f #+end_src - #+CALL: generate_c_header(table=qmckl_compute_factor_een_rescaled_n_deriv_e_args,rettyp=get_value("CRetType"),fname=get_value("Name")) +# #+CALL: generate_c_header(table=qmckl_compute_factor_een_rescaled_n_deriv_e_args,rettyp=get_value("CRetType"),fname=get_value("Name")) - #+RESULTS: - #+begin_src c :tangle (eval h_func) :comments org + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_factor_een_rescaled_n_deriv_e ( - const qmckl_context context, - const int64_t walk_num, - const int64_t elec_num, - const int64_t nucl_num, - const int64_t cord_num, - const double rescale_factor_kappa_en, - const double* coord_new, - const double* coord, - const double* en_distance, - const double* een_rescaled_n, - double* const een_rescaled_n_deriv_e ); + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t nucl_num, + const int64_t cord_num, + const double rescale_factor_kappa_en, + const double* coord_new, + const double* coord, + const double* en_distance, + const double* een_rescaled_n, + double* const een_rescaled_n_deriv_e ); #+end_src #+CALL: generate_c_interface(table=qmckl_compute_factor_een_rescaled_n_deriv_e_args,rettyp=get_value("CRetType"),fname=get_value("Name")) @@ -4329,18 +4747,18 @@ end function qmckl_compute_factor_een_rescaled_n_deriv_e_f #+RESULTS: #+begin_src f90 :tangle (eval f) :comments org :exports none integer(c_int32_t) function qmckl_compute_factor_een_rescaled_n_deriv_e & - (context, & - walk_num, & - elec_num, & - nucl_num, & - cord_num, & - rescale_factor_kappa_en, & - coord_new, & - coord, & - en_distance, & - een_rescaled_n, & - een_rescaled_n_deriv_e) & - bind(C) result(info) + (context, & + walk_num, & + elec_num, & + nucl_num, & + cord_num, & + rescale_factor_kappa_en, & + coord_new, & + coord, & + en_distance, & + een_rescaled_n, & + een_rescaled_n_deriv_e) & + bind(C) result(info) use, intrinsic :: iso_c_binding implicit none @@ -4354,22 +4772,22 @@ end function qmckl_compute_factor_een_rescaled_n_deriv_e_f real (c_double ) , intent(in) :: coord_new(elec_num,3,walk_num) real (c_double ) , intent(in) :: coord(nucl_num,3) real (c_double ) , intent(in) :: en_distance(nucl_num,elec_num,walk_num) - real (c_double ) , intent(in) :: een_rescaled_n(0:cord_num,nucl_num,elec_num,walk_num) + real (c_double ) , intent(in) :: een_rescaled_n(elec_num,nucl_num,0:cord_num,walk_num) real (c_double ) , intent(out) :: een_rescaled_n_deriv_e(elec_num,4,nucl_num,0:cord_num,walk_num) integer(c_int32_t), external :: qmckl_compute_factor_een_rescaled_n_deriv_e_f info = qmckl_compute_factor_een_rescaled_n_deriv_e_f & - (context, & - walk_num, & - elec_num, & - nucl_num, & - cord_num, & - rescale_factor_kappa_en, & - coord_new, & - coord, & - en_distance, & - een_rescaled_n, & - een_rescaled_n_deriv_e) + (context, & + walk_num, & + elec_num, & + nucl_num, & + cord_num, & + rescale_factor_kappa_en, & + coord_new, & + coord, & + en_distance, & + een_rescaled_n, & + een_rescaled_n_deriv_e) end function qmckl_compute_factor_een_rescaled_n_deriv_e #+end_src @@ -4486,7 +4904,7 @@ qmckl_exit_code qmckl_get_jastrow_dim_cord_vect(qmckl_context context, int64_t* rc = qmckl_provide_dim_cord_vect(context); if (rc != QMCKL_SUCCESS) return rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); *dim_cord_vect = ctx->jastrow.dim_cord_vect; @@ -4508,7 +4926,7 @@ qmckl_exit_code qmckl_get_jastrow_cord_vect_full(qmckl_context context, double* rc = qmckl_provide_cord_vect_full(context); if (rc != QMCKL_SUCCESS) return rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); size_t sze = ctx->jastrow.dim_cord_vect * ctx->nucleus.num; @@ -4531,7 +4949,7 @@ qmckl_exit_code qmckl_get_jastrow_lkpm_combined_index(qmckl_context context, int rc = qmckl_provide_cord_vect_full(context); if (rc != QMCKL_SUCCESS) return rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); size_t sze = ctx->jastrow.dim_cord_vect * 4; @@ -4557,7 +4975,7 @@ qmckl_exit_code qmckl_get_jastrow_tmp_c(qmckl_context context, double* const tmp rc = qmckl_provide_tmp_c(context); if (rc != QMCKL_SUCCESS) return rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); size_t sze = (ctx->jastrow.cord_num) * (ctx->jastrow.cord_num + 1) @@ -4584,7 +5002,7 @@ qmckl_exit_code qmckl_get_jastrow_dtmp_c(qmckl_context context, double* const dt rc = qmckl_provide_dtmp_c(context); if (rc != QMCKL_SUCCESS) return rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); size_t sze = (ctx->jastrow.cord_num) * (ctx->jastrow.cord_num + 1) @@ -4613,7 +5031,7 @@ qmckl_exit_code qmckl_provide_dim_cord_vect(qmckl_context context) return QMCKL_NULL_CONTEXT; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); /* Compute if necessary */ @@ -4640,7 +5058,7 @@ qmckl_exit_code qmckl_provide_cord_vect_full(qmckl_context context) return QMCKL_NULL_CONTEXT; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); /* Check if dim_cord_vect is provided */ @@ -4666,14 +5084,13 @@ qmckl_exit_code qmckl_provide_cord_vect_full(qmckl_context context) ctx->jastrow.cord_vect_full = cord_vect_full; } - qmckl_exit_code rc = - qmckl_compute_cord_vect_full(context, - ctx->nucleus.num, - ctx->jastrow.dim_cord_vect, - ctx->jastrow.type_nucl_num, - ctx->jastrow.type_nucl_vector, - ctx->jastrow.cord_vector, - ctx->jastrow.cord_vect_full); + rc = qmckl_compute_cord_vect_full(context, + ctx->nucleus.num, + ctx->jastrow.dim_cord_vect, + ctx->jastrow.type_nucl_num, + ctx->jastrow.type_nucl_vector, + ctx->jastrow.cord_vector, + ctx->jastrow.cord_vect_full); if (rc != QMCKL_SUCCESS) { return rc; } @@ -4691,7 +5108,7 @@ qmckl_exit_code qmckl_provide_lkpm_combined_index(qmckl_context context) return QMCKL_NULL_CONTEXT; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); /* Check if dim_cord_vect is provided */ @@ -4717,11 +5134,10 @@ qmckl_exit_code qmckl_provide_lkpm_combined_index(qmckl_context context) ctx->jastrow.lkpm_combined_index = lkpm_combined_index; } - qmckl_exit_code rc = - qmckl_compute_lkpm_combined_index(context, - ctx->jastrow.cord_num, - ctx->jastrow.dim_cord_vect, - ctx->jastrow.lkpm_combined_index); + rc = qmckl_compute_lkpm_combined_index(context, + ctx->jastrow.cord_num, + ctx->jastrow.dim_cord_vect, + ctx->jastrow.lkpm_combined_index); if (rc != QMCKL_SUCCESS) { return rc; } @@ -4734,12 +5150,11 @@ qmckl_exit_code qmckl_provide_lkpm_combined_index(qmckl_context context) qmckl_exit_code qmckl_provide_tmp_c(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; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); /* Check if dim_cord_vect is provided */ @@ -4766,19 +5181,57 @@ qmckl_exit_code qmckl_provide_tmp_c(qmckl_context context) ctx->jastrow.tmp_c = tmp_c; } - qmckl_exit_code rc = - qmckl_compute_tmp_c(context, - ctx->jastrow.cord_num, - ctx->electron.num, - ctx->nucleus.num, - ctx->electron.walk_num, - ctx->jastrow.een_rescaled_e, - ctx->jastrow.een_rescaled_n, - ctx->jastrow.tmp_c); - if (rc != QMCKL_SUCCESS) { - return rc; + + /* Choose the correct compute function (depending on offload type) */ +#ifdef HAVE_HPC + const bool gpu_offload = ctx->jastrow.gpu_offload; +#else + const bool gpu_offload = false; +#endif + + if (gpu_offload) { +#ifdef HAVE_CUBLAS_OFFLOAD + rc = qmckl_compute_tmp_c_cublas_offload(context, + ctx->jastrow.cord_num, + ctx->electron.num, + ctx->nucleus.num, + ctx->electron.walk_num, + ctx->jastrow.een_rescaled_e, + ctx->jastrow.een_rescaled_n, + ctx->jastrow.tmp_c); +#elif HAVE_OPENACC_OFFLOAD + rc = qmckl_compute_tmp_c_acc_offload(context, + ctx->jastrow.cord_num, + ctx->electron.num, + ctx->nucleus.num, + ctx->electron.walk_num, + ctx->jastrow.een_rescaled_e, + ctx->jastrow.een_rescaled_n, + ctx->jastrow.tmp_c); +#elif HAVE_OPENMP_OFFLOAD + rc = qmckl_compute_tmp_c_omp_offload(context, + ctx->jastrow.cord_num, + ctx->electron.num, + ctx->nucleus.num, + ctx->electron.walk_num, + ctx->jastrow.een_rescaled_e, + ctx->jastrow.een_rescaled_n, + ctx->jastrow.tmp_c); +#else + rc = QMCKL_FAILURE; +#endif + } else { + rc = qmckl_compute_tmp_c(context, + ctx->jastrow.cord_num, + ctx->electron.num, + ctx->nucleus.num, + ctx->electron.walk_num, + ctx->jastrow.een_rescaled_e, + ctx->jastrow.een_rescaled_n, + ctx->jastrow.tmp_c); } + ctx->jastrow.tmp_c_date = ctx->date; } @@ -4787,12 +5240,11 @@ qmckl_exit_code qmckl_provide_tmp_c(qmckl_context context) qmckl_exit_code qmckl_provide_dtmp_c(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; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); /* Check if dim_cord_vect is provided */ @@ -4807,7 +5259,7 @@ qmckl_exit_code qmckl_provide_dtmp_c(qmckl_context context) qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = (ctx->jastrow.cord_num) * (ctx->jastrow.cord_num + 1) - * 4 * ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num * sizeof(double); + ,* 4 * ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num * sizeof(double); double* dtmp_c = (double*) qmckl_malloc(context, mem_info); if (dtmp_c == NULL) { @@ -4819,19 +5271,60 @@ qmckl_exit_code qmckl_provide_dtmp_c(qmckl_context context) ctx->jastrow.dtmp_c = dtmp_c; } - qmckl_exit_code rc = - qmckl_compute_dtmp_c(context, - ctx->jastrow.cord_num, - ctx->electron.num, - ctx->nucleus.num, - ctx->electron.walk_num, - ctx->jastrow.een_rescaled_e_deriv_e, - ctx->jastrow.een_rescaled_n, - ctx->jastrow.dtmp_c); + +#ifdef HAVE_HPC + const bool gpu_offload = ctx->jastrow.gpu_offload; +#else + const bool gpu_offload = false; +#endif + + if (gpu_offload) { +#ifdef HAVE_CUBLAS_OFFLOAD + rc = qmckl_compute_dtmp_c_cublas_offload(context, + ctx->jastrow.cord_num, + ctx->electron.num, + ctx->nucleus.num, + ctx->electron.walk_num, + ctx->jastrow.een_rescaled_e_deriv_e, + ctx->jastrow.een_rescaled_n, + ctx->jastrow.dtmp_c); +#elif HAVE_OPENACC_OFFLOAD + rc = qmckl_compute_dtmp_c_acc_offload(context, + ctx->jastrow.cord_num, + ctx->electron.num, + ctx->nucleus.num, + ctx->electron.walk_num, + ctx->jastrow.een_rescaled_e_deriv_e, + ctx->jastrow.een_rescaled_n, + ctx->jastrow.dtmp_c); +#elif HAVE_OPENMP_OFFLOAD + rc = qmckl_compute_dtmp_c_omp_offload(context, + ctx->jastrow.cord_num, + ctx->electron.num, + ctx->nucleus.num, + ctx->electron.walk_num, + ctx->jastrow.een_rescaled_e_deriv_e, + ctx->jastrow.een_rescaled_n, + ctx->jastrow.dtmp_c); +#else + rc = QMCKL_FAILURE; +#endif + } else { + rc = qmckl_compute_dtmp_c(context, + ctx->jastrow.cord_num, + ctx->electron.num, + ctx->nucleus.num, + ctx->electron.walk_num, + ctx->jastrow.een_rescaled_e_deriv_e, + ctx->jastrow.een_rescaled_n, + ctx->jastrow.dtmp_c); + } + if (rc != QMCKL_SUCCESS) { return rc; } + ctx->jastrow.dtmp_c_date = ctx->date; } @@ -4854,7 +5347,8 @@ qmckl_exit_code qmckl_provide_dtmp_c(qmckl_context context) | ~dim_cord_vect~ | ~int64_t~ | out | dimension of cord_vect_full table | #+begin_src f90 :comments org :tangle (eval f) :noweb yes -integer function qmckl_compute_dim_cord_vect_f(context, cord_num, dim_cord_vect) & +integer function qmckl_compute_dim_cord_vect_f( & + context, cord_num, dim_cord_vect) & result(info) use qmckl implicit none @@ -4895,39 +5389,53 @@ integer function qmckl_compute_dim_cord_vect_f(context, cord_num, dim_cord_vect) end function qmckl_compute_dim_cord_vect_f #+end_src - #+CALL: generate_c_header(table=qmckl_factor_dim_cord_vect_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + #+begin_src c :comments org :tangle (eval c) :noweb yes +qmckl_exit_code qmckl_compute_dim_cord_vect ( + const qmckl_context context, + const int64_t cord_num, + int64_t* const dim_cord_vect){ - #+RESULTS: - #+begin_src c :tangle (eval h_func) :comments org + int lmax; + + + if (context == QMCKL_NULL_CONTEXT) { + return QMCKL_INVALID_CONTEXT; + } + + if (cord_num <= 0) { + return QMCKL_INVALID_ARG_2; + } + + *dim_cord_vect = 0; + + for (int p=2; p <= cord_num; ++p){ + for (int k=p-1; k >= 0; --k) { + if (k != 0) { + lmax = p - k; + } else { + lmax = p - k - 2; + } + for (int l = lmax; l >= 0; --l) { + if ( ((p - k - l) & 1)==1) continue; + *dim_cord_vect=*dim_cord_vect+1; + } + } + } + + return QMCKL_SUCCESS; +} + #+end_src + +# #+CALL: generate_c_header(table=qmckl_factor_dim_cord_vect_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_dim_cord_vect ( - const qmckl_context context, - const int64_t cord_num, - int64_t* const dim_cord_vect ); + const qmckl_context context, + const int64_t cord_num, + int64_t* const dim_cord_vect ); #+end_src - #+CALL: generate_c_interface(table=qmckl_factor_dim_cord_vect_args,rettyp=get_value("CRetType"),fname=get_value("Name")) - - #+RESULTS: - #+begin_src f90 :tangle (eval f) :comments org :exports none - integer(c_int32_t) function qmckl_compute_dim_cord_vect & - (context, cord_num, dim_cord_vect) & - bind(C) result(info) - - use, intrinsic :: iso_c_binding - implicit none - - integer (c_int64_t) , intent(in) , value :: context - integer (c_int64_t) , intent(in) , value :: cord_num - integer (c_int64_t) , intent(out) :: dim_cord_vect - - integer(c_int32_t), external :: qmckl_compute_dim_cord_vect_f - info = qmckl_compute_dim_cord_vect_f & - (context, cord_num, dim_cord_vect) - - end function qmckl_compute_dim_cord_vect - #+end_src - *** Compute cord_vect_full :PROPERTIES: :Name: qmckl_compute_cord_vect_full @@ -4947,7 +5455,8 @@ end function qmckl_compute_dim_cord_vect_f | ~cord_vect_full~ | ~double[dim_cord_vect][nucl_num]~ | out | Full list of coefficients | #+begin_src f90 :comments org :tangle (eval f) :noweb yes -integer function qmckl_compute_cord_vect_full_f(context, nucl_num, dim_cord_vect, type_nucl_num, & +integer function qmckl_compute_cord_vect_full_doc_f( & + context, nucl_num, dim_cord_vect, type_nucl_num, & type_nucl_vector, cord_vector, cord_vect_full) & result(info) use qmckl @@ -4989,29 +5498,14 @@ integer function qmckl_compute_cord_vect_full_f(context, nucl_num, dim_cord_vect cord_vect_full(a,1:dim_cord_vect) = cord_vector(type_nucl_vector(a),1:dim_cord_vect) end do -end function qmckl_compute_cord_vect_full_f +end function qmckl_compute_cord_vect_full_doc_f #+end_src - #+CALL: generate_c_header(table=qmckl_factor_cord_vect_full_args,rettyp=get_value("CRetType"),fname=get_value("Name")) - - #+RESULTS: - #+begin_src c :tangle (eval h_func) :comments org - qmckl_exit_code qmckl_compute_cord_vect_full ( - const qmckl_context context, - const int64_t nucl_num, - const int64_t dim_cord_vect, - const int64_t type_nucl_num, - const int64_t* type_nucl_vector, - const double* cord_vector, - double* const cord_vect_full ); - #+end_src - - - #+CALL: generate_c_interface(table=qmckl_factor_cord_vect_full_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + #+CALL: generate_c_interface(table=qmckl_factor_cord_vect_full_args,rettyp=get_value("CRetType"),fname="qmckl_compute_cord_vect_full_doc") #+RESULTS: #+begin_src f90 :tangle (eval f) :comments org :exports none - integer(c_int32_t) function qmckl_compute_cord_vect_full & + integer(c_int32_t) function qmckl_compute_cord_vect_full_doc & (context, nucl_num, dim_cord_vect, type_nucl_num, type_nucl_vector, cord_vector, cord_vect_full) & bind(C) result(info) @@ -5026,13 +5520,105 @@ end function qmckl_compute_cord_vect_full_f real (c_double ) , intent(in) :: cord_vector(type_nucl_num,dim_cord_vect) real (c_double ) , intent(out) :: cord_vect_full(nucl_num,dim_cord_vect) - integer(c_int32_t), external :: qmckl_compute_cord_vect_full_f - info = qmckl_compute_cord_vect_full_f & + integer(c_int32_t), external :: qmckl_compute_cord_vect_full_doc_f + info = qmckl_compute_cord_vect_full_doc_f & (context, nucl_num, dim_cord_vect, type_nucl_num, type_nucl_vector, cord_vector, cord_vect_full) - end function qmckl_compute_cord_vect_full + end function qmckl_compute_cord_vect_full_doc #+end_src + #+begin_src c :comments org :tangle (eval c) :noweb yes +qmckl_exit_code qmckl_compute_cord_vect_full_hpc ( + const qmckl_context context, + const int64_t nucl_num, + const int64_t dim_cord_vect, + const int64_t type_nucl_num, + const int64_t* type_nucl_vector, + const double* cord_vector, + double* const cord_vect_full ) { + + if (context == QMCKL_NULL_CONTEXT) { + return QMCKL_INVALID_CONTEXT; + } + + if (nucl_num <= 0) { + return QMCKL_INVALID_ARG_2; + } + + if (type_nucl_num <= 0) { + return QMCKL_INVALID_ARG_4; + } + + if (dim_cord_vect <= 0) { + return QMCKL_INVALID_ARG_5; + } + + for (int i=0; i < dim_cord_vect; ++i) { + for (int a=0; a < nucl_num; ++a){ + cord_vect_full[a + i*nucl_num] = cord_vector[(type_nucl_vector[a]-1)+i*type_nucl_num]; + } + } + + return QMCKL_SUCCESS; + } + #+end_src + + +# #+CALL: generate_c_header(table=qmckl_factor_cord_vect_full_args,rettyp=get_value("CRetType"),fname="qmckl_compute_cord_vect_full_doc") + + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none + qmckl_exit_code qmckl_compute_cord_vect_full ( + const qmckl_context context, + const int64_t nucl_num, + const int64_t dim_cord_vect, + const int64_t type_nucl_num, + const int64_t* type_nucl_vector, + const double* cord_vector, + double* const cord_vect_full ); + #+end_src + + #+begin_src c :tangle (eval h_private_func) :comments org + qmckl_exit_code qmckl_compute_cord_vect_full_doc ( + const qmckl_context context, + const int64_t nucl_num, + const int64_t dim_cord_vect, + const int64_t type_nucl_num, + const int64_t* type_nucl_vector, + const double* cord_vector, + double* const cord_vect_full ); + #+end_src + + #+begin_src c :tangle (eval h_private_func) :comments org + qmckl_exit_code qmckl_compute_cord_vect_full_hpc ( + const qmckl_context context, + const int64_t nucl_num, + const int64_t dim_cord_vect, + const int64_t type_nucl_num, + const int64_t* type_nucl_vector, + const double* cord_vector, + double* const cord_vect_full ); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes +qmckl_exit_code qmckl_compute_cord_vect_full ( + const qmckl_context context, + const int64_t nucl_num, + const int64_t dim_cord_vect, + const int64_t type_nucl_num, + const int64_t* type_nucl_vector, + const double* cord_vector, + double* const cord_vect_full ) { + + #ifdef HAVE_HPC + return qmckl_compute_cord_vect_full_hpc(context, nucl_num, dim_cord_vect, type_nucl_num, type_nucl_vector, cord_vector, cord_vect_full); + #else + return qmckl_compute_cord_vect_full_doc(context, nucl_num, dim_cord_vect, type_nucl_num, type_nucl_vector, cord_vector, cord_vect_full); + #endif + } + #+end_src + + + *** Compute lkpm_combined_index :PROPERTIES: :Name: qmckl_compute_lkpm_combined_index @@ -5046,11 +5632,11 @@ end function qmckl_compute_cord_vect_full_f | ~context~ | ~qmckl_context~ | in | Global state | | ~cord_num~ | ~int64_t~ | in | Order of polynomials | | ~dim_cord_vect~ | ~int64_t~ | in | dimension of cord full table | - | ~lpkm_combined_index~ | ~int64_t[4][dim_cord_vect]~ | out | Full list of combined indices | + | ~lkpm_combined_index~ | ~int64_t[4][dim_cord_vect]~ | out | Full list of combined indices | #+begin_src f90 :comments org :tangle (eval f) :noweb yes -integer function qmckl_compute_lkpm_combined_index_f(context, cord_num, dim_cord_vect, & - lkpm_combined_index) & +integer function qmckl_compute_lkpm_combined_index_f( & + context, cord_num, dim_cord_vect, lkpm_combined_index) & result(info) use qmckl implicit none @@ -5102,40 +5688,64 @@ integer function qmckl_compute_lkpm_combined_index_f(context, cord_num, dim_cord end function qmckl_compute_lkpm_combined_index_f #+end_src - #+CALL: generate_c_header(table=qmckl_factor_lkpm_combined_index_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + #+begin_src c :comments org :tangle (eval c) :noweb yes +qmckl_exit_code qmckl_compute_lkpm_combined_index ( + const qmckl_context context, + const int64_t cord_num, + const int64_t dim_cord_vect, + int64_t* const lkpm_combined_index ) { - #+RESULTS: - #+begin_src c :tangle (eval h_func) :comments org + int kk, lmax, m; + + if (context == QMCKL_NULL_CONTEXT) { + return QMCKL_INVALID_CONTEXT; + } + + if (cord_num <= 0) { + return QMCKL_INVALID_ARG_2; + } + + if (dim_cord_vect <= 0) { + return QMCKL_INVALID_ARG_3; + } + +/* +*/ + kk = 0; + for (int p = 2; p <= cord_num; ++p) { + for (int k=(p-1); k >= 0; --k) { + if (k != 0) { + lmax = p - k; + } else { + lmax = p - k - 2; + } + for (int l=lmax; l >= 0; --l) { + if (((p - k - l) & 1) == 1) continue; + m = (p - k - l)/2; + lkpm_combined_index[kk ] = l; + lkpm_combined_index[kk + dim_cord_vect] = k; + lkpm_combined_index[kk + 2*dim_cord_vect] = p; + lkpm_combined_index[kk + 3*dim_cord_vect] = m; + kk = kk + 1; + } + } + } + + return QMCKL_SUCCESS; +} + #+end_src + +# #+CALL: generate_c_header(table=qmckl_factor_lkpm_combined_index_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_lkpm_combined_index ( - const qmckl_context context, - const int64_t cord_num, - const int64_t dim_cord_vect, - int64_t* const lpkm_combined_index ); + const qmckl_context context, + const int64_t cord_num, + const int64_t dim_cord_vect, + int64_t* const lkpm_combined_index ); #+end_src - #+CALL: generate_c_interface(table=qmckl_factor_lkpm_combined_index_args,rettyp=get_value("CRetType"),fname=get_value("Name")) - - #+RESULTS: - #+begin_src f90 :tangle (eval f) :comments org :exports none - integer(c_int32_t) function qmckl_compute_lkpm_combined_index & - (context, cord_num, dim_cord_vect, lpkm_combined_index) & - bind(C) result(info) - - use, intrinsic :: iso_c_binding - implicit none - - integer (c_int64_t) , intent(in) , value :: context - integer (c_int64_t) , intent(in) , value :: cord_num - integer (c_int64_t) , intent(in) , value :: dim_cord_vect - integer (c_int64_t) , intent(out) :: lpkm_combined_index(dim_cord_vect,4) - - integer(c_int32_t), external :: qmckl_compute_lkpm_combined_index_f - info = qmckl_compute_lkpm_combined_index_f & - (context, cord_num, dim_cord_vect, lpkm_combined_index) - - end function qmckl_compute_lkpm_combined_index - #+end_src *** Compute tmp_c :PROPERTIES: @@ -5156,8 +5766,41 @@ end function qmckl_compute_lkpm_combined_index_f | ~een_rescaled_n~ | ~double[walk_num][0:cord_num][nucl_num][elec_num]~ | in | Electron-nucleus rescaled factor | | ~tmp_c~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]~ | out | vector of non-zero coefficients | + #+begin_src c :comments org :tangle (eval c) :noweb yes +qmckl_exit_code qmckl_compute_tmp_c (const qmckl_context context, + const int64_t cord_num, + const int64_t elec_num, + const int64_t nucl_num, + const int64_t walk_num, + const double* een_rescaled_e, + const double* een_rescaled_n, + double* const tmp_c ) +{ +#ifdef HAVE_HPC + return qmckl_compute_tmp_c_hpc(context, cord_num, elec_num, nucl_num, walk_num, een_rescaled_e, een_rescaled_n, tmp_c); +#else + return qmckl_compute_tmp_c_doc(context, cord_num, elec_num, nucl_num, walk_num, een_rescaled_e, een_rescaled_n, tmp_c); +#endif +} + #+end_src + +# #+CALL: generate_c_header(table=qmckl_factor_tmp_c_args,rettyp=get_value("CRetType"),fname="qmckl_compute_tmp_c") + + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none + qmckl_exit_code qmckl_compute_tmp_c ( + const qmckl_context context, + const int64_t cord_num, + const int64_t elec_num, + const int64_t nucl_num, + const int64_t walk_num, + const double* een_rescaled_e, + const double* een_rescaled_n, + double* const tmp_c ); + #+end_src + #+begin_src f90 :comments org :tangle (eval f) :noweb yes -integer function qmckl_compute_tmp_c_f(context, cord_num, elec_num, nucl_num, & +integer function qmckl_compute_tmp_c_doc_f( & + context, cord_num, elec_num, nucl_num, & walk_num, een_rescaled_e, een_rescaled_n, tmp_c) & result(info) use qmckl @@ -5212,7 +5855,7 @@ integer function qmckl_compute_tmp_c_f(context, cord_num, elec_num, nucl_num, & do nw=1, walk_num do i=0, cord_num-1 - info = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & + info = qmckl_dgemm(context, TransA, TransB, M, N, K, alpha, & een_rescaled_e(1,1,i,nw),LDA*1_8, & een_rescaled_n(1,1,0,nw),LDB*1_8, & beta, & @@ -5220,52 +5863,450 @@ integer function qmckl_compute_tmp_c_f(context, cord_num, elec_num, nucl_num, & end do end do -end function qmckl_compute_tmp_c_f +end function qmckl_compute_tmp_c_doc_f #+end_src - #+CALL: generate_c_header(table=qmckl_factor_tmp_c_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + #+begin_src c :tangle (eval h_private_func) :comments org +qmckl_exit_code qmckl_compute_tmp_c_doc ( + const qmckl_context context, + const int64_t cord_num, + const int64_t elec_num, + const int64_t nucl_num, + const int64_t walk_num, + const double* een_rescaled_e, + const double* een_rescaled_n, + double* const tmp_c ); + #+end_src + + #+CALL: generate_c_interface(table=qmckl_factor_tmp_c_args,rettyp=get_value("FRetType"),fname="qmckl_compute_tmp_c_doc") + +#+RESULTS: +#+begin_src f90 :tangle (eval f) :comments org :exports none +integer(c_int32_t) function qmckl_compute_tmp_c_doc & + (context, cord_num, elec_num, nucl_num, walk_num, een_rescaled_e, een_rescaled_n, tmp_c) & + bind(C) result(info) + + use, intrinsic :: iso_c_binding + implicit none + + integer (c_int64_t) , intent(in) , value :: context + integer (c_int64_t) , intent(in) , value :: cord_num + integer (c_int64_t) , intent(in) , value :: elec_num + integer (c_int64_t) , intent(in) , value :: nucl_num + integer (c_int64_t) , intent(in) , value :: walk_num + real (c_double ) , intent(in) :: een_rescaled_e(elec_num,elec_num,0:cord_num,walk_num) + real (c_double ) , intent(in) :: een_rescaled_n(elec_num,nucl_num,0:cord_num,walk_num) + real (c_double ) , intent(out) :: tmp_c(elec_num,nucl_num,0:cord_num,0:cord_num-1,walk_num) + + integer(c_int32_t), external :: qmckl_compute_tmp_c_doc_f + info = qmckl_compute_tmp_c_doc_f & + (context, cord_num, elec_num, nucl_num, walk_num, een_rescaled_e, een_rescaled_n, tmp_c) + +end function qmckl_compute_tmp_c_doc +#+end_src + +**** CPU :noexport: + + #+begin_src c :comments org :tangle (eval c) :noweb yes +qmckl_exit_code qmckl_compute_tmp_c_hpc ( + const qmckl_context context, + const int64_t cord_num, + const int64_t elec_num, + const int64_t nucl_num, + const int64_t walk_num, + const double* een_rescaled_e, + const double* een_rescaled_n, + double* const tmp_c ) { + + if (context == QMCKL_NULL_CONTEXT) { + return QMCKL_INVALID_CONTEXT; + } + + if (cord_num <= 0) { + return QMCKL_INVALID_ARG_2; + } + + if (elec_num <= 0) { + return QMCKL_INVALID_ARG_3; + } + + if (nucl_num <= 0) { + return QMCKL_INVALID_ARG_4; + } + + if (walk_num <= 0) { + return QMCKL_INVALID_ARG_5; + } + + qmckl_exit_code info = QMCKL_SUCCESS; + + const char TransA = 'N'; + const char TransB = 'N'; + const double alpha = 1.0; + const double beta = 0.0; + + const int64_t M = elec_num; + const int64_t N = nucl_num*(cord_num + 1); + const int64_t K = elec_num; + + const int64_t LDA = elec_num; + const int64_t LDB = elec_num; + const int64_t LDC = elec_num; + + const int64_t af = elec_num*elec_num; + const int64_t bf = elec_num*nucl_num*(cord_num+1); + const int64_t cf = bf; + +#ifdef HAVE_OPENMP +#pragma omp parallel for collapse(2) +#endif + for (int64_t nw=0; nw < walk_num; ++nw) { + for (int64_t i=0; ielectron.walk_num; @@ -5540,7 +6994,7 @@ qmckl_exit_code qmckl_provide_factor_een(qmckl_context context) return QMCKL_NULL_CONTEXT; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); /* Check if en rescaled distance is provided */ @@ -5582,18 +7036,17 @@ qmckl_exit_code qmckl_provide_factor_een(qmckl_context context) ctx->jastrow.factor_een = factor_een; } - qmckl_exit_code rc = - qmckl_compute_factor_een(context, - ctx->electron.walk_num, - ctx->electron.num, - ctx->nucleus.num, - ctx->jastrow.cord_num, - ctx->jastrow.dim_cord_vect, - ctx->jastrow.cord_vect_full, - ctx->jastrow.lkpm_combined_index, - ctx->jastrow.tmp_c, - ctx->jastrow.een_rescaled_n, - ctx->jastrow.factor_een); + rc = qmckl_compute_factor_een(context, + ctx->electron.walk_num, + ctx->electron.num, + ctx->nucleus.num, + ctx->jastrow.cord_num, + ctx->jastrow.dim_cord_vect, + ctx->jastrow.cord_vect_full, + ctx->jastrow.lkpm_combined_index, + ctx->jastrow.tmp_c, + ctx->jastrow.een_rescaled_n, + ctx->jastrow.factor_een); if (rc != QMCKL_SUCCESS) { return rc; } @@ -5628,9 +7081,10 @@ qmckl_exit_code qmckl_provide_factor_een(qmckl_context context) | ~factor_een~ | ~double[walk_num]~ | out | Electron-nucleus jastrow | #+begin_src f90 :comments org :tangle (eval f) :noweb yes -integer function qmckl_compute_factor_een_naive_f(context, walk_num, elec_num, nucl_num, cord_num,& - dim_cord_vect, cord_vect_full, lkpm_combined_index, & - een_rescaled_e, een_rescaled_n, factor_een) & +integer function qmckl_compute_factor_een_naive_f( & + context, walk_num, elec_num, nucl_num, cord_num,& + dim_cord_vect, cord_vect_full, lkpm_combined_index, & + een_rescaled_e, een_rescaled_n, factor_een) & result(info) use qmckl implicit none @@ -5704,22 +7158,21 @@ integer function qmckl_compute_factor_een_naive_f(context, walk_num, elec_num, n end function qmckl_compute_factor_een_naive_f #+end_src - #+CALL: generate_c_header(table=qmckl_factor_een_naive_args,rettyp=get_value("CRetType"),fname=get_value("Name")) +# #+CALL: generate_c_header(table=qmckl_factor_een_naive_args,rettyp=get_value("CRetType"),fname=get_value("Name")) - #+RESULTS: - #+begin_src c :tangle (eval h_func) :comments org + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_factor_een_naive ( - const qmckl_context context, - const int64_t walk_num, - const int64_t elec_num, - const int64_t nucl_num, - const int64_t cord_num, - const int64_t dim_cord_vect, - const double* cord_vect_full, - const int64_t* lkpm_combined_index, - const double* een_rescaled_e, - const double* een_rescaled_n, - double* const factor_een ); + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t nucl_num, + const int64_t cord_num, + const int64_t dim_cord_vect, + const double* cord_vect_full, + const int64_t* lkpm_combined_index, + const double* een_rescaled_e, + const double* een_rescaled_n, + double* const factor_een ); #+end_src #+CALL: generate_c_interface(table=qmckl_factor_een_naive_args,rettyp=get_value("CRetType"),fname=get_value("Name")) @@ -5727,18 +7180,18 @@ end function qmckl_compute_factor_een_naive_f #+RESULTS: #+begin_src f90 :tangle (eval f) :comments org :exports none integer(c_int32_t) function qmckl_compute_factor_een_naive & - (context, & - walk_num, & - elec_num, & - nucl_num, & - cord_num, & - dim_cord_vect, & - cord_vect_full, & - lkpm_combined_index, & - een_rescaled_e, & - een_rescaled_n, & - factor_een) & - bind(C) result(info) + (context, & + walk_num, & + elec_num, & + nucl_num, & + cord_num, & + dim_cord_vect, & + cord_vect_full, & + lkpm_combined_index, & + een_rescaled_e, & + een_rescaled_n, & + factor_een) & + bind(C) result(info) use, intrinsic :: iso_c_binding implicit none @@ -5757,17 +7210,17 @@ end function qmckl_compute_factor_een_naive_f integer(c_int32_t), external :: qmckl_compute_factor_een_naive_f info = qmckl_compute_factor_een_naive_f & - (context, & - walk_num, & - elec_num, & - nucl_num, & - cord_num, & - dim_cord_vect, & - cord_vect_full, & - lkpm_combined_index, & - een_rescaled_e, & - een_rescaled_n, & - factor_een) + (context, & + walk_num, & + elec_num, & + nucl_num, & + cord_num, & + dim_cord_vect, & + cord_vect_full, & + lkpm_combined_index, & + een_rescaled_e, & + een_rescaled_n, & + factor_een) end function qmckl_compute_factor_een_naive #+end_src @@ -5795,9 +7248,10 @@ end function qmckl_compute_factor_een_naive_f | ~factor_een~ | ~double[walk_num]~ | out | Electron-nucleus jastrow | #+begin_src f90 :comments org :tangle (eval f) :noweb yes -integer function qmckl_compute_factor_een_f(context, walk_num, elec_num, nucl_num, cord_num, & - dim_cord_vect, cord_vect_full, lkpm_combined_index, & - tmp_c, een_rescaled_n, factor_een) & +integer function qmckl_compute_factor_een_f( & + context, walk_num, elec_num, nucl_num, cord_num, & + dim_cord_vect, cord_vect_full, lkpm_combined_index, & + tmp_c, een_rescaled_n, factor_een) & result(info) use qmckl implicit none @@ -5864,22 +7318,21 @@ integer function qmckl_compute_factor_een_f(context, walk_num, elec_num, nucl_nu end function qmckl_compute_factor_een_f #+end_src - #+CALL: generate_c_header(table=qmckl_factor_een_args,rettyp=get_value("CRetType"),fname=get_value("Name")) +# #+CALL: generate_c_header(table=qmckl_factor_een_args,rettyp=get_value("CRetType"),fname=get_value("Name")) - #+RESULTS: - #+begin_src c :tangle (eval h_func) :comments org + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_factor_een ( - const qmckl_context context, - const int64_t walk_num, - const int64_t elec_num, - const int64_t nucl_num, - const int64_t cord_num, - const int64_t dim_cord_vect, - const double* cord_vect_full, - const int64_t* lkpm_combined_index, - const double* een_rescaled_e, - const double* een_rescaled_n, - double* const factor_een ); + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t nucl_num, + const int64_t cord_num, + const int64_t dim_cord_vect, + const double* cord_vect_full, + const int64_t* lkpm_combined_index, + const double* een_rescaled_e, + const double* een_rescaled_n, + double* const factor_een ); #+end_src #+CALL: generate_c_interface(table=qmckl_factor_een_args,rettyp=get_value("CRetType"),fname=get_value("Name")) @@ -5887,18 +7340,18 @@ end function qmckl_compute_factor_een_f #+RESULTS: #+begin_src f90 :tangle (eval f) :comments org :exports none integer(c_int32_t) function qmckl_compute_factor_een & - (context, & - walk_num, & - elec_num, & - nucl_num, & - cord_num, & - dim_cord_vect, & - cord_vect_full, & - lkpm_combined_index, & - een_rescaled_e, & - een_rescaled_n, & - factor_een) & - bind(C) result(info) + (context, & + walk_num, & + elec_num, & + nucl_num, & + cord_num, & + dim_cord_vect, & + cord_vect_full, & + lkpm_combined_index, & + een_rescaled_e, & + een_rescaled_n, & + factor_een) & + bind(C) result(info) use, intrinsic :: iso_c_binding implicit none @@ -5917,17 +7370,17 @@ end function qmckl_compute_factor_een_f integer(c_int32_t), external :: qmckl_compute_factor_een_f info = qmckl_compute_factor_een_f & - (context, & - walk_num, & - elec_num, & - nucl_num, & - cord_num, & - dim_cord_vect, & - cord_vect_full, & - lkpm_combined_index, & - een_rescaled_e, & - een_rescaled_n, & - factor_een) + (context, & + walk_num, & + elec_num, & + nucl_num, & + cord_num, & + dim_cord_vect, & + cord_vect_full, & + lkpm_combined_index, & + een_rescaled_e, & + een_rescaled_n, & + factor_een) end function qmckl_compute_factor_een #+end_src @@ -5976,6 +7429,7 @@ double factor_een[walk_num]; rc = qmckl_get_jastrow_factor_een(context, &(factor_een[0]),walk_num); assert(fabs(factor_een[0] + 0.37407972141304213) < 1e-12); +return QMCKL_SUCCESS; #+end_src ** Electron-electron-nucleus Jastrow \(f_{een}\) derivative @@ -6008,7 +7462,7 @@ qmckl_get_jastrow_factor_een_deriv_e(qmckl_context context, rc = qmckl_provide_factor_een_deriv_e(context); if (rc != QMCKL_SUCCESS) return rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int64_t sze = ctx->electron.walk_num * 4 * ctx->electron.num; @@ -6039,7 +7493,7 @@ qmckl_exit_code qmckl_provide_factor_een_deriv_e(qmckl_context context) return QMCKL_NULL_CONTEXT; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); /* Check if en rescaled distance is provided */ @@ -6093,20 +7547,19 @@ qmckl_exit_code qmckl_provide_factor_een_deriv_e(qmckl_context context) ctx->jastrow.factor_een_deriv_e = factor_een_deriv_e; } - qmckl_exit_code rc = - qmckl_compute_factor_een_deriv_e(context, - ctx->electron.walk_num, - ctx->electron.num, - ctx->nucleus.num, - ctx->jastrow.cord_num, - ctx->jastrow.dim_cord_vect, - ctx->jastrow.cord_vect_full, - ctx->jastrow.lkpm_combined_index, - ctx->jastrow.tmp_c, - ctx->jastrow.dtmp_c, - ctx->jastrow.een_rescaled_n, - ctx->jastrow.een_rescaled_n_deriv_e, - ctx->jastrow.factor_een_deriv_e); + rc = qmckl_compute_factor_een_deriv_e(context, + ctx->electron.walk_num, + ctx->electron.num, + ctx->nucleus.num, + ctx->jastrow.cord_num, + ctx->jastrow.dim_cord_vect, + ctx->jastrow.cord_vect_full, + ctx->jastrow.lkpm_combined_index, + ctx->jastrow.tmp_c, + ctx->jastrow.dtmp_c, + ctx->jastrow.een_rescaled_n, + ctx->jastrow.een_rescaled_n_deriv_e, + ctx->jastrow.factor_een_deriv_e); if (rc != QMCKL_SUCCESS) { return rc; } @@ -6143,10 +7596,10 @@ qmckl_exit_code qmckl_provide_factor_een_deriv_e(qmckl_context context) | ~factor_een_deriv_e~ | ~double[walk_num][4][elec_num]~ | out | Electron-nucleus jastrow | #+begin_src f90 :comments org :tangle (eval f) :noweb yes -integer function qmckl_compute_factor_een_deriv_e_naive_f(context, walk_num, elec_num, nucl_num, cord_num, & - dim_cord_vect, cord_vect_full, lkpm_combined_index, & - een_rescaled_e, een_rescaled_n, & - een_rescaled_e_deriv_e, een_rescaled_n_deriv_e, factor_een_deriv_e)& +integer function qmckl_compute_factor_een_deriv_e_naive_f( & + context, walk_num, elec_num, nucl_num, cord_num, dim_cord_vect, & + cord_vect_full, lkpm_combined_index, een_rescaled_e, een_rescaled_n, & + een_rescaled_e_deriv_e, een_rescaled_n_deriv_e, factor_een_deriv_e)& result(info) use qmckl implicit none @@ -6239,24 +7692,23 @@ integer function qmckl_compute_factor_een_deriv_e_naive_f(context, walk_num, ele end function qmckl_compute_factor_een_deriv_e_naive_f #+end_src - #+CALL: generate_c_header(table=qmckl_factor_een_deriv_e_naive_args,rettyp=get_value("CRetType"),fname=get_value("Name")) +# #+CALL: generate_c_header(table=qmckl_factor_een_deriv_e_naive_args,rettyp=get_value("CRetType"),fname=get_value("Name")) - #+RESULTS: - #+begin_src c :tangle (eval h_func) :comments org + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_factor_een_deriv_e_naive ( - const qmckl_context context, - const int64_t walk_num, - const int64_t elec_num, - const int64_t nucl_num, - const int64_t cord_num, - const int64_t dim_cord_vect, - const double* cord_vect_full, - const int64_t* lkpm_combined_index, - const double* een_rescaled_e, - const double* een_rescaled_n, - const double* een_rescaled_e_deriv_e, - const double* een_rescaled_n_deriv_e, - double* const factor_een_deriv_e ); + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t nucl_num, + const int64_t cord_num, + const int64_t dim_cord_vect, + const double* cord_vect_full, + const int64_t* lkpm_combined_index, + const double* een_rescaled_e, + const double* een_rescaled_n, + const double* een_rescaled_e_deriv_e, + const double* een_rescaled_n_deriv_e, + double* const factor_een_deriv_e ); #+end_src @@ -6265,20 +7717,20 @@ end function qmckl_compute_factor_een_deriv_e_naive_f #+RESULTS: #+begin_src f90 :tangle (eval f) :comments org :exports none integer(c_int32_t) function qmckl_compute_factor_een_deriv_e_naive & - (context, & - walk_num, & - elec_num, & - nucl_num, & - cord_num, & - dim_cord_vect, & - cord_vect_full, & - lkpm_combined_index, & - een_rescaled_e, & - een_rescaled_n, & - een_rescaled_e_deriv_e, & - een_rescaled_n_deriv_e, & - factor_een_deriv_e) & - bind(C) result(info) + (context, & + walk_num, & + elec_num, & + nucl_num, & + cord_num, & + dim_cord_vect, & + cord_vect_full, & + lkpm_combined_index, & + een_rescaled_e, & + een_rescaled_n, & + een_rescaled_e_deriv_e, & + een_rescaled_n_deriv_e, & + factor_een_deriv_e) & + bind(C) result(info) use, intrinsic :: iso_c_binding implicit none @@ -6299,19 +7751,19 @@ end function qmckl_compute_factor_een_deriv_e_naive_f integer(c_int32_t), external :: qmckl_compute_factor_een_deriv_e_naive_f info = qmckl_compute_factor_een_deriv_e_naive_f & - (context, & - walk_num, & - elec_num, & - nucl_num, & - cord_num, & - dim_cord_vect, & - cord_vect_full, & - lkpm_combined_index, & - een_rescaled_e, & - een_rescaled_n, & - een_rescaled_e_deriv_e, & - een_rescaled_n_deriv_e, & - factor_een_deriv_e) + (context, & + walk_num, & + elec_num, & + nucl_num, & + cord_num, & + dim_cord_vect, & + cord_vect_full, & + lkpm_combined_index, & + een_rescaled_e, & + een_rescaled_n, & + een_rescaled_e_deriv_e, & + een_rescaled_n_deriv_e, & + factor_een_deriv_e) end function qmckl_compute_factor_een_deriv_e_naive #+end_src @@ -6342,9 +7794,10 @@ end function qmckl_compute_factor_een_deriv_e_naive_f #+begin_src f90 :comments org :tangle (eval f) :noweb yes -integer function qmckl_compute_factor_een_deriv_e_f(context, walk_num, elec_num, nucl_num, cord_num, dim_cord_vect, & - cord_vect_full, lkpm_combined_index, & - tmp_c, dtmp_c, een_rescaled_n, een_rescaled_n_deriv_e, factor_een_deriv_e) & +integer function qmckl_compute_factor_een_deriv_e_f( & + context, walk_num, elec_num, nucl_num, & + cord_num, dim_cord_vect, cord_vect_full, lkpm_combined_index, & + tmp_c, dtmp_c, een_rescaled_n, een_rescaled_n_deriv_e, factor_een_deriv_e)& result(info) use qmckl implicit none @@ -6430,24 +7883,23 @@ integer function qmckl_compute_factor_een_deriv_e_f(context, walk_num, elec_num, end function qmckl_compute_factor_een_deriv_e_f #+end_src - #+CALL: generate_c_header(table=qmckl_factor_een_deriv_e_args,rettyp=get_value("CRetType"),fname=get_value("Name")) +# #+CALL: generate_c_header(table=qmckl_factor_een_deriv_e_args,rettyp=get_value("CRetType"),fname=get_value("Name")) - #+RESULTS: - #+begin_src c :tangle (eval h_func) :comments org + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_factor_een_deriv_e ( - const qmckl_context context, - const int64_t walk_num, - const int64_t elec_num, - const int64_t nucl_num, - const int64_t cord_num, - const int64_t dim_cord_vect, - const double* cord_vect_full, - const int64_t* lkpm_combined_index, - const double* tmp_c, - const double* dtmp_c, - const double* een_rescaled_n, - const double* een_rescaled_n_deriv_e, - double* const factor_een_deriv_e ); + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t nucl_num, + const int64_t cord_num, + const int64_t dim_cord_vect, + const double* cord_vect_full, + const int64_t* lkpm_combined_index, + const double* tmp_c, + const double* dtmp_c, + const double* een_rescaled_n, + const double* een_rescaled_n_deriv_e, + double* const factor_een_deriv_e ); #+end_src @@ -6456,20 +7908,20 @@ end function qmckl_compute_factor_een_deriv_e_f #+RESULTS: #+begin_src f90 :tangle (eval f) :comments org :exports none integer(c_int32_t) function qmckl_compute_factor_een_deriv_e & - (context, & - walk_num, & - elec_num, & - nucl_num, & - cord_num, & - dim_cord_vect, & - cord_vect_full, & - lkpm_combined_index, & - tmp_c, & - dtmp_c, & - een_rescaled_n, & - een_rescaled_n_deriv_e, & - factor_een_deriv_e) & - bind(C) result(info) + (context, & + walk_num, & + elec_num, & + nucl_num, & + cord_num, & + dim_cord_vect, & + cord_vect_full, & + lkpm_combined_index, & + tmp_c, & + dtmp_c, & + een_rescaled_n, & + een_rescaled_n_deriv_e, & + factor_een_deriv_e) & + bind(C) result(info) use, intrinsic :: iso_c_binding implicit none @@ -6490,19 +7942,19 @@ end function qmckl_compute_factor_een_deriv_e_f integer(c_int32_t), external :: qmckl_compute_factor_een_deriv_e_f info = qmckl_compute_factor_een_deriv_e_f & - (context, & - walk_num, & - elec_num, & - nucl_num, & - cord_num, & - dim_cord_vect, & - cord_vect_full, & - lkpm_combined_index, & - tmp_c, & - dtmp_c, & - een_rescaled_n, & - een_rescaled_n_deriv_e, & - factor_een_deriv_e) + (context, & + walk_num, & + elec_num, & + nucl_num, & + cord_num, & + dim_cord_vect, & + cord_vect_full, & + lkpm_combined_index, & + tmp_c, & + dtmp_c, & + een_rescaled_n, & + een_rescaled_n_deriv_e, & + factor_een_deriv_e) end function qmckl_compute_factor_een_deriv_e #+end_src diff --git a/org/qmckl_local_energy.org b/org/qmckl_local_energy.org index aaad9f2..1ee61e0 100644 --- a/org/qmckl_local_energy.org +++ b/org/qmckl_local_energy.org @@ -226,7 +226,7 @@ qmckl_exit_code qmckl_get_kinetic_energy(qmckl_context context, double * const k rc = qmckl_provide_kinetic_energy(context); if (rc != QMCKL_SUCCESS) return rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); size_t sze = ctx->electron.walk_num * sizeof(double); @@ -250,7 +250,7 @@ qmckl_exit_code qmckl_provide_kinetic_energy(qmckl_context context) { return QMCKL_NULL_CONTEXT; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); if(!(ctx->nucleus.provided)) { @@ -549,37 +549,28 @@ end function qmckl_compute_kinetic_energy_f *** Test #+begin_src c :tangle (eval c_test) :exports none - -#define walk_num chbrclf_walk_num -#define elec_num chbrclf_elec_num -#define shell_num chbrclf_shell_num -#define ao_num chbrclf_ao_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]); -const int64_t nucl_num = chbrclf_nucl_num; const double* nucl_charge = chbrclf_charge; const double* nucl_coord = &(chbrclf_nucl_coord[0][0]); -rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num); +rc = qmckl_set_electron_num (context, chbrclf_elec_up_num, chbrclf_elec_dn_num); assert (rc == QMCKL_SUCCESS); -rc = qmckl_set_electron_walk_num (context, walk_num); +rc = qmckl_set_electron_walk_num (context, chbrclf_walk_num); assert (rc == QMCKL_SUCCESS); assert(qmckl_electron_provided(context)); -rc = qmckl_set_electron_coord (context, 'N', elec_coord, walk_num*elec_num*3); +rc = qmckl_set_electron_coord (context, 'N', elec_coord, chbrclf_walk_num*chbrclf_elec_num*3); assert(rc == QMCKL_SUCCESS); -rc = qmckl_set_nucleus_num (context, nucl_num); +rc = qmckl_set_nucleus_num (context, chbrclf_nucl_num); assert(rc == QMCKL_SUCCESS); -rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), nucl_num*3); +rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), chbrclf_nucl_num*3); assert(rc == QMCKL_SUCCESS); -rc = qmckl_set_nucleus_charge(context, nucl_charge, nucl_num); +rc = qmckl_set_nucleus_charge(context, nucl_charge, chbrclf_nucl_num); assert(rc == QMCKL_SUCCESS); assert(qmckl_nucleus_provided(context)); @@ -611,11 +602,11 @@ rc = qmckl_set_ao_basis_prim_num (context, chbrclf_prim_num); assert(rc == QMCKL_SUCCESS); assert(!qmckl_ao_basis_provided(context)); -rc = qmckl_set_ao_basis_nucleus_index (context, nucleus_index, nucl_num); +rc = qmckl_set_ao_basis_nucleus_index (context, nucleus_index, chbrclf_nucl_num); assert(rc == QMCKL_SUCCESS); assert(!qmckl_ao_basis_provided(context)); -rc = qmckl_set_ao_basis_nucleus_shell_num (context, nucleus_shell_num, nucl_num); +rc = qmckl_set_ao_basis_nucleus_shell_num (context, nucleus_shell_num, chbrclf_nucl_num); assert(rc == QMCKL_SUCCESS); assert(!qmckl_ao_basis_provided(context)); @@ -655,10 +646,10 @@ assert(rc == QMCKL_SUCCESS); assert(qmckl_ao_basis_provided(context)); -double ao_vgl[walk_num*elec_num][5][chbrclf_ao_num]; +double ao_vgl[chbrclf_walk_num*chbrclf_elec_num][5][chbrclf_ao_num]; rc = qmckl_get_ao_basis_ao_vgl(context, &(ao_vgl[0][0][0]), - (int64_t) 5*walk_num*elec_num*chbrclf_ao_num); + (int64_t) 5*chbrclf_walk_num*chbrclf_elec_num*chbrclf_ao_num); assert (rc == QMCKL_SUCCESS); /* Set up MO data */ @@ -673,31 +664,31 @@ assert (rc == QMCKL_SUCCESS); assert(qmckl_mo_basis_provided(context)); -double mo_vgl[walk_num*elec_num][5][chbrclf_mo_num]; -rc = qmckl_get_mo_basis_mo_vgl(context, &(mo_vgl[0][0][0]), 5*walk_num*elec_num*chbrclf_mo_num); +double mo_vgl[chbrclf_walk_num*chbrclf_elec_num][5][chbrclf_mo_num]; +rc = qmckl_get_mo_basis_mo_vgl(context, &(mo_vgl[0][0][0]), 5*chbrclf_walk_num*chbrclf_elec_num*chbrclf_mo_num); assert (rc == QMCKL_SUCCESS); /* Set up determinant data */ -const int64_t det_num_alpha = 1; -const int64_t det_num_beta = 1; -int64_t mo_index_alpha[det_num_alpha][walk_num][elec_up_num]; -int64_t mo_index_beta[det_num_alpha][walk_num][elec_dn_num]; +#define det_num_alpha 1 +#define det_num_beta 1 +int64_t mo_index_alpha[det_num_alpha][chbrclf_walk_num][chbrclf_elec_up_num]; +int64_t mo_index_beta[det_num_alpha][chbrclf_walk_num][chbrclf_elec_dn_num]; int i, j, k; for(k = 0; k < det_num_alpha; ++k) - for(i = 0; i < walk_num; ++i) - for(j = 0; j < elec_up_num; ++j) + for(i = 0; i < chbrclf_walk_num; ++i) + for(j = 0; j < chbrclf_elec_up_num; ++j) mo_index_alpha[k][i][j] = j + 1; for(k = 0; k < det_num_beta; ++k) - for(i = 0; i < walk_num; ++i) - for(j = 0; j < elec_up_num; ++j) + for(i = 0; i < chbrclf_walk_num; ++i) + for(j = 0; j < chbrclf_elec_up_num; ++j) mo_index_beta[k][i][j] = j + 1; rc = qmckl_set_determinant_type (context, typ); assert(rc == QMCKL_SUCCESS); -rc = qmckl_set_determinant_walk_num (context, walk_num); +rc = qmckl_set_determinant_walk_num (context, chbrclf_walk_num); assert (rc == QMCKL_SUCCESS); rc = qmckl_set_determinant_det_num_alpha (context, det_num_alpha); @@ -714,8 +705,8 @@ assert (rc == QMCKL_SUCCESS); // Get alpha determinant -double det_vgl_alpha[det_num_alpha][walk_num][5][elec_up_num][elec_up_num]; -double det_vgl_beta[det_num_beta][walk_num][5][elec_dn_num][elec_dn_num]; +double det_vgl_alpha[det_num_alpha][chbrclf_walk_num][5][chbrclf_elec_up_num][chbrclf_elec_up_num]; +double det_vgl_beta[det_num_beta][chbrclf_walk_num][5][chbrclf_elec_dn_num][chbrclf_elec_dn_num]; rc = qmckl_get_det_vgl_alpha(context, &(det_vgl_alpha[0][0][0][0][0])); assert (rc == QMCKL_SUCCESS); @@ -725,8 +716,8 @@ assert (rc == QMCKL_SUCCESS); // Get adjoint of the slater-determinant -double det_inv_matrix_alpha[det_num_alpha][walk_num][elec_up_num][elec_up_num]; -double det_inv_matrix_beta[det_num_beta][walk_num][elec_dn_num][elec_dn_num]; +double det_inv_matrix_alpha[det_num_alpha][chbrclf_walk_num][chbrclf_elec_up_num][chbrclf_elec_up_num]; +double det_inv_matrix_beta[det_num_beta][chbrclf_walk_num][chbrclf_elec_dn_num][chbrclf_elec_dn_num]; rc = qmckl_get_det_inv_matrix_alpha(context, &(det_inv_matrix_alpha[0][0][0][0])); assert (rc == QMCKL_SUCCESS); @@ -736,7 +727,7 @@ assert (rc == QMCKL_SUCCESS); // Calculate the Kinetic energy -double kinetic_energy[walk_num]; +double kinetic_energy[chbrclf_walk_num]; rc = qmckl_get_kinetic_energy(context, &(kinetic_energy[0])); assert (rc == QMCKL_SUCCESS); @@ -799,7 +790,7 @@ qmckl_exit_code qmckl_get_potential_energy(qmckl_context context, double * const rc = qmckl_provide_potential_energy(context); if (rc != QMCKL_SUCCESS) return rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); size_t sze = ctx->electron.walk_num * sizeof(double); @@ -822,7 +813,7 @@ qmckl_exit_code qmckl_provide_potential_energy(qmckl_context context) { return QMCKL_NULL_CONTEXT; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); qmckl_exit_code rc; @@ -1034,7 +1025,7 @@ end function qmckl_compute_potential_energy_f #+begin_src c :tangle (eval c_test) :exports none // Calculate the Potential energy -double potential_energy[walk_num]; +double potential_energy[chbrclf_walk_num]; rc = qmckl_get_potential_energy(context, &(potential_energy[0])); assert (rc == QMCKL_SUCCESS); @@ -1083,7 +1074,7 @@ qmckl_exit_code qmckl_get_local_energy(qmckl_context context, double * const loc rc = qmckl_provide_local_energy(context); if (rc != QMCKL_SUCCESS) return rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); const int64_t sze = ctx->electron.walk_num; @@ -1112,7 +1103,7 @@ qmckl_exit_code qmckl_provide_local_energy(qmckl_context context) { return QMCKL_NULL_CONTEXT; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); if(!(ctx->nucleus.provided)) { @@ -1296,7 +1287,7 @@ end function qmckl_compute_local_energy_f #+begin_src c :tangle (eval c_test) :exports none // Calculate the Local energy -double local_energy[walk_num]; +double local_energy[chbrclf_walk_num]; rc = qmckl_get_local_energy(context, &(local_energy[0]), walk_num); assert (rc == QMCKL_SUCCESS); @@ -1345,7 +1336,7 @@ qmckl_exit_code qmckl_get_drift_vector(qmckl_context context, double * const dri rc = qmckl_provide_drift_vector(context); if (rc != QMCKL_SUCCESS) return rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); size_t sze = ctx->electron.walk_num * ctx->electron.num * 3 * sizeof(double); @@ -1368,7 +1359,7 @@ qmckl_exit_code qmckl_provide_drift_vector(qmckl_context context) { return QMCKL_NULL_CONTEXT; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); if(!(ctx->nucleus.provided)) { @@ -1651,7 +1642,7 @@ end function qmckl_compute_drift_vector_f #+begin_src c :tangle (eval c_test) :exports none // Calculate the Drift vector -double drift_vector[walk_num][elec_num][3]; +double drift_vector[chbrclf_walk_num][chbrclf_elec_num][3]; rc = qmckl_get_drift_vector(context, &(drift_vector[0][0][0])); assert (rc == QMCKL_SUCCESS); diff --git a/org/qmckl_memory.org b/org/qmckl_memory.org index b24de98..1bc7006 100644 --- a/org/qmckl_memory.org +++ b/org/qmckl_memory.org @@ -116,7 +116,7 @@ void* qmckl_malloc(qmckl_context context, const qmckl_memory_info_struct info) { assert (qmckl_context_check(context) != QMCKL_NULL_CONTEXT); - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; /* Allocate memory and zero it */ void * pointer = malloc(info.size); @@ -217,7 +217,7 @@ qmckl_exit_code qmckl_free(qmckl_context context, void * const ptr) { "NULL pointer"); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; qmckl_lock(context); { diff --git a/org/qmckl_mo.org b/org/qmckl_mo.org index 88fe69c..4a148c6 100644 --- a/org/qmckl_mo.org +++ b/org/qmckl_mo.org @@ -92,10 +92,12 @@ int main() { Computed data: - |---------------+--------------------------+-------------------------------------------------------------------------------------| - | ~mo_vgl~ | ~[point_num][5][mo_num]~ | Value, gradients, Laplacian of the MOs at point positions | - | ~mo_vgl_date~ | ~uint64_t~ | Late modification date of Value, gradients, Laplacian of the MOs at point positions | - |---------------+--------------------------+-------------------------------------------------------------------------------------| + |-----------------+--------------------------+-------------------------------------------------------------------------------------| + | ~mo_value~ | ~[point_num][mo_num]~ | Value of the MOs at point positions | + | ~mo_value_date~ | ~uint64_t~ | Late modification date of the value of the MOs at point positions | + | ~mo_vgl~ | ~[point_num][5][mo_num]~ | Value, gradients, Laplacian of the MOs at point positions | + | ~mo_vgl_date~ | ~uint64_t~ | Late modification date of Value, gradients, Laplacian of the MOs at point positions | + |-----------------+--------------------------+-------------------------------------------------------------------------------------| ** Data structure @@ -106,7 +108,9 @@ typedef struct qmckl_mo_basis_struct { double * restrict coefficient_t; double * restrict mo_vgl; + double * restrict mo_value; uint64_t mo_vgl_date; + uint64_t mo_value_date; int32_t uninitialized; bool provided; @@ -131,7 +135,7 @@ qmckl_exit_code qmckl_init_mo_basis(qmckl_context context) { return false; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); ctx->mo_basis.uninitialized = (1 << 2) - 1; @@ -158,10 +162,9 @@ qmckl_get_mo_basis_mo_num (const qmckl_context context, QMCKL_INVALID_CONTEXT, "qmckl_get_mo_basis_mo_num", NULL); - return (int64_t) 0; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int32_t mask = 1; @@ -200,7 +203,7 @@ qmckl_get_mo_basis_coefficient (const qmckl_context context, NULL); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int32_t mask = 1 << 1; @@ -248,7 +251,7 @@ bool qmckl_mo_basis_provided(const qmckl_context context) { return false; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); return ctx->mo_basis.provided; @@ -257,10 +260,8 @@ bool qmckl_mo_basis_provided(const qmckl_context context) { #+end_src - *** Fortran interfaces - #+begin_src f90 :tangle (eval fh_func) :comments org interface integer(c_int32_t) function qmckl_get_mo_basis_mo_num (context, & mo_num) bind(C) @@ -280,7 +281,7 @@ interface implicit none integer (c_int64_t) , intent(in) , value :: context double precision, intent(out) :: coefficient(*) - integer (c_int64_t) , intent(in) , value :: size_max + integer (c_int64_t) , intent(int), value :: size_max end function qmckl_get_mo_basis_coefficient end interface @@ -302,7 +303,7 @@ if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } -qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; +qmckl_context_struct* const ctx = (qmckl_context_struct*) context; #+end_src #+NAME:post @@ -383,7 +384,7 @@ qmckl_exit_code qmckl_finalize_mo_basis(qmckl_context context) { NULL); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; @@ -421,7 +422,464 @@ qmckl_exit_code qmckl_finalize_mo_basis(qmckl_context context) { * Computation -** Computation of MOs +** Computation of MOs: values only + +*** Get + + #+begin_src c :comments org :tangle (eval h_func) :noweb yes +qmckl_exit_code +qmckl_get_mo_basis_mo_value(qmckl_context context, + double* const mo_value, + const int64_t size_max); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code +qmckl_get_mo_basis_mo_value(qmckl_context context, + double* const mo_value, + const int64_t size_max) +{ + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + + qmckl_exit_code rc; + + rc = qmckl_provide_ao_value(context); + if (rc != QMCKL_SUCCESS) return rc; + + rc = qmckl_provide_mo_value(context); + if (rc != QMCKL_SUCCESS) return rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; + assert (ctx != NULL); + + const int64_t sze = ctx->point.num * ctx->mo_basis.mo_num; + if (size_max < sze) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_get_mo_basis_mo_value", + "input array too small"); + } + memcpy(mo_value, ctx->mo_basis.mo_value, sze * sizeof(double)); + + return QMCKL_SUCCESS; +} + #+end_src + + #+begin_src f90 :tangle (eval fh_func) :comments org :exports none + interface + integer(c_int32_t) function qmckl_get_mo_basis_mo_value (context, & + mo_value, size_max) bind(C) + use, intrinsic :: iso_c_binding + import + implicit none + + integer (c_int64_t) , intent(in) , value :: context + double precision, intent(out) :: mo_value(*) + integer (c_int64_t) , intent(in) , value :: size_max + end function qmckl_get_mo_basis_mo_value + end interface + #+end_src + + Uses the given array to compute the values. + + #+begin_src c :comments org :tangle (eval h_func) :noweb yes +qmckl_exit_code +qmckl_get_mo_basis_mo_value_inplace (qmckl_context context, + double* const mo_value, + const int64_t size_max); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code +qmckl_get_mo_basis_mo_value_inplace (qmckl_context context, + double* const mo_value, + const int64_t size_max) +{ + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return qmckl_failwith( context, + QMCKL_INVALID_CONTEXT, + "qmckl_get_mo_basis_mo_value", + NULL); + } + + qmckl_exit_code rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; + assert (ctx != NULL); + + const int64_t sze = ctx->mo_basis.mo_num * ctx->point.num; + if (size_max < sze) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_get_mo_basis_mo_value", + "input array too small"); + } + + rc = qmckl_context_touch(context); + if (rc != QMCKL_SUCCESS) return rc; + + double* old_array = ctx->mo_basis.mo_value; + + ctx->mo_basis.mo_value = mo_value; + + rc = qmckl_provide_mo_value(context); + if (rc != QMCKL_SUCCESS) return rc; + + ctx->mo_basis.mo_value = old_array; + + return QMCKL_SUCCESS; +} + #+end_src + + #+begin_src f90 :tangle (eval fh_func) :comments org :exports none + interface + integer(c_int32_t) function qmckl_get_mo_basis_mo_value_inplace (context, & + mo_value, size_max) bind(C) + use, intrinsic :: iso_c_binding + import + implicit none + integer (c_int64_t) , intent(in) , value :: context + double precision, intent(out) :: mo_value(*) + integer (c_int64_t) , intent(in) , value :: size_max + end function qmckl_get_mo_basis_mo_value_inplace + end interface + #+end_src + +*** Provide + + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none +qmckl_exit_code qmckl_provide_mo_value(qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_provide_mo_value(qmckl_context context) +{ + + qmckl_exit_code rc; + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; + assert (ctx != NULL); + + if (!ctx->ao_basis.provided) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_ao_basis", + NULL); + } + + rc = qmckl_provide_ao_value(context); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_ao_value", + NULL); + } + + if (!ctx->mo_basis.provided) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_mo_basis", + NULL); + } + + /* Compute if necessary */ + if (ctx->point.date > ctx->mo_basis.mo_value_date) { + + /* Allocate array */ + if (ctx->mo_basis.mo_value == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->point.num * ctx->mo_basis.mo_num * sizeof(double); + double* mo_value = (double*) qmckl_malloc(context, mem_info); + + if (mo_value == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_mo_basis_mo_value", + NULL); + } + ctx->mo_basis.mo_value = mo_value; + } + + if (ctx->mo_basis.mo_vgl_date == ctx->point.date) { + + // mo_vgl has been computed at this step: Just copy the data. + + double * v = &(ctx->mo_basis.mo_value[0]); + double * vgl = &(ctx->mo_basis.mo_vgl[0]); + for (int i=0 ; ipoint.num ; ++i) { + for (int k=0 ; kmo_basis.mo_num ; ++k) { + v[k] = vgl[k]; + } + v += ctx->mo_basis.mo_num; + vgl += ctx->mo_basis.mo_num * 5; + } + + } else { + + rc = qmckl_compute_mo_basis_mo_value(context, + ctx->ao_basis.ao_num, + ctx->mo_basis.mo_num, + ctx->point.num, + ctx->mo_basis.coefficient_t, + ctx->ao_basis.ao_value, + ctx->mo_basis.mo_value); + + if (rc != QMCKL_SUCCESS) { + return rc; + } + + } + + ctx->mo_basis.mo_value_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + #+end_src + +*** Compute + :PROPERTIES: + :Name: qmckl_compute_mo_basis_mo_value + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_mo_basis_mo_value_args + | Variable | Type | In/Out | Description | + |---------------------+-----------------------------+--------+-------------------------------------------------| + | ~context~ | ~qmckl_context~ | in | Global state | + | ~ao_num~ | ~int64_t~ | in | Number of AOs | + | ~mo_num~ | ~int64_t~ | in | Number of MOs | + | ~point_num~ | ~int64_t~ | in | Number of points | + | ~coef_normalized_t~ | ~double[mo_num][ao_num]~ | in | Transpose of the AO to MO transformation matrix | + | ~ao_value~ | ~double[point_num][ao_num]~ | in | Value of the AOs | + | ~mo_value~ | ~double[point_num][mo_num]~ | out | Value of the MOs | + + + The matrix of AO values is very sparse, so we use a sparse-dense + matrix multiplication instead of a dgemm, as exposed in + https://dx.doi.org/10.1007/978-3-642-38718-0_14. + + + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer function qmckl_compute_mo_basis_mo_value_doc_f(context, & + ao_num, mo_num, point_num, & + coef_normalized_t, ao_value, mo_value) & + result(info) + use qmckl + implicit none + integer(qmckl_context), intent(in) :: context + integer*8 , intent(in) :: ao_num, mo_num + integer*8 , intent(in) :: point_num + double precision , intent(in) :: ao_value(ao_num,point_num) + double precision , intent(in) :: coef_normalized_t(mo_num,ao_num) + double precision , intent(out) :: mo_value(mo_num,point_num) + integer*8 :: i,j,k + double precision :: c1, c2, c3, c4, c5 + + integer*8 :: LDA, LDB, LDC + + info = QMCKL_SUCCESS + if (.True.) then ! fast algorithm + do j=1,point_num + mo_value(:,j) = 0.d0 + do k=1,ao_num + if (ao_value(k,j) /= 0.d0) then + c1 = ao_value(k,j) + do i=1,mo_num + mo_value(i,j) = mo_value(i,j) + coef_normalized_t(i,k) * c1 + end do + end if + end do + end do + + else ! dgemm + + LDA = size(coef_normalized_t,1) + LDB = size(ao_value,1) + LDC = size(mo_value,1) + + info = qmckl_dgemm(context,'N', 'N', mo_num, point_num, ao_num, 1.d0, & + coef_normalized_t, LDA, ao_value, LDB, & + 0.d0, mo_value, LDC) + + end if + +end function qmckl_compute_mo_basis_mo_value_doc_f + #+end_src + + #+CALL: generate_c_header(table=qmckl_mo_basis_mo_value_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_value")) + + #+RESULTS: + #+begin_src c :tangle (eval h_func) :comments org + qmckl_exit_code qmckl_compute_mo_basis_mo_value ( + const qmckl_context context, + const int64_t ao_num, + const int64_t mo_num, + const int64_t point_num, + const double* coef_normalized_t, + const double* ao_value, + double* const mo_value ); + #+end_src + + #+CALL: generate_c_header(table=qmckl_mo_basis_mo_value_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_value_doc")) + + #+RESULTS: + #+begin_src c :tangle (eval h_func) :comments org + qmckl_exit_code qmckl_compute_mo_basis_mo_value_doc ( + const qmckl_context context, + const int64_t ao_num, + const int64_t mo_num, + const int64_t point_num, + const double* coef_normalized_t, + const double* ao_value, + double* const mo_value ); + #+end_src + + #+CALL: generate_c_interface(table=qmckl_mo_basis_mo_value_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_value_doc")) + + #+RESULTS: + #+begin_src f90 :tangle (eval f) :comments org :exports none + integer(c_int32_t) function qmckl_compute_mo_basis_mo_value_doc & + (context, ao_num, mo_num, point_num, coef_normalized_t, ao_value, mo_value) & + 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 :: ao_num + integer (c_int64_t) , intent(in) , value :: mo_num + integer (c_int64_t) , intent(in) , value :: point_num + real (c_double ) , intent(in) :: coef_normalized_t(ao_num,mo_num) + real (c_double ) , intent(in) :: ao_value(ao_num,point_num) + real (c_double ) , intent(out) :: mo_value(mo_num,point_num) + + integer(c_int32_t), external :: qmckl_compute_mo_basis_mo_value_doc_f + info = qmckl_compute_mo_basis_mo_value_doc_f & + (context, ao_num, mo_num, point_num, coef_normalized_t, ao_value, mo_value) + + end function qmckl_compute_mo_basis_mo_value_doc + #+end_src + + #+begin_src c :tangle (eval c) :comments org +qmckl_exit_code +qmckl_compute_mo_basis_mo_value (const qmckl_context context, + const int64_t ao_num, + const int64_t mo_num, + const int64_t point_num, + const double* coef_normalized_t, + const double* ao_value, + double* const mo_value ) +{ +#ifdef HAVE_HPC + return qmckl_compute_mo_basis_mo_value_hpc (context, ao_num, mo_num, point_num, coef_normalized_t, ao_value, mo_value); +#else + return qmckl_compute_mo_basis_mo_value_doc (context, ao_num, mo_num, point_num, coef_normalized_t, ao_value, mo_value); +#endif +} + #+end_src + +*** HPC version + + + #+begin_src c :tangle (eval h_func) :comments org +#ifdef HAVE_HPC +qmckl_exit_code +qmckl_compute_mo_basis_mo_value_hpc (const qmckl_context context, + const int64_t ao_num, + const int64_t mo_num, + const int64_t point_num, + const double* coef_normalized_t, + const double* ao_value, + double* const mo_value ); +#endif + #+end_src + + #+begin_src c :tangle (eval c) :comments org +#ifdef HAVE_HPC +qmckl_exit_code +qmckl_compute_mo_basis_mo_value_hpc (const qmckl_context context, + const int64_t ao_num, + const int64_t mo_num, + const int64_t point_num, + const double* restrict coef_normalized_t, + const double* restrict ao_value, + double* restrict const mo_value ) +{ + assert (context != QMCKL_NULL_CONTEXT); + +#ifdef HAVE_OPENMP + #pragma omp parallel for +#endif + for (int64_t ipoint=0 ; ipoint < point_num ; ++ipoint) { + double* restrict const vgl1 = &(mo_value[ipoint*mo_num]); + const double* restrict avgl1 = &(ao_value[ipoint*ao_num]); + + for (int64_t i=0 ; ipoint.num * 5 * ctx->mo_basis.mo_num; @@ -507,7 +965,7 @@ qmckl_get_mo_basis_mo_vgl_inplace (qmckl_context context, qmckl_exit_code rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); const int64_t sze = ctx->mo_basis.mo_num * 5 * ctx->point.num; @@ -563,7 +1021,7 @@ qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context) return QMCKL_NULL_CONTEXT; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); if (!ctx->ao_basis.provided) { @@ -666,25 +1124,41 @@ integer function qmckl_compute_mo_basis_mo_vgl_doc_f(context, & integer*8 :: i,j,k double precision :: c1, c2, c3, c4, c5 - do j=1,point_num - mo_vgl(:,:,j) = 0.d0 - do k=1,ao_num - if (ao_vgl(k,1,j) /= 0.d0) then - c1 = ao_vgl(k,1,j) - c2 = ao_vgl(k,2,j) - c3 = ao_vgl(k,3,j) - c4 = ao_vgl(k,4,j) - c5 = ao_vgl(k,5,j) - do i=1,mo_num - mo_vgl(i,1,j) = mo_vgl(i,1,j) + coef_normalized_t(i,k) * c1 - mo_vgl(i,2,j) = mo_vgl(i,2,j) + coef_normalized_t(i,k) * c2 - mo_vgl(i,3,j) = mo_vgl(i,3,j) + coef_normalized_t(i,k) * c3 - mo_vgl(i,4,j) = mo_vgl(i,4,j) + coef_normalized_t(i,k) * c4 - mo_vgl(i,5,j) = mo_vgl(i,5,j) + coef_normalized_t(i,k) * c5 - end do - end if + integer*8 :: LDA, LDB, LDC + + info = QMCKL_SUCCESS + if (.True.) then ! fast algorithm + do j=1,point_num + mo_vgl(:,:,j) = 0.d0 + do k=1,ao_num + if (ao_vgl(k,1,j) /= 0.d0) then + c1 = ao_vgl(k,1,j) + c2 = ao_vgl(k,2,j) + c3 = ao_vgl(k,3,j) + c4 = ao_vgl(k,4,j) + c5 = ao_vgl(k,5,j) + do i=1,mo_num + mo_vgl(i,1,j) = mo_vgl(i,1,j) + coef_normalized_t(i,k) * c1 + mo_vgl(i,2,j) = mo_vgl(i,2,j) + coef_normalized_t(i,k) * c2 + mo_vgl(i,3,j) = mo_vgl(i,3,j) + coef_normalized_t(i,k) * c3 + mo_vgl(i,4,j) = mo_vgl(i,4,j) + coef_normalized_t(i,k) * c4 + mo_vgl(i,5,j) = mo_vgl(i,5,j) + coef_normalized_t(i,k) * c5 + end do + end if + end do end do - end do + + else ! dgemm + + LDA = size(coef_normalized_t,1) + LDB = size(ao_vgl,1) + LDC = size(mo_vgl,1) + + info = qmckl_dgemm(context,'N', 'N', mo_num, point_num*5_8, ao_num*1_8, 1.d0, & + coef_normalized_t, LDA, ao_vgl, LDB, & + 0.d0, mo_vgl, LDC) + + end if end function qmckl_compute_mo_basis_mo_vgl_doc_f #+end_src @@ -761,7 +1235,6 @@ qmckl_compute_mo_basis_mo_vgl (const qmckl_context context, } #+end_src - *** HPC version @@ -789,16 +1262,19 @@ qmckl_compute_mo_basis_mo_vgl_hpc (const qmckl_context context, const double* restrict ao_vgl, double* restrict const mo_vgl ) { + assert (context != QMCKL_NULL_CONTEXT); + #ifdef HAVE_OPENMP #pragma omp parallel for #endif for (int64_t ipoint=0 ; ipoint < point_num ; ++ipoint) { double* restrict const vgl1 = &(mo_vgl[ipoint*5*mo_num]); - const double* restrict avgl1 = &(ao_vgl[ipoint*5*ao_num]); double* restrict const vgl2 = vgl1 + mo_num; double* restrict const vgl3 = vgl1 + (mo_num << 1); double* restrict const vgl4 = vgl1 + (mo_num << 1) + mo_num; double* restrict const vgl5 = vgl1 + (mo_num << 2); + + const double* restrict avgl1 = &(ao_vgl[ipoint*5*ao_num]); const double* restrict avgl2 = avgl1 + ao_num; const double* restrict avgl3 = avgl1 + (ao_num << 1); const double* restrict avgl4 = avgl1 + (ao_num << 1) + ao_num; @@ -820,7 +1296,6 @@ qmckl_compute_mo_basis_mo_vgl_hpc (const qmckl_context context, double av4[ao_num]; double av5[ao_num]; for (int64_t k=0 ; knucleus.uninitialized = (1 << 3) - 1; @@ -167,7 +167,7 @@ qmckl_get_nucleus_num (const qmckl_context context, int64_t* const num) { "num is a null pointer"); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int32_t mask = 1 << 0; @@ -226,7 +226,7 @@ qmckl_get_nucleus_charge (const qmckl_context context, "charge is a null pointer"); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int32_t mask = 1 << 1; @@ -293,7 +293,7 @@ qmckl_get_nucleus_rescale_factor (const qmckl_context context, "rescale_factor_kappa is a null pointer"); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); assert (ctx->nucleus.rescale_factor_kappa > 0.0); @@ -351,7 +351,7 @@ qmckl_get_nucleus_coord (const qmckl_context context, "coord is a null pointer"); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int32_t mask = 1 << 2; @@ -410,7 +410,7 @@ bool qmckl_nucleus_provided(const qmckl_context context) { return false; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); return ctx->nucleus.provided; @@ -425,7 +425,7 @@ if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } -qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; +qmckl_context_struct* const ctx = (qmckl_context_struct*) context; #+end_src #+NAME:post2 @@ -672,7 +672,6 @@ end interface ** Test #+begin_src c :tangle (eval c_test) -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 = 2.0; @@ -688,13 +687,13 @@ rc = qmckl_get_nucleus_num (context, &n); assert(rc == QMCKL_NOT_PROVIDED); -rc = qmckl_set_nucleus_num (context, nucl_num); +rc = qmckl_set_nucleus_num (context, chbrclf_nucl_num); assert(rc == QMCKL_SUCCESS); assert(!qmckl_nucleus_provided(context)); rc = qmckl_get_nucleus_num (context, &n); assert(rc == QMCKL_SUCCESS); -assert(n == nucl_num); +assert(n == chbrclf_nucl_num); double k; rc = qmckl_get_nucleus_rescale_factor (context, &k); @@ -709,41 +708,41 @@ rc = qmckl_get_nucleus_rescale_factor (context, &k); assert(rc == QMCKL_SUCCESS); assert(k == nucl_rescale_factor_kappa); -double nucl_coord2[3*nucl_num]; +double nucl_coord2[3*chbrclf_nucl_num]; -rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2, 3*nucl_num); +rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2, 3*chbrclf_nucl_num); assert(rc == QMCKL_NOT_PROVIDED); -rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), 3*nucl_num); +rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), 3*chbrclf_nucl_num); assert(rc == QMCKL_SUCCESS); assert(!qmckl_nucleus_provided(context)); -rc = qmckl_get_nucleus_coord (context, 'N', nucl_coord2, 3*nucl_num); +rc = qmckl_get_nucleus_coord (context, 'N', nucl_coord2, 3*chbrclf_nucl_num); assert(rc == QMCKL_SUCCESS); for (size_t k=0 ; k<3 ; ++k) { - for (int64_t i=0 ; inucleus.num * ctx->nucleus.num; @@ -828,7 +827,7 @@ qmckl_exit_code qmckl_provide_nn_distance(qmckl_context context) return (char) 0; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); if (!ctx->nucleus.provided) return QMCKL_NOT_PROVIDED; @@ -940,10 +939,10 @@ qmckl_exit_code qmckl_compute_nn_distance ( assert(qmckl_nucleus_provided(context)); -double distance[nucl_num*nucl_num]; -rc = qmckl_get_nucleus_nn_distance(context, distance, nucl_num*nucl_num); +double distance[chbrclf_nucl_num*chbrclf_nucl_num]; +rc = qmckl_get_nucleus_nn_distance(context, distance, chbrclf_nucl_num*chbrclf_nucl_num); assert(distance[0] == 0.); -assert(distance[1] == distance[nucl_num]); +assert(distance[1] == distance[chbrclf_nucl_num]); assert(fabs(distance[1]-2.070304721365169) < 1.e-12); #+end_src @@ -973,7 +972,7 @@ qmckl_get_nucleus_nn_distance_rescaled(qmckl_context context, qmckl_exit_code rc = qmckl_provide_nn_distance_rescaled(context); if (rc != QMCKL_SUCCESS) return rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); const int64_t sze = ctx->nucleus.num * ctx->nucleus.num; @@ -1019,7 +1018,7 @@ qmckl_exit_code qmckl_provide_nn_distance_rescaled(qmckl_context context) return (char) 0; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); if (!ctx->nucleus.provided) return QMCKL_NOT_PROVIDED; @@ -1167,7 +1166,7 @@ qmckl_exit_code qmckl_get_nucleus_repulsion(qmckl_context context, double* const qmckl_exit_code rc = qmckl_provide_nucleus_repulsion(context); if (rc != QMCKL_SUCCESS) return rc; - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); *energy = ctx->nucleus.repulsion; @@ -1203,7 +1202,7 @@ qmckl_exit_code qmckl_provide_nucleus_repulsion(qmckl_context context) return (char) 0; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); qmckl_exit_code rc; diff --git a/org/qmckl_numprec.org b/org/qmckl_numprec.org index 7c5f26f..3174464 100644 --- a/org/qmckl_numprec.org +++ b/org/qmckl_numprec.org @@ -141,7 +141,7 @@ qmckl_exit_code qmckl_set_numprec_precision(const qmckl_context context, const i "precision > 53"); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; /* This should be always true because the context is valid */ assert (ctx != NULL); @@ -185,7 +185,7 @@ int qmckl_get_numprec_precision(const qmckl_context context) { ""); } - const qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + const qmckl_context_struct* const ctx = (qmckl_context_struct*) context; return ctx->numprec.precision; } #+end_src @@ -232,7 +232,7 @@ qmckl_exit_code qmckl_set_numprec_range(const qmckl_context context, const int r "range > 11"); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; /* This should be always true because the context is valid */ assert (ctx != NULL); @@ -275,7 +275,7 @@ int qmckl_get_numprec_range(const qmckl_context context) { ""); } - const qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + const qmckl_context_struct* const ctx = (qmckl_context_struct*) context; return ctx->numprec.range; } #+end_src diff --git a/org/qmckl_point.org b/org/qmckl_point.org index 07b0863..74a1e6a 100644 --- a/org/qmckl_point.org +++ b/org/qmckl_point.org @@ -81,7 +81,7 @@ int main() { |----------+----------------+-------------------------------------------| | ~num~ | ~int64_t~ | Total number of points | | ~date~ | ~uint64_t~ | Last modification date of the coordinates | - | ~coord~ | ~qmckl_matrix~ | ~num~ \times 3 matrix3 | + | ~coord~ | ~qmckl_matrix~ | ~num~ \times 3 matrix | We consider that the matrix is stored 'transposed' and 'normal' corresponds to the 3 \times ~num~ matrix. @@ -108,7 +108,7 @@ qmckl_exit_code qmckl_init_point(qmckl_context context) { return false; } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); memset(&(ctx->point), 0, sizeof(qmckl_point_struct)); @@ -148,7 +148,7 @@ qmckl_get_point_num (const qmckl_context context, int64_t* const num) { "num is a null pointer"); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); assert (ctx->point.num >= (int64_t) 0); @@ -202,7 +202,7 @@ qmckl_get_point(const qmckl_context context, "coord is a null pointer"); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int64_t point_num = ctx->point.num; @@ -295,7 +295,7 @@ qmckl_set_point (qmckl_context context, "coord is a NULL pointer"); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); qmckl_exit_code rc; diff --git a/org/qmckl_sherman_morrison_woodbury.org b/org/qmckl_sherman_morrison_woodbury.org index 598ad32..ae358e8 100644 --- a/org/qmckl_sherman_morrison_woodbury.org +++ b/org/qmckl_sherman_morrison_woodbury.org @@ -965,7 +965,7 @@ qmckl_exit_code qmckl_sherman_morrison_smw32s(const qmckl_context context, rc = qmckl_woodbury_3(context, LDS, Dim, Updates_3block, Updates_index_3block, breakdown, Slater_inv, determinant); if (rc != 0) { // Send the entire block to slagel_splitting uint64_t l = 0; - rc = qmckl_slagel_splitting(LDS, Dim, 3, Updates_3block, Updates_index_3block, + (void) qmckl_slagel_splitting(LDS, Dim, 3, Updates_3block, Updates_index_3block, breakdown, Slater_inv, later_updates + (Dim * later), later_index + later, &l, determinant); later = later + l; }