diff --git a/.github/workflows/test-build.yml b/.github/workflows/test-build.yml index a4dd421..41b59e0 100644 --- a/.github/workflows/test-build.yml +++ b/.github/workflows/test-build.yml @@ -40,7 +40,9 @@ jobs: - name: Build QMCkl run: | ./autogen.sh - ./configure --enable-silent-rules --enable-debug + mkdir _build + cd _build + ../configure --enable-silent-rules --enable-debug make -j 4 - name: Run test @@ -72,7 +74,7 @@ jobs: - uses: actions/checkout@v2 - name: install dependencies run: brew install emacs hdf5 automake pkg-config - + - name: Symlink gfortran (macOS) if: runner.os == 'macOS' run: | @@ -89,6 +91,28 @@ jobs: git clone https://github.com/TREX-CoE/trexio.git cd trexio ./autogen.sh + mkdir _build + cd _build + ../configure --enable-silent-rules + make -j 4 + + - name: Run test + run: make -j 4 check + + - name: Archive test log file + if: failure() + uses: actions/upload-artifact@v2 + with: + name: test-report-ubuntu + path: test-suite.log + + - name: Dist test + run: make distcheck + + - name: Archive test log file + if: failure() + uses: actions/upload-artifact@v2 + with: ./configure --prefix=${PWD}/_install --enable-silent-rules make -j 4 make install diff --git a/Makefile.am b/Makefile.am index 1137766..8eff2f1 100644 --- a/Makefile.am +++ b/Makefile.am @@ -82,18 +82,17 @@ htmldir_loc = share/doc/qmckl/html textdir_loc = share/doc/qmckl/text htmlize_el = $(htmldir_loc)/htmlize.el dist_html_DATA = $(HTML_FILES) \ - share/doc/qmckl/html/qmckl.css \ - share/doc/qmckl/html/index.html + $(htmldir_loc)/qmckl.css \ + $(htmldir_loc)/index.html dist_text_DATA = $(TEXT_FILES) -html-local: $(dist_html_DATA) -text: $(dist_text_DATA) +html-local: $(htmlize_el) $(dist_html_DATA) +text: $(htmlize_el) $(dist_text_DATA) doc: html text - if QMCKL_DEVEL if VFC_CI @@ -102,9 +101,9 @@ endif EXTRA_DIST += $(ORG_FILES) $(TANGLED_FILES) $(EXPORTED_FILES) $(test_qmckl_f) -BUILT_SOURCES = $(C_FILES) $(F_FILES) $(FH_FUNC_FILES) $(FH_TYPE_FILES) $(H_FUNC_FILES) $(H_TYPE_FILES) $(H_PRIVATE_FUNC_FILES) $(H_PRIVATE_TYPE_FILES) $(src_qmckl_f) $(test_qmckl_f) $(qmckl_h) $(header_tests) +BUILT_SOURCES = $(C_FILES) $(F_FILES) $(FH_FUNC_FILES) $(FH_TYPE_FILES) $(H_FUNC_FILES) $(H_TYPE_FILES) $(H_PRIVATE_FUNC_FILES) $(H_PRIVATE_TYPE_FILES) $(src_qmckl_f) $(test_qmckl_f) $(qmckl_h) $(header_tests) $(htmlize_el) -CLEANFILES += $(BUILT_SOURCES) $(C_TEST_FILES) $(F_TEST_FILES) $(TANGLED_FILES) $(C_TEST_FILES) $(F_TEST_FILES) $(src_qmckl_f) $(test_qmckl_f) $(qmckl_h) $(HTML_FILES) $(TEXT_FILES) $(EXPORTED_FILES) $(header_tests) +CLEANFILES += $(BUILT_SOURCES) $(C_TEST_FILES) $(F_TEST_FILES) $(TANGLED_FILES) $(C_TEST_FILES) $(F_TEST_FILES) $(src_qmckl_f) $(test_qmckl_f) $(qmckl_h) $(HTML_FILES) $(TEXT_FILES) $(EXPORTED_FILES) $(header_tests) $(htmlize_el) EXTRA_DIST += \ tools/build_doc.sh \ @@ -152,8 +151,7 @@ $(src_qmckl_f): $(FH_FUNC_FILES) $(FH_TYPE_FILES) $(cat_h_verbose)top_builddir=$(abs_top_builddir) srcdir=$(abs_srcdir) src_qmckl_f=$(src_qmckl_f) $(srcdir)/tools/build_qmckl_f.sh $(htmlize_el): - $(MKDIR_P) $(htmldir_loc) - $(MKDIR_P) $(textdir_loc) + $(MKDIR_P) $(htmldir_loc) $(textdir_loc) abs_srcdir=$(abs_srcdir) $(srcdir)/tools/install_htmlize.sh $(htmlize_el) tests/chbrclf.h: $(qmckl_h) diff --git a/README.md b/README.md index fb1281b..4ab5bae 100644 --- a/README.md +++ b/README.md @@ -36,12 +36,39 @@ make check ## For users -Obtain a source distribution and run +Obtain a source distribution. + +To build the documentation version: ``` -./configure -make -make check +./configure +``` + +To build an optimized version with Intel compilers: +``` +./configure \ + --with-icc \ + --with-ifort \ + --enable-hpc \ + --with-openmp +``` + +To build an optimized version with GCC: +``` +./configure \ + CC=gcc \ + CFLAGS="-g -O2 -march=native -flto -fno-trapping-math -fno-math-errno -ftree-vectorize" \ + FC=gfortran \ + FCFLAGS="-g -O2 -march=native -flto -ftree-vectorize" \ + --enable-hpc \ + --with-openmp +``` + + +Then, compile with: +``` +make -j +make -j check sudo make install sudo make installcheck ``` @@ -54,7 +81,12 @@ by the preprocessor otherwise). To enable vfc_ci support, the library should be configured with the following command : ``` -./configure --prefix=$PWD/_install \ --enable-vfc_ci --host=x86_64 \ CC="verificarlo-f" FC="verificarlo-f" +./configure \ + CC="verificarlo-f" \ + FC="verificarlo-f" \ + --prefix=$PWD/_install \ + --enable-vfc_ci \ + --host=x86_64 \ ``` where CC and FC are set to verificarlo-f, and support is explicitely enabled diff --git a/configure.ac b/configure.ac index 6654fd7..fcc1dc5 100644 --- a/configure.ac +++ b/configure.ac @@ -46,6 +46,22 @@ AS_IF([test -d ${srcdir}/.git], [enable_maintainer_mode="no"] ) +AC_ARG_WITH(ifort, [AS_HELP_STRING([--with-ifort],[Use Intel Fortran compiler])], with_ifort=$withval, with_ifort=no) +AS_IF([test "$with_ifort" == "yes"], [ + FC=ifort + FCFLAGS="-xHost -ip -O2 -ftz -finline -g -mkl=sequential" ]) + +AC_ARG_WITH(icc, [AS_HELP_STRING([--with-icc],[Use Intel C compiler])], with_icc=$withval, with_icc=no) +AS_IF([test "$with_icc" == "yes"], [ + CC=icc + CFLAGS="-xHost -ip -O2 -ftz -finline -g -mkl=sequential" ]) + +AS_IF([test "$with_icc"."$with_ifort" == "yes.yes"], [ + ax_blas_ok="yes" + ax_lapack_ok="yes" + BLAS_LIBS="" + LAPACK_LIBS=""]) + AM_PROG_AR AM_MAINTAINER_MODE() LT_INIT @@ -100,13 +116,14 @@ AC_CHECK_HEADERS([assert.h errno.h math.h pthread.h stdbool.h stdint.h stdio.h s AC_CHECK_LIB([pthread], [pthread_create]) # OpenMP -#AC_ARG_WITH(openmp, [AS_HELP_STRING([--with-openmp],[enable OpenMP])], with_omp=$withval, with_omp=no) -#if test "x$with_omp" = xyes; then -# AC_DEFINE([HAVE_OPENMP], [1], [Define to use OpenMP threading.]) -# AX_OPENMP([], -# [AC_MSG_ERROR([Could not find OpenMP flags; configure with --without-openmp])]) -# CFLAGS="${CFLAGS} ${OPENMP_CFLAGS}" -#fi +AC_ARG_WITH(openmp, [AS_HELP_STRING([--with-openmp],[activate OpenMP])], with_omp=$withval, with_omp=no) +if test "x$with_omp" = xyes; then + AC_DEFINE([HAVE_OPENMP], [1], [Define to use OpenMP threading.]) + AX_OPENMP([], + [AC_MSG_ERROR([Could not find OpenMP flags; configure with --without-openmp])]) + CFLAGS="${CFLAGS} ${OPENMP_CFLAGS}" + FCFLAGS="${CFLAGS} ${OPENMP_FCFLAGS}" +fi # CHAMELEON AC_ARG_WITH(chameleon, @@ -177,6 +194,7 @@ AX_BLAS([], [AC_MSG_ERROR([BLAS was not found.])]) ## LAPACK AX_LAPACK([], [AC_MSG_ERROR([LAPACK was not found.])]) +AS_IF([test "$BLAS_LIBS" == "$LAPACK_LIBS"], [BLAS_LIBS=""]) # Specific options required with some compilers @@ -186,10 +204,10 @@ case $FC in FCFLAGS="$FCFLAGS -nofor-main" ;; - gfortran*) - # Order is important here - FCFLAGS="-cpp $FCFLAGS" - ;; +#gfortran*) +# # Order is important here +# FCFLAGS="-cpp $FCFLAGS" +# ;; esac @@ -204,12 +222,13 @@ AC_ARG_ENABLE(debug, [AS_HELP_STRING([--enable-debug],[compile for debugging])], if test "$ok" = "yes"; then if test "$GCC" = "yes"; then CFLAGS="$CFLAGS \ --Wall -W -Wbad-function-cast -Wcast-qual \ --Wpointer-arith -Wcast-align -Wpedantic -Wextra -fmax-errors=3" +-g -Wall -W -Wbad-function-cast -Wcast-qual -Warray-bounds -Wdisabled-optimization \ +-Wpointer-arith -Wcast-align -Wpedantic -Wextra -Walloc-zero -Werror \ +" fi if test "$GFC" = "yes"; then FCFLAGS="$FCFLAGS \ --fcheck=all -Waliasing -Wampersand -Wconversion \ +-g -fcheck=all -Waliasing -Wampersand -Wconversion \ -Wsurprising -ffpe-trap=zero,overflow,underflow \ -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation \ -Wreal-q-constant -Wuninitialized -fbacktrace -finit-real=nan" @@ -314,7 +333,7 @@ fi #mkl-dynamic-lp64-seq PKG_LIBS="$PKG_LIBS $LIBS" -LIBS="$BLAS_LIBS $LAPACK_LIBS $BLAS_LIBS $PKG_LIBS" +LIBS="$BLAS_LIBS $LAPACK_LIBS $PKG_LIBS" CFLAGS="$CFLAGS $PKG_CFLAGS" AC_SUBST([PKG_LIBS]) AC_SUBST([PKG_CFLAGS]) diff --git a/org/qmckl_ao.org b/org/qmckl_ao.org index d40e478..55057bd 100644 --- a/org/qmckl_ao.org +++ b/org/qmckl_ao.org @@ -73,6 +73,8 @@ gradients and Laplacian of the atomic basis functions. #include #include +#include + #include "chbrclf.h" #include "qmckl_ao_private_func.h" @@ -96,6 +98,8 @@ int main() { #include #include #include +#include +#include #include "qmckl.h" #include "qmckl_context_private_type.h" @@ -2397,6 +2401,7 @@ assert (rc == QMCKL_SUCCESS); for (int64_t i=0 ; i < shell_num ; ++i) { assert(shell_ang_mom_test[i] == shell_ang_mom[i]); } +free(shell_ang_mom_test); shell_factor_test = (double*) malloc ( shell_num * sizeof(double)); rc = qmckl_get_ao_basis_shell_factor (context, shell_factor_test, shell_num); @@ -2405,6 +2410,7 @@ assert (rc == QMCKL_SUCCESS); for (int64_t i=0 ; i < shell_num ; ++i) { assert(shell_factor_test[i] == shell_factor[i]); } +free(shell_factor_test); shell_prim_num_test = (int64_t*) malloc ( shell_num * sizeof(int64_t)); rc = qmckl_get_ao_basis_shell_prim_num (context, shell_prim_num_test, shell_num); @@ -2413,6 +2419,7 @@ assert (rc == QMCKL_SUCCESS); for (int64_t i=0 ; i < shell_num ; ++i) { assert(shell_prim_num_test[i] == shell_prim_num[i]); } +free(shell_prim_num_test); shell_prim_index_test = (int64_t*) malloc ( shell_num * sizeof(int64_t)); rc = qmckl_get_ao_basis_shell_prim_index (context, shell_prim_index_test, shell_num); @@ -2421,6 +2428,7 @@ assert (rc == QMCKL_SUCCESS); for (int64_t i=0 ; i < shell_num ; ++i) { assert(shell_prim_index_test[i] == shell_prim_index[i]); } +free(shell_prim_index_test); exponent_test = (double*) malloc ( prim_num * sizeof(double)); rc = qmckl_get_ao_basis_exponent(context, exponent_test, prim_num); @@ -2429,6 +2437,7 @@ assert (rc == QMCKL_SUCCESS); for (int64_t i=0 ; i < prim_num ; ++i) { assert(exponent_test[i] == exponent[i]); } +free(exponent_test); coefficient_test = (double*) malloc ( prim_num * sizeof(double)); rc = qmckl_get_ao_basis_coefficient(context, coefficient_test, prim_num); @@ -2437,6 +2446,7 @@ assert (rc == QMCKL_SUCCESS); for (int64_t i=0 ; i < prim_num ; ++i) { assert(coefficient_test[i] == coefficient[i]); } +free(coefficient_test); prim_factor_test = (double*) malloc ( prim_num * sizeof(double)); rc = qmckl_get_ao_basis_prim_factor (context, prim_factor_test, prim_num); @@ -2445,6 +2455,7 @@ assert (rc == QMCKL_SUCCESS); for (int64_t i=0 ; i < prim_num ; ++i) { assert(prim_factor_test[i] == prim_factor[i]); } +free(prim_factor_test); rc = qmckl_get_ao_basis_ao_num(context, &ao_num_test); assert(ao_num == ao_num_test); @@ -2456,6 +2467,7 @@ assert (rc == QMCKL_SUCCESS); for (int64_t i=0 ; i < ao_num ; ++i) { assert(ao_factor_test[i] == ao_factor[i]); } +free(ao_factor_test); #+end_src @@ -2466,11 +2478,11 @@ for (int64_t i=0 ; i < ao_num ; ++i) { |----------------------+-----------------------------------+----------------------------------------------------------------------------------------------| | Variable | Type | Description | |----------------------+-----------------------------------+----------------------------------------------------------------------------------------------| - | ~primitive_vgl~ | ~double[point_num][5][prim_num]~ | Value, gradients, Laplacian of the primitives at current positions | + | ~primitive_vgl~ | ~double[point_num][5][prim_num]~ | Value, gradients, Laplacian of the primitives at current positions | | ~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 primitives at current positions | | ~ao_vgl_date~ | ~uint64_t~ | Last modification date of Value, gradients, Laplacian of the AOs at current positions | @@ -2810,6 +2822,72 @@ qmckl_get_ao_basis_ao_vgl (qmckl_context context, end interface #+end_src + Uses the give array to compute the VGL. + + #+begin_src c :comments org :tangle (eval h_func) :noweb yes +qmckl_exit_code +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, + double* const ao_vgl, + 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_vgl", + NULL); + } + + qmckl_exit_code rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + int64_t sze = ctx->ao_basis.ao_num * 5 * ctx->point.num; + if (size_max < sze) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_get_ao_basis_ao_vgl", + "input array too small"); + } + + rc = qmckl_context_touch(context); + if (rc != QMCKL_SUCCESS) return rc; + + double* old_array = ctx->ao_basis.ao_vgl; + + ctx->ao_basis.ao_vgl = ao_vgl; + + rc = qmckl_provide_ao_vgl(context); + if (rc != QMCKL_SUCCESS) return rc; + + ctx->ao_basis.ao_vgl = 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_vgl_inplace (context, & + ao_vgl, 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_vgl(*) + integer (c_int64_t) , intent(in) , value :: size_max + end function qmckl_get_ao_basis_ao_vgl_inplace + end interface + #+end_src + * Radial part ** General functions for Gaussian basis functions @@ -3076,7 +3154,7 @@ assert(0 == test_qmckl_ao_gaussian_vgl(context)); const double* coord, const double* nucl_coord, const double* expo, - double* const primitive_vgl ); + double* const primitive_vgl ); #+end_src @@ -3387,6 +3465,7 @@ for (j=0 ; j cutoff*nucleus_range(inucl)) then + cycle + end if + + ! C is zero-based, so shift bounds by one + ishell_start = nucleus_index(inucl) + 1 + ishell_end = nucleus_index(inucl) + nucleus_shell_num(inucl) + do ishell=ishell_start, ishell_end shell_vgl(ishell, 1, ipoint) = 0.d0 @@ -3523,6 +3608,7 @@ end function qmckl_compute_ao_basis_shell_gaussian_vgl_f nucl_num, & nucleus_shell_num, & nucleus_index, & + nucleus_range, & shell_prim_index, & shell_prim_num, & coord, & @@ -3542,6 +3628,7 @@ end function qmckl_compute_ao_basis_shell_gaussian_vgl_f integer (c_int64_t) , intent(in) , value :: nucl_num integer (c_int64_t) , intent(in) :: nucleus_shell_num(nucl_num) integer (c_int64_t) , intent(in) :: nucleus_index(nucl_num) + real (c_double ) , intent(in) :: nucleus_range(nucl_num) integer (c_int64_t) , intent(in) :: shell_prim_index(shell_num) integer (c_int64_t) , intent(in) :: shell_prim_num(shell_num) real (c_double ) , intent(in) :: coord(point_num,3) @@ -3559,6 +3646,7 @@ end function qmckl_compute_ao_basis_shell_gaussian_vgl_f nucl_num, & nucleus_shell_num, & nucleus_index, & + nucleus_range, & shell_prim_index, & shell_prim_num, & coord, & @@ -3625,6 +3713,7 @@ qmckl_exit_code qmckl_provide_ao_basis_shell_vgl(qmckl_context context) ctx->nucleus.num, ctx->ao_basis.nucleus_shell_num, ctx->ao_basis.nucleus_index, + ctx->ao_basis.nucleus_range, ctx->ao_basis.shell_prim_index, ctx->ao_basis.shell_prim_num, ctx->point.coord.data, @@ -4047,8 +4136,45 @@ assert(0 == test_qmckl_ao_power(context)); const int64_t ldv ); #+end_src + #+CALL: generate_c_header(table=qmckl_ao_polynomial_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_ao_polynomial_vgl_doc") + + #+RESULTS: + #+begin_src c :tangle (eval h_func) :comments org + qmckl_exit_code qmckl_ao_polynomial_vgl_doc ( + const qmckl_context context, + const double* X, + const double* R, + const int32_t lmax, + int64_t* n, + int32_t* const L, + const int64_t ldl, + double* const VGL, + const int64_t ldv ); + #+end_src + + #+begin_src c :tangle (eval c) :comments org +qmckl_exit_code +qmckl_ao_polynomial_vgl (const qmckl_context context, + const double* X, + const double* R, + const int32_t lmax, + int64_t* n, + int32_t* const L, + const int64_t ldl, + double* const VGL, + const int64_t ldv ) +{ +#ifdef HAVE_HPC + //return qmckl_ao_polynomial_vgl_hpc (context, X, R, lmax, n, L, ldl, VGL, ldv); + return qmckl_ao_polynomial_vgl_doc (context, X, R, lmax, n, L, ldl, VGL, ldv); +#else + return qmckl_ao_polynomial_vgl_doc (context, X, R, lmax, n, L, ldl, VGL, ldv); +#endif +} + #+end_src + #+begin_src f90 :tangle (eval f) -integer function qmckl_ao_polynomial_vgl_f (context, & +integer function qmckl_ao_polynomial_vgl_doc_f (context, & X, R, lmax, n, L, ldl, VGL, ldv) result(info) use qmckl implicit none @@ -4064,7 +4190,6 @@ integer function qmckl_ao_polynomial_vgl_f (context, & integer*8 :: i,j integer :: a,b,c,d real*8 :: Y(3) - integer :: lmax_array(3) real*8 :: pows(-2:lmax,3) double precision :: xy, yz, xz double precision :: da, db, dc, dd @@ -4096,7 +4221,6 @@ integer function qmckl_ao_polynomial_vgl_f (context, & Y(i) = X(i) - R(i) end do - lmax_array(1:3) = lmax if (lmax == 0) then VGL(1,1) = 1.d0 VGL(2:5,1) = 0.d0 @@ -4174,14 +4298,14 @@ integer function qmckl_ao_polynomial_vgl_f (context, & info = QMCKL_SUCCESS -end function qmckl_ao_polynomial_vgl_f +end function qmckl_ao_polynomial_vgl_doc_f #+end_src - #+CALL: generate_c_interface(table=qmckl_ao_polynomial_vgl_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + #+CALL: generate_c_interface(table=qmckl_ao_polynomial_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_ao_polynomial_vgl_doc" ) #+RESULTS: #+begin_src f90 :tangle (eval f) :comments org :exports none - integer(c_int32_t) function qmckl_ao_polynomial_vgl & + integer(c_int32_t) function qmckl_ao_polynomial_vgl_doc & (context, X, R, lmax, n, L, ldl, VGL, ldv) & bind(C) result(info) @@ -4198,14 +4322,40 @@ end function qmckl_ao_polynomial_vgl_f real (c_double ) , intent(out) :: VGL(ldv,n) integer (c_int64_t) , intent(in) , value :: ldv - integer(c_int32_t), external :: qmckl_ao_polynomial_vgl_f - info = qmckl_ao_polynomial_vgl_f & + integer(c_int32_t), external :: qmckl_ao_polynomial_vgl_doc_f + info = qmckl_ao_polynomial_vgl_doc_f & (context, X, R, lmax, n, L, ldl, VGL, ldv) - end function qmckl_ao_polynomial_vgl + end function qmckl_ao_polynomial_vgl_doc #+end_src - #+CALL: generate_f_interface(table=qmckl_ao_polynomial_vgl_args,rettyp=get_value("FRetType"),fname=get_value("Name")) + #+CALL: generate_f_interface(table=qmckl_ao_polynomial_vgl_args,rettyp=get_value("FRetType"),fname="qmckl_ao_polynomial_vgl_doc" ) + + #+RESULTS: + #+begin_src f90 :tangle (eval fh_func) :comments org :exports none + interface + integer(c_int32_t) function qmckl_ao_polynomial_vgl_doc & + (context, X, R, lmax, n, L, ldl, VGL, ldv) & + bind(C) + use, intrinsic :: iso_c_binding + import + implicit none + + integer (c_int64_t) , intent(in) , value :: context + real (c_double ) , intent(in) :: X(3) + real (c_double ) , intent(in) :: R(3) + integer (c_int32_t) , intent(in) , value :: lmax + integer (c_int64_t) , intent(inout) :: n + integer (c_int32_t) , intent(out) :: L(ldl,n) + integer (c_int64_t) , intent(in) , value :: ldl + real (c_double ) , intent(out) :: VGL(ldv,n) + integer (c_int64_t) , intent(in) , value :: ldv + + end function qmckl_ao_polynomial_vgl_doc + end interface + #+end_src + + #+CALL: generate_f_interface(table=qmckl_ao_polynomial_vgl_args,rettyp=get_value("FRetType"),fname="qmckl_ao_polynomial_vgl" ) #+RESULTS: #+begin_src f90 :tangle (eval fh_func) :comments org :exports none @@ -4231,7 +4381,9 @@ end function qmckl_ao_polynomial_vgl_f end interface #+end_src + #+CALL: generate_c_header(table=qmckl_ao_polynomial_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_ao_polynomial_transp_vgl") + #+RESULTS: #+begin_src c :tangle (eval h_func) :comments org qmckl_exit_code qmckl_ao_polynomial_transp_vgl ( const qmckl_context context, @@ -4245,8 +4397,60 @@ end function qmckl_ao_polynomial_vgl_f const int64_t ldv ); #+end_src + #+CALL: generate_c_header(table=qmckl_ao_polynomial_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_ao_polynomial_transp_vgl_doc") + + #+RESULTS: + #+begin_src c :tangle (eval h_func) :comments org + qmckl_exit_code qmckl_ao_polynomial_transp_vgl_doc ( + const qmckl_context context, + const double* X, + const double* R, + const int32_t lmax, + int64_t* n, + int32_t* const L, + const int64_t ldl, + double* const VGL, + const int64_t ldv ); + #+end_src + + #+CALL: generate_c_header(table=qmckl_ao_polynomial_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_ao_polynomial_transp_vgl_hpc") + + #+RESULTS: + #+begin_src c :tangle (eval h_func) :comments org + qmckl_exit_code qmckl_ao_polynomial_transp_vgl_hpc ( + const qmckl_context context, + const double* X, + const double* R, + const int32_t lmax, + int64_t* n, + int32_t* const L, + const int64_t ldl, + double* const VGL, + const int64_t ldv ); + #+end_src + + #+begin_src c :tangle (eval c) :comments org +qmckl_exit_code +qmckl_ao_polynomial_transp_vgl (const qmckl_context context, + const double* X, + const double* R, + const int32_t lmax, + int64_t* n, + int32_t* const L, + const int64_t ldl, + double* const VGL, + const int64_t ldv ) +{ +#ifdef HAVE_HPC + return qmckl_ao_polynomial_transp_vgl_hpc (context, X, R, lmax, n, L, ldl, VGL, ldv); +#else + return qmckl_ao_polynomial_transp_vgl_doc (context, X, R, lmax, n, L, ldl, VGL, ldv); +#endif +} + #+end_src + #+begin_src f90 :tangle (eval f) -integer function qmckl_ao_polynomial_transp_vgl_f (context, & +integer function qmckl_ao_polynomial_transp_vgl_doc_f (context, & X, R, lmax, n, L, ldl, VGL, ldv) result(info) use qmckl implicit none @@ -4262,8 +4466,7 @@ integer function qmckl_ao_polynomial_transp_vgl_f (context, & integer*8 :: i,j integer :: a,b,c,d real*8 :: Y(3) - integer :: lmax_array(3) - real*8 :: pows(-2:lmax,3) + real*8 :: pows(-2:21,3) ! lmax < 22 double precision :: xy, yz, xz double precision :: da, db, dc, dd @@ -4290,17 +4493,11 @@ integer function qmckl_ao_polynomial_transp_vgl_f (context, & endif - do i=1,3 - Y(i) = X(i) - R(i) - end do + if (lmax > 0) then - lmax_array(1:3) = lmax - if (lmax == 0) then - VGL(1,1) = 1.d0 - VGL(1,2:5) = 0.d0 - l(1:3,1) = 0 - n=1 - else if (lmax > 0) then + do i=1,3 + Y(i) = X(i) - R(i) + end do pows(-2:0,1:3) = 1.d0 do i=1,lmax pows(i,1) = pows(i-1,1) * Y(1) @@ -4312,21 +4509,26 @@ integer function qmckl_ao_polynomial_transp_vgl_f (context, & VGL(1:4,1:5) = 0.d0 VGL(1 ,1 ) = 1.d0 - VGL(2:4,1:5) = 0.d0 l (1,2) = 1 - VGL(2,1) = pows(1,1) + VGL(2,1) = Y(1) VGL(2,2) = 1.d0 l (2,3) = 1 - VGL(3,1) = pows(1,2) + VGL(3,1) = Y(2) VGL(3,3) = 1.d0 l (3,4) = 1 - VGL(4,1) = pows(1,3) + VGL(4,1) = Y(3) VGL(4,4) = 1.d0 n=4 + else + VGL(1,1) = 1.d0 + VGL(1,2:5) = 0.d0 + l(1:3,1) = 0 + n=1 + return endif ! l>=2 @@ -4340,14 +4542,14 @@ integer function qmckl_ao_polynomial_transp_vgl_f (context, & dc = dd - da - db n = n+1 - l(1,n) = a - l(2,n) = b - l(3,n) = c - xy = pows(a,1) * pows(b,2) yz = pows(b,2) * pows(c,3) xz = pows(a,1) * pows(c,3) + l(1,n) = a + l(2,n) = b + l(3,n) = c + VGL(n,1) = xy * pows(c,3) xy = dc * xy @@ -4372,14 +4574,292 @@ integer function qmckl_ao_polynomial_transp_vgl_f (context, & info = QMCKL_SUCCESS -end function qmckl_ao_polynomial_transp_vgl_f +end function qmckl_ao_polynomial_transp_vgl_doc_f #+end_src - #+CALL: generate_c_interface(table=qmckl_ao_polynomial_vgl_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + #+begin_src c :tangle (eval c) :comments org +qmckl_exit_code +qmckl_ao_polynomial_transp_vgl_hpc (const qmckl_context context, + const double* X, + const double* R, + const int32_t lmax, + int64_t* n, + int32_t* const L, + const int64_t ldl, + double* const VGL, + const int64_t ldv ) +{ + const qmckl_context_struct* ctx = (qmckl_context_struct* const) 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; + + int32_t m; + + switch (lmax) { + case 0: + { + if (ldv < 1) return QMCKL_INVALID_ARG_9; + L[0] = 0; L[1] = 0; L[2] = 0; + + VGL[0 ] = 1.0; + VGL[ldv ] = 0.0; + VGL[ldv<<1 ] = 0.0; + VGL[(ldv<<1)+ldv] = 0.0; + VGL[ldv<<2 ] = 0.0; + m=1; + break; + } + case 1: + { + if (ldv < 4) return QMCKL_INVALID_ARG_9; + + if (ldl == 3) { + + const int32_t ltmp[12] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1}; + for (int i=0 ; i<12 ; ++i) + L[i] = ltmp[i]; + + } else { + + int32_t* const l[4] = {L, L+ldl, L+(ldl<<1), L+ldl+(ldl<<1)}; + l[0][0] = 0; l[0][1] = 0; l[0][2] = 0; + l[1][0] = 1; l[1][1] = 0; l[1][2] = 0; + l[2][0] = 0; l[2][1] = 1; l[2][2] = 0; + l[3][0] = 0; l[3][1] = 0; l[3][2] = 1; + } + + + if (ldv == 4) { + + const double tmp[20] = {1.0, X[0]-R[0], X[1]-R[1], X[2]-R[2], + 0.0, 1.0, 0.0, 0.0, + 0.0, 0.0, 1.0, 0.0, + 0.0, 0.0, 0.0, 1.0, + 0.0, 0.0, 0.0, 0.0}; + + for (int i=0 ; i<20 ; ++i) + VGL[i] = tmp[i]; + + } else { + + double* const vgl1 = VGL; + double* const vgl2 = VGL + ldv; + double* const vgl3 = VGL + (ldv << 1); + double* const vgl4 = VGL + ldv + (ldv << 1); + double* const vgl5 = VGL + (ldv << 2); + + for (int32_t k=0 ; k<4 ; ++k) { + vgl2[k] = 0.0; + vgl3[k] = 0.0; + vgl4[k] = 0.0; + vgl5[k] = 0.0; + } + vgl1[0] = 1.0; + vgl1[1] = X[0]-R[0]; + vgl1[2] = X[1]-R[1]; + vgl1[3] = X[2]-R[2]; + vgl2[1] = 1.0; + vgl3[2] = 1.0; + vgl4[3] = 1.0; + } + m=4; + break; + } + case 2: + { + if (ldv < 10) return QMCKL_INVALID_ARG_9; + if (ldl == 3) { + const int32_t ltmp[30] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, + 2, 0, 0, 1, 1, 0, 1, 0, 1, 0, 2, 0, + 0, 1, 1, 0, 0, 2}; + for (int i=0 ; i<30 ; ++i) + L[i] = ltmp[i]; + + } else { + int32_t* l[10]; + for (int32_t i=0 ; i<10 ; ++i) { + l[i] = L + i*ldl; + } + + l[0][0] = 0; l[0][1] = 0; l[0][2] = 0; + l[1][0] = 1; l[1][1] = 0; l[1][2] = 0; + l[2][0] = 0; l[2][1] = 1; l[2][2] = 0; + l[3][0] = 0; l[3][1] = 0; l[3][2] = 1; + l[4][0] = 2; l[4][1] = 0; l[4][2] = 0; + l[5][0] = 1; l[5][1] = 1; l[5][2] = 0; + l[6][0] = 1; l[6][1] = 0; l[6][2] = 1; + l[7][0] = 0; l[7][1] = 2; l[7][2] = 0; + l[8][0] = 0; l[8][1] = 1; l[8][2] = 1; + l[9][0] = 0; l[9][1] = 0; l[9][2] = 2; + } + + const double Y[3] = { X[0]-R[0], + X[1]-R[1], + X[2]-R[2] }; + + if (ldv == 50) { + const double tmp[50] = { + 1.0, Y[0], Y[1], Y[2], Y[0] * Y[0], + Y[0] * Y[1], Y[0] * Y[2], Y[1] * Y[1], Y[1] * Y[2], Y[2] * Y[2], + 0.0, 1.0, 0.0, 0.0, Y[0] + Y[0], + Y[1], Y[2], 0.0, 0.0, 0.0, + 0.0, 0.0, 1.0, 0.0, 0.0, + Y[0], 0.0, Y[1] + Y[1], Y[2], 0.0, + 0.0, 0.0, 0.0, 1.0, 0.0, + 0.0, Y[0], 0.0, Y[1], Y[2] + Y[2], + 0.0, 0.0, 0.0, 0.0, 2.0, + 0.0, 0.0, 2.0, 0., 2.0 }; + for (int i=0 ; i<50 ; ++i) + VGL[i] = tmp[i]; + } else { + double* const vgl1 = VGL; + double* const vgl2 = VGL + ldv; + double* const vgl3 = VGL + (ldv << 1); + double* const vgl4 = VGL + 3*ldv; + double* const vgl5 = VGL + (ldv << 2); + + vgl1[0] = 1.0 ; vgl1[1] = Y[0] ; vgl1[2] = Y[1]; + vgl1[3] = Y[2] ; vgl1[4] = Y[0]*Y[0]; vgl1[5] = Y[0]*Y[1]; + vgl1[6] = Y[0]*Y[2]; vgl1[7] = Y[1]*Y[1]; vgl1[8] = Y[1]*Y[2]; + vgl1[9] = Y[2]*Y[2]; + + vgl2[0] = 0.0 ; vgl2[1] = 1.0 ; vgl2[2] = 0.0 ; + vgl2[3] = 0.0 ; vgl2[4] = Y[0]+Y[0]; vgl2[5] = Y[1]; + vgl2[6] = Y[2]; vgl2[7] = 0.0 ; vgl2[8] = 0.0 ; + vgl2[9] = 0.0 ; + + vgl3[0] = 0.0; vgl3[1] = 0.0 ; vgl3[2] = 1.0 ; + vgl3[3] = 0.0; vgl3[4] = 0.0 ; vgl3[5] = Y[0]; + vgl3[6] = 0.0; vgl3[7] = Y[1]+Y[1]; vgl3[8] = Y[2]; + vgl3[9] = 0.0; + + vgl4[0] = 0.0 ; vgl4[1] = 0.0; vgl4[2] = 0.0 ; + vgl4[3] = 1.0 ; vgl4[4] = 0.0; vgl4[5] = 0.0 ; + vgl4[6] = Y[0] ; vgl4[7] = 0.0; vgl4[8] = Y[1]; + vgl4[9] = Y[2]+Y[2]; + + vgl5[0] = 0.0; vgl5[1] = 0.0; vgl5[2] = 0.0; + vgl5[3] = 0.0; vgl5[4] = 2.0; vgl5[5] = 0.0; + vgl5[6] = 0.0; vgl5[7] = 2.0; vgl5[8] = 0.0; + vgl5[9] = 2.0; + } + m=10; + break; + } + default: + { + const int32_t size_max = (lmax+1)*(lmax+2)*(lmax+3)/6; + if (ldv < size_max) return QMCKL_INVALID_ARG_9; + + double* const vgl1 = VGL; + double* const vgl2 = VGL + ldv; + double* const vgl3 = VGL + (ldv<<1); + double* const vgl4 = VGL + ldv + (ldv<<1); + double* const vgl5 = VGL + (ldv<<2); + + const double Y[3] = { X[0]-R[0], + X[1]-R[1], + X[2]-R[2] }; + + assert(size_max > lmax+3); + double pows[3][size_max]; + + for (int32_t i=0 ; i<3 ; ++i) { + pows[0][i] = 1.0; + pows[1][i] = 1.0; + pows[2][i] = 1.0; + } + + for (int32_t i=3 ; i<=lmax+2 ; ++i) { + pows[0][i] = pows[0][i-1] * Y[0]; + pows[1][i] = pows[1][i-1] * Y[1]; + pows[2][i] = pows[2][i-1] * Y[2]; + } + + int32_t* l[size_max]; + for (int32_t i=0 ; i=0 ; --a) { + double db = dd-da; + + for (int32_t b=d-a ; b>=0 ; --b) { + const int32_t c = d - a - b; + const double dc = dd - da - db; + + double xy = pows[0][a+2] * pows[1][b+2]; + double yz = pows[1][b+2] * pows[2][c+2]; + double xz = pows[0][a+2] * pows[2][c+2]; + + l[m][0] = a; + l[m][1] = b; + l[m][2] = c; + + vgl1[m] = xy * pows[2][c+2]; + + xy *= dc; + xz *= db; + yz *= da; + + vgl2[m] = pows[0][a+1] * yz; + vgl3[m] = pows[1][b+1] * xz; + vgl4[m] = pows[2][c+1] * xy; + + vgl5[m] = (da-1.) * pows[0][a] * yz + + (db-1.) * pows[1][b] * xz + + (dc-1.) * pows[2][c] * xy; + + db -= 1.0; + ++m; + } + da -= 1.0; + } + dd += 1.0; + } + } + } + ,*n = m; + return QMCKL_SUCCESS; +} + #+end_src + + #+CALL: generate_c_interface(table=qmckl_ao_polynomial_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_ao_polynomial_transp_vgl_doc") #+RESULTS: #+begin_src f90 :tangle (eval f) :comments org :exports none - integer(c_int32_t) function qmckl_ao_polynomial_transp_vgl & + integer(c_int32_t) function qmckl_ao_polynomial_transp_vgl_doc & (context, X, R, lmax, n, L, ldl, VGL, ldv) & bind(C) result(info) @@ -4393,17 +4873,43 @@ end function qmckl_ao_polynomial_transp_vgl_f integer (c_int64_t) , intent(inout) :: n integer (c_int32_t) , intent(out) :: L(ldl,n) integer (c_int64_t) , intent(in) , value :: ldl - real (c_double ) , intent(out) :: VGL(ldv,5) + real (c_double ) , intent(out) :: VGL(ldv,n) integer (c_int64_t) , intent(in) , value :: ldv - integer(c_int32_t), external :: qmckl_ao_polynomial_transp_vgl_f - info = qmckl_ao_polynomial_transp_vgl_f & + integer(c_int32_t), external :: qmckl_ao_polynomial_transp_vgl_doc_f + info = qmckl_ao_polynomial_transp_vgl_doc_f & (context, X, R, lmax, n, L, ldl, VGL, ldv) - end function qmckl_ao_polynomial_transp_vgl + end function qmckl_ao_polynomial_transp_vgl_doc #+end_src - #+CALL: generate_f_interface(table=qmckl_ao_polynomial_vgl_args,rettyp=get_value("FRetType"),fname=get_value("Name")) + #+CALL: generate_f_interface(table=qmckl_ao_polynomial_vgl_args,rettyp=get_value("FRetType"),fname="qmckl_ao_polynomial_transp_vgl_doc") + + #+RESULTS: + #+begin_src f90 :tangle (eval fh_func) :comments org :exports none + interface + integer(c_int32_t) function qmckl_ao_polynomial_transp_vgl_doc & + (context, X, R, lmax, n, L, ldl, VGL, ldv) & + bind(C) + use, intrinsic :: iso_c_binding + import + implicit none + + integer (c_int64_t) , intent(in) , value :: context + real (c_double ) , intent(in) :: X(3) + real (c_double ) , intent(in) :: R(3) + integer (c_int32_t) , intent(in) , value :: lmax + integer (c_int64_t) , intent(inout) :: n + integer (c_int32_t) , intent(out) :: L(ldl,n) + integer (c_int64_t) , intent(in) , value :: ldl + real (c_double ) , intent(out) :: VGL(ldv,n) + integer (c_int64_t) , intent(in) , value :: ldv + + end function qmckl_ao_polynomial_transp_vgl_doc + end interface + #+end_src + + #+CALL: generate_f_interface(table=qmckl_ao_polynomial_vgl_args,rettyp=get_value("FRetType"),fname="qmckl_ao_polynomial_transp_vgl") #+RESULTS: #+begin_src f90 :tangle (eval fh_func) :comments org :exports none @@ -4422,11 +4928,11 @@ end function qmckl_ao_polynomial_transp_vgl_f integer (c_int64_t) , intent(inout) :: n integer (c_int32_t) , intent(out) :: L(ldl,n) integer (c_int64_t) , intent(in) , value :: ldl - real (c_double ) , intent(out) :: VGL(ldv,5) + real (c_double ) , intent(out) :: VGL(ldv,n) integer (c_int64_t) , intent(in) , value :: ldv end function qmckl_ao_polynomial_transp_vgl - end interface + end interface #+end_src *** Test :noexport: @@ -4524,8 +5030,44 @@ end function test_qmckl_ao_polynomial_vgl #+end_src #+begin_src c :tangle (eval c_test) - int test_qmckl_ao_polynomial_vgl(qmckl_context context); - assert(0 == test_qmckl_ao_polynomial_vgl(context)); +int test_qmckl_ao_polynomial_vgl(qmckl_context context); +assert(0 == test_qmckl_ao_polynomial_vgl(context)); + +double X[3] = { 1.1, 2.2, 3.3 }; +double R[3] = { 0.2, 1.1, 3.0 }; +int32_t ldv[8] = {1, 4, 10, 20, 35, 56, 84, 120}; +for (int32_t ldl=3 ; ldl<=5 ; ++ldl) { + int64_t n; + int32_t L0[200][ldl]; + int32_t L1[200][ldl]; + printf("ldl=%d\n", ldl); + for (int32_t lmax=0 ; lmax<=7 ; lmax++) { + double VGL0[5][ldv[lmax]]; + double VGL1[5][ldv[lmax]]; + memset(&L0[0][0], 0, sizeof(L0)); + memset(&L1[0][0], 0, sizeof(L1)); + memset(&VGL0[0][0], 0, sizeof(VGL0)); + memset(&VGL1[0][0], 0, sizeof(VGL1)); + rc = qmckl_ao_polynomial_transp_vgl_doc (context, X, R, lmax, &n, &(L0[0][0]), ldl, &(VGL0[0][0]), ldv[lmax]); + assert (rc == QMCKL_SUCCESS); + rc = qmckl_ao_polynomial_transp_vgl_hpc (context, X, R, lmax, &n, &(L1[0][0]), ldl, &(VGL1[0][0]), ldv[lmax]); + assert (rc == QMCKL_SUCCESS); + printf("lmax=%d\n", lmax); + for (int32_t l=0 ; l cutoff*nucleus_range(inucl)) then cycle end if - ! Compute polynomials - info = qmckl_ao_polynomial_vgl_f(context, e_coord, n_coord, & + ! 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) @@ -4635,6 +5186,7 @@ integer function qmckl_compute_ao_vgl_doc_f(context, & 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 @@ -4680,147 +5232,302 @@ end function qmckl_compute_ao_vgl_doc_f #+end_src ** HPC version - #+begin_src f90 :comments org :tangle (eval f) :noweb yes -integer function qmckl_compute_ao_vgl_hpc_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_vgl) & - 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_vgl(ao_num,5,point_num) + #+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 | - 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_transp_vgl_f + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code +qmckl_compute_ao_vgl_hpc_gaussian ( + const qmckl_context context, + const int64_t ao_num, + const int64_t shell_num, + const int64_t prim_num, + const int64_t point_num, + const int64_t nucl_num, + const double* coord, + const double* nucl_coord, + const int64_t* nucleus_index, + const int64_t* nucleus_shell_num, + const double* nucleus_range, + const int32_t* nucleus_max_ang_mom, + const int32_t* shell_ang_mom, + const int64_t* shell_prim_index, + const int64_t* shell_prim_num, + const double* ao_factor, + const double* expo, + const double* coef_normalized, + double* const ao_vgl ) +{ + int32_t lstart[32]; + for (int32_t l=0 ; l<32 ; ++l) { + lstart[l] = l*(l+1)*(l+2)/6; + } - double precision, allocatable :: poly_vgl(:,:) - integer , allocatable :: powers(:,:) + int64_t ao_index[shell_num+1]; + int64_t size_max = 0; + { + int64_t k=0; + for (int inucl=0 ; inucl < nucl_num ; ++inucl) { + 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; + } - allocate(poly_vgl(ao_num,5), powers(3,ao_num)) + /* Don't compute polynomials when the radial part is zero. */ + double cutoff = -log(1.e-12); - ! Pre-computed data - do l=0,20 - lstart(l) = l*(l+1)*(l+2)/6 +1 - end do +#ifdef HAVE_OPENMP + #pragma omp parallel +#endif + { + qmckl_exit_code rc; + double c_[prim_num]; + double expo_[prim_num]; + double ar2[prim_num]; + int32_t powers[prim_num]; + 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}, + {0.0, 0.0, 0.0, 1.0}}; + double poly_vgl_l2[5][10] = {{1., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, + {0., 1., 0., 0., 0., 0., 0., 0., 0., 0.}, + {0., 0., 1., 0., 0., 0., 0., 0., 0., 0.}, + {0., 0., 0., 1., 0., 0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 2., 0., 0., 2., 0., 2.}}; + double poly_vgl[5][size_max]; - info = QMCKL_SUCCESS +#ifdef HAVE_OPENMP + #pragma omp for +#endif + for (int64_t ipoint=0 ; ipoint < point_num ; ++ipoint) { + const double e_coord[3] = { coord[ipoint], + coord[ipoint + point_num], + coord[ipoint + 2*point_num] }; - ! Don't compute polynomials when the radial part is zero. - cutoff = -dlog(1.d-12) + for (int64_t inucl=0 ; inucl < nucl_num ; ++inucl) { + const double n_coord[3] = { nucl_coord[inucl], + nucl_coord[inucl + nucl_num], + nucl_coord[inucl + 2*nucl_num] }; - do ipoint = 1, point_num - e_coord(1) = coord(ipoint,1) - e_coord(2) = coord(ipoint,2) - e_coord(3) = coord(ipoint,3) - k=1 - 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 */ + const double x = e_coord[0] - n_coord[0]; + const double y = e_coord[1] - n_coord[1]; + const double z = e_coord[2] - n_coord[2]; - ! 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) + const double r2 = x*x + y*y + z*z; - r2 = x*x + z*z + z*z + if (r2 > cutoff * nucleus_range[inucl]) { + continue; + } - if (r2 > cutoff*nucleus_range(inucl)) then - cycle - end if + int64_t n_poly; + switch (nucleus_max_ang_mom[inucl]) { + case 0: + break; - ! Compute polynomials - info = qmckl_ao_polynomial_transp_vgl_f(context, e_coord, n_coord, & - nucleus_max_ang_mom(inucl), n_poly, powers, 3_8, & - poly_vgl, int(ao_num,8)) + case 1: + poly_vgl_l1[0][1] = x; + poly_vgl_l1[0][2] = y; + poly_vgl_l1[0][3] = z; + break; - ! Loop over shells - 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) - if (shell_vgl(ishell,1,ipoint) /= 0.d0) then - do il = lstart(l), lstart(l+1)-1 - ! Value - ao_vgl(k,1,ipoint) = & - poly_vgl(il,1) * shell_vgl(ishell,1,ipoint) * ao_factor(k) - - ! Grad_x - ao_vgl(k,2,ipoint) = ( & - poly_vgl(il,2) * shell_vgl(ishell,1,ipoint) + & - poly_vgl(il,1) * shell_vgl(ishell,2,ipoint) & - ) * ao_factor(k) - - ! Grad_y - ao_vgl(k,3,ipoint) = ( & - poly_vgl(il,3) * shell_vgl(ishell,1,ipoint) + & - poly_vgl(il,1) * shell_vgl(ishell,3,ipoint) & - ) * ao_factor(k) - - ! Grad_z - ao_vgl(k,4,ipoint) = ( & - poly_vgl(il,4) * shell_vgl(ishell,1,ipoint) + & - poly_vgl(il,1) * shell_vgl(ishell,4,ipoint) & - ) * ao_factor(k) - - ! Lapl_z - ao_vgl(k,5,ipoint) = ( & - poly_vgl(il,5) * shell_vgl(ishell,1,ipoint) + & - poly_vgl(il,1) * shell_vgl(ishell,5,ipoint) + & - 2.d0 * ( & - poly_vgl(il,2) * shell_vgl(ishell,2,ipoint) + & - poly_vgl(il,3) * shell_vgl(ishell,3,ipoint) + & - poly_vgl(il,4) * shell_vgl(ishell,4,ipoint) ) & - ) * ao_factor(k) - k = k+1 - end do - else - do il = lstart(l), lstart(l+1)-1 - ao_vgl(k,1,ipoint) = 0.d0 - ao_vgl(k,2,ipoint) = 0.d0 - ao_vgl(k,3,ipoint) = 0.d0 - ao_vgl(k,4,ipoint) = 0.d0 - ao_vgl(k,5,ipoint) = 0.d0 - k = k+1 - end do - end if - end do - end do - end do + case 2: + poly_vgl_l2[0][1] = x; + poly_vgl_l2[0][2] = y; + poly_vgl_l2[0][3] = z; + poly_vgl_l2[0][4] = x*x; + poly_vgl_l2[0][5] = x*y; + poly_vgl_l2[0][6] = x*z; + poly_vgl_l2[0][7] = y*y; + poly_vgl_l2[0][8] = y*z; + poly_vgl_l2[0][9] = z*z; + poly_vgl_l2[1][4] = x+x; + poly_vgl_l2[1][5] = y; + poly_vgl_l2[1][6] = z; + poly_vgl_l2[2][5] = x; + poly_vgl_l2[2][7] = y+y; + poly_vgl_l2[2][8] = z; + poly_vgl_l2[3][6] = x; + poly_vgl_l2[3][8] = y; + poly_vgl_l2[3][9] = z+z; + break; - deallocate(poly_vgl, powers) -end function qmckl_compute_ao_vgl_hpc_f + 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[0][0]), size_max); + + assert (rc == QMCKL_SUCCESS); + break; + } + + 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) { + + int64_t nidx = 0; + const int64_t iprim_start = shell_prim_index[ishell]; + const int64_t iprim_end = shell_prim_index[ishell] + shell_prim_num[ishell]; + for (int64_t iprim = iprim_start ; iprim < iprim_end ; ++iprim) { + const double v = expo[iprim]*r2; + if (v <= cutoff) { + ar2[nidx] = v; + c_[nidx] = coef_normalized[iprim]; + expo_[nidx] = expo[iprim]; + ++nidx; + } + } + + double s1_ = 0.; + double s5_ = 0.; + double s6_ = 0.; + for (int idx=0 ; idx 0) { + const double* f = &(ao_factor[k]); + const int64_t idx = lstart[l]; + + switch (nucleus_max_ang_mom[inucl]) { + case 0: + ao_vgl_1[0] = s1 * f[0]; + ao_vgl_2[0] = s2 * f[0]; + ao_vgl_3[0] = s3 * f[0]; + ao_vgl_4[0] = s4 * f[0]; + ao_vgl_5[0] = s5; + break; + case 1: + poly_vgl_1 = &(poly_vgl_l1[0][idx]); + poly_vgl_2 = &(poly_vgl_l1[1][idx]); + poly_vgl_3 = &(poly_vgl_l1[2][idx]); + poly_vgl_4 = &(poly_vgl_l1[3][idx]); + for (int64_t il=0 ; ilao_basis.ao_vgl == NULL) { @@ -4958,22 +5685,83 @@ qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context) } ctx->ao_basis.ao_vgl = ao_vgl; } - - rc = qmckl_compute_ao_vgl(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_vgl); +#ifdef HAVE_HPC + if (ctx->ao_basis.type == 'G') { + rc = qmckl_compute_ao_vgl_hpc_gaussian(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_vgl); + /* + } else if (ctx->ao_basis.type == 'S') { + rc = qmck_compute_ao_vgl_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_vgl); + ,*/ + } else { + rc = qmckl_compute_ao_vgl_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_vgl); + } +#else + rc = qmckl_compute_ao_vgl_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_vgl); +#endif if (rc != QMCKL_SUCCESS) { return rc; } diff --git a/org/qmckl_blas.org b/org/qmckl_blas.org index 998e6ef..274d269 100644 --- a/org/qmckl_blas.org +++ b/org/qmckl_blas.org @@ -111,26 +111,26 @@ qmckl_vector_alloc( qmckl_context context, #+begin_src c :comments org :tangle (eval h_private_func) qmckl_exit_code qmckl_vector_free( qmckl_context context, - qmckl_vector vector); + qmckl_vector* vector); #+end_src #+begin_src c :comments org :tangle (eval c) qmckl_exit_code qmckl_vector_free( qmckl_context context, - qmckl_vector vector) + qmckl_vector* vector) { /* Always true */ - assert (vector.data != NULL); + assert (vector->data != NULL); qmckl_exit_code rc; - rc = qmckl_free(context, vector.data); + rc = qmckl_free(context, vector->data); if (rc != QMCKL_SUCCESS) { return rc; } - vector.size = (int64_t) 0; - vector.data = NULL; + vector->size = (int64_t) 0; + vector->data = NULL; return QMCKL_SUCCESS; } #+end_src @@ -192,26 +192,26 @@ qmckl_matrix_alloc( qmckl_context context, #+begin_src c :comments org :tangle (eval h_private_func) qmckl_exit_code qmckl_matrix_free( qmckl_context context, - qmckl_matrix matrix); + qmckl_matrix* matrix); #+end_src #+begin_src c :comments org :tangle (eval c) qmckl_exit_code qmckl_matrix_free( qmckl_context context, - qmckl_matrix matrix) + qmckl_matrix* matrix) { /* Always true */ - assert (matrix.data != NULL); + assert (matrix->data != NULL); qmckl_exit_code rc; - rc = qmckl_free(context, matrix.data); + rc = qmckl_free(context, matrix->data); if (rc != QMCKL_SUCCESS) { return rc; } - matrix.data = NULL; - matrix.size[0] = (int64_t) 0; - matrix.size[1] = (int64_t) 0; + matrix->data = NULL; + matrix->size[0] = (int64_t) 0; + matrix->size[1] = (int64_t) 0; return QMCKL_SUCCESS; } @@ -286,25 +286,25 @@ qmckl_tensor_alloc( qmckl_context context, #+begin_src c :comments org :tangle (eval h_private_func) qmckl_exit_code qmckl_tensor_free( qmckl_context context, - qmckl_tensor tensor); + qmckl_tensor* tensor); #+end_src #+begin_src c :comments org :tangle (eval c) qmckl_exit_code qmckl_tensor_free( qmckl_context context, - qmckl_tensor tensor) + qmckl_tensor* tensor) { /* Always true */ - assert (tensor.data != NULL); + assert (tensor->data != NULL); qmckl_exit_code rc; - rc = qmckl_free(context, tensor.data); + rc = qmckl_free(context, tensor->data); if (rc != QMCKL_SUCCESS) { return rc; } - memset(&tensor, 0, sizeof(qmckl_tensor)); + memset(tensor, 0, sizeof(qmckl_tensor)); return QMCKL_SUCCESS; } @@ -712,7 +712,7 @@ qmckl_tensor_of_double(const qmckl_context context, for (int64_t i=0 ; i

date += 1UL; + ctx->point.date += 1UL; + return QMCKL_SUCCESS; +} + #+end_src + ** Creation To create a new context, ~qmckl_context_create()~ should be used. diff --git a/org/qmckl_determinant.org b/org/qmckl_determinant.org index f7ca97b..a28f49c 100644 --- a/org/qmckl_determinant.org +++ b/org/qmckl_determinant.org @@ -316,7 +316,7 @@ int64_t* qmckl_get_determinant_mo_index_alpha (const qmckl_context context) { int32_t mask = 1 << 4; if ( (ctx->det.uninitialized & mask) != 0) { - return (int64_t) 0; + return NULL; } assert (ctx->det.mo_index_alpha != NULL); diff --git a/org/qmckl_electron.org b/org/qmckl_electron.org index 5e0d5dd..33ac366 100644 --- a/org/qmckl_electron.org +++ b/org/qmckl_electron.org @@ -499,18 +499,15 @@ ctx->electron.provided = (ctx->electron.uninitialized == 0); if (ctx->electron.provided) { if (ctx->electron.coord_new.data != NULL) { - const qmckl_exit_code rc = qmckl_matrix_free(context, ctx->electron.coord_new); + const qmckl_exit_code rc = qmckl_matrix_free(context, &(ctx->electron.coord_new)); assert (rc == QMCKL_SUCCESS); } if (ctx->electron.coord_old.data != NULL) { - const qmckl_exit_code rc = qmckl_matrix_free(context, ctx->electron.coord_old); + const qmckl_exit_code rc = qmckl_matrix_free(context, &(ctx->electron.coord_old)); assert (rc == QMCKL_SUCCESS); } - const int64_t walk_num = ctx->electron.walk_num; - const int64_t elec_num = ctx->electron.num; - - const int64_t point_num = walk_num * elec_num; + const int64_t point_num = ctx->electron.walk_num * ctx->electron.num; qmckl_matrix coord_new = qmckl_matrix_alloc(context, point_num, 3); @@ -1667,8 +1664,7 @@ qmckl_exit_code qmckl_provide_ee_potential(qmckl_context context) ctx->electron.ee_pot = ee_pot; } - qmckl_exit_code rc = - qmckl_compute_ee_potential(context, + rc = qmckl_compute_ee_potential(context, ctx->electron.num, ctx->electron.walk_num, ctx->electron.ee_distance, @@ -2712,8 +2708,7 @@ qmckl_exit_code qmckl_provide_en_potential(qmckl_context context) ctx->electron.en_pot = en_pot; } - qmckl_exit_code rc = - qmckl_compute_en_potential(context, + rc = qmckl_compute_en_potential(context, ctx->electron.num, ctx->nucleus.num, ctx->electron.walk_num, diff --git a/org/qmckl_jastrow.org b/org/qmckl_jastrow.org index 9867010..aac48b8 100644 --- a/org/qmckl_jastrow.org +++ b/org/qmckl_jastrow.org @@ -9,10 +9,10 @@ nuclear ($\mathbf{R}$) coordinates. Its defined as $\exp(J(\mathbf{r},\mathbf{R}))$, where \[ - J(\mathbf{r},\mathbf{R}) = J_{\text{eN}}(\mathbf{r},\mathbf{R}) + J_{\text{ee}}(\mathbf{r}) + J_{\text{eeN}}(\mathbf{r},\mathbf{R}) + 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 + 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|$. $J_{\text{eN}}$ contains electron-nucleus terms: @@ -23,9 +23,9 @@ \sum_{p=2}^{N_\text{ord}^a} a_{p+1}\, [f(R_{i\alpha})]^p - J_{eN}^\infty \] - $J_{\text{ee}}$ contains electron-electron terms: + $J_{\text{ee}}$ contains electron-electron terms: \[ - J_{\text{ee}}(\mathbf{r}) = + J_{\text{ee}}(\mathbf{r}) = \sum_{i=1}^{N_\text{elec}} \sum_{j=1}^{i-1} \frac{b_1\, f(r_{ij})}{1+b_2\, f(r_{ij})} + \sum_{p=2}^{N_\text{ord}^b} a_{p+1}\, [f(r_{ij})]^p - J_{ee}^\infty @@ -43,7 +43,7 @@ \sum_{l=0}^{p-k-2\delta_{k,0}} c_{lkp\alpha} \left[ g({r}_{ij}) \right]^k \left[ \left[ g({R}_{i\alpha}) \right]^l + \left[ g({R}_{j\alpha}) \right]^l \right] - \left[ g({R}_{i\,\alpha}) \, g({R}_{j\alpha}) \right]^{(p-k-l)/2} + \left[ g({R}_{i\,\alpha}) \, g({R}_{j\alpha}) \right]^{(p-k-l)/2} \] $c_{lkp\alpha}$ are non-zero only when $p-k-l$ is even. @@ -55,10 +55,10 @@ g(r) = e^{-\kappa\, r}. \] - The terms $J_{\text{ee}}^\infty$ and $J_{\text{eN}}^\infty$ are shifts to ensure that + 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") @@ -412,10 +412,10 @@ qmckl_exit_code qmckl_get_jastrow_aord_num (qmckl_context context, int qmckl_exit_code qmckl_get_jastrow_bord_num (qmckl_context context, int64_t* const bord_num); qmckl_exit_code qmckl_get_jastrow_cord_num (qmckl_context context, int64_t* const bord_num); qmckl_exit_code qmckl_get_jastrow_type_nucl_num (qmckl_context context, int64_t* const type_nucl_num); -qmckl_exit_code qmckl_get_jastrow_type_nucl_vector (qmckl_context context, int64_t* const type_nucl_num, int64_t* size_max); -qmckl_exit_code qmckl_get_jastrow_aord_vector (qmckl_context context, double * const aord_vector, int64_t* size_max); -qmckl_exit_code qmckl_get_jastrow_bord_vector (qmckl_context context, double * const bord_vector, int64_t* size_max); -qmckl_exit_code qmckl_get_jastrow_cord_vector (qmckl_context context, double * const cord_vector, int64_t* size_max); +qmckl_exit_code qmckl_get_jastrow_type_nucl_vector (qmckl_context context, int64_t* const type_nucl_num, const int64_t size_max); +qmckl_exit_code qmckl_get_jastrow_aord_vector (qmckl_context context, double * const aord_vector, const int64_t size_max); +qmckl_exit_code qmckl_get_jastrow_bord_vector (qmckl_context context, double * const bord_vector, const int64_t size_max); +qmckl_exit_code qmckl_get_jastrow_cord_vector (qmckl_context context, double * const cord_vector, const int64_t size_max); #+end_src Along with these core functions, calculation of the jastrow factor @@ -474,7 +474,7 @@ qmckl_exit_code qmckl_get_jastrow_aord_num (const qmckl_context context, int64_t } assert (ctx->jastrow.aord_num > 0); - *aord_num = ctx->jastrow.aord_num; + ,*aord_num = ctx->jastrow.aord_num; return QMCKL_SUCCESS; } @@ -501,7 +501,7 @@ qmckl_exit_code qmckl_get_jastrow_bord_num (const qmckl_context context, int64_t } assert (ctx->jastrow.bord_num > 0); - *bord_num = ctx->jastrow.bord_num; + ,*bord_num = ctx->jastrow.bord_num; return QMCKL_SUCCESS; } @@ -528,7 +528,7 @@ qmckl_exit_code qmckl_get_jastrow_cord_num (const qmckl_context context, int64_t } assert (ctx->jastrow.cord_num > 0); - *cord_num = ctx->jastrow.cord_num; + ,*cord_num = ctx->jastrow.cord_num; return QMCKL_SUCCESS; } @@ -555,11 +555,15 @@ qmckl_exit_code qmckl_get_jastrow_type_nucl_num (const qmckl_context context, in } assert (ctx->jastrow.type_nucl_num > 0); - *type_nucl_num = ctx->jastrow.type_nucl_num; + ,*type_nucl_num = ctx->jastrow.type_nucl_num; return QMCKL_SUCCESS; } -qmckl_exit_code qmckl_get_jastrow_type_nucl_vector (const qmckl_context context, int64_t * const type_nucl_vector, int64_t* size_max) { +qmckl_exit_code +qmckl_get_jastrow_type_nucl_vector (const qmckl_context context, + int64_t* const type_nucl_vector, + const int64_t size_max) +{ if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return (char) 0; @@ -582,12 +586,21 @@ qmckl_exit_code qmckl_get_jastrow_type_nucl_vector (const qmckl_context context, } assert (ctx->jastrow.type_nucl_vector != NULL); + if (size_max < ctx->jastrow.type_nucl_num) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_get_jastrow_type_nucl_vector", + "Array too small. Expected jastrow.type_nucl_num"); + } + memcpy(type_nucl_vector, ctx->jastrow.type_nucl_vector, ctx->jastrow.type_nucl_num*sizeof(int64_t)); - (*size_max) = ctx->jastrow.type_nucl_num; return QMCKL_SUCCESS; } -qmckl_exit_code qmckl_get_jastrow_aord_vector (const qmckl_context context, double * const aord_vector, int64_t* size_max) { +qmckl_exit_code +qmckl_get_jastrow_aord_vector (const qmckl_context context, + double * const aord_vector, + const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return (char) 0; @@ -611,12 +624,20 @@ qmckl_exit_code qmckl_get_jastrow_aord_vector (const qmckl_context context, doub assert (ctx->jastrow.aord_vector != NULL); int64_t sze = (ctx->jastrow.aord_num + 1)*ctx->jastrow.type_nucl_num; + if (size_max < sze) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_get_jastrow_aord_vector", + "Array too small. Expected (ctx->jastrow.aord_num + 1)*ctx->jastrow.type_nucl_num"); + } memcpy(aord_vector, ctx->jastrow.aord_vector, sze*sizeof(double)); - (*size_max) = sze; return QMCKL_SUCCESS; } -qmckl_exit_code qmckl_get_jastrow_bord_vector (const qmckl_context context, double * const bord_vector, int64_t* size_max) { +qmckl_exit_code +qmckl_get_jastrow_bord_vector (const qmckl_context context, + double * const bord_vector, + const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return (char) 0; @@ -640,12 +661,20 @@ qmckl_exit_code qmckl_get_jastrow_bord_vector (const qmckl_context context, doub assert (ctx->jastrow.bord_vector != NULL); int64_t sze=ctx->jastrow.bord_num +1; + if (size_max < sze) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_get_jastrow_bord_vector", + "Array too small. Expected (ctx->jastrow.bord_num + 1)"); + } memcpy(bord_vector, ctx->jastrow.bord_vector, sze*sizeof(double)); - (*size_max) = sze; return QMCKL_SUCCESS; } -qmckl_exit_code qmckl_get_jastrow_cord_vector (const qmckl_context context, double * const cord_vector, int64_t* size_max) { +qmckl_exit_code +qmckl_get_jastrow_cord_vector (const qmckl_context context, + double * const cord_vector, + const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return (char) 0; @@ -668,14 +697,19 @@ qmckl_exit_code qmckl_get_jastrow_cord_vector (const qmckl_context context, doub } assert (ctx->jastrow.cord_vector != NULL); - + int64_t dim_cord_vect; qmckl_exit_code rc = qmckl_get_jastrow_dim_cord_vect(context, &dim_cord_vect); if (rc != QMCKL_SUCCESS) return rc; int64_t sze=dim_cord_vect * ctx->jastrow.type_nucl_num; + if (size_max < sze) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_get_jastrow_cord_vector", + "Array too small. Expected dim_cord_vect * jastrow.type_nucl_num"); + } memcpy(cord_vector, ctx->jastrow.cord_vector, sze*sizeof(double)); - (*size_max) = sze; return QMCKL_SUCCESS; } @@ -690,9 +724,9 @@ qmckl_exit_code qmckl_get_jastrow_cord_vector (const qmckl_context context, doub qmckl_exit_code qmckl_set_jastrow_ord_num (qmckl_context context, const int64_t aord_num, const int64_t bord_num, const int64_t cord_num); qmckl_exit_code qmckl_set_jastrow_type_nucl_num (qmckl_context context, const int64_t type_nucl_num); qmckl_exit_code qmckl_set_jastrow_type_nucl_vector (qmckl_context context, const int64_t* type_nucl_vector, const int64_t nucl_num); -qmckl_exit_code qmckl_set_jastrow_aord_vector (qmckl_context context, const double * aord_vector, int64_t size_max); -qmckl_exit_code qmckl_set_jastrow_bord_vector (qmckl_context context, const double * bord_vector, int64_t size_max); -qmckl_exit_code qmckl_set_jastrow_cord_vector (qmckl_context context, const double * cord_vector, int64_t size_max); +qmckl_exit_code qmckl_set_jastrow_aord_vector (qmckl_context context, const double * aord_vector, const int64_t size_max); +qmckl_exit_code qmckl_set_jastrow_bord_vector (qmckl_context context, const double * bord_vector, const int64_t size_max); +qmckl_exit_code qmckl_set_jastrow_cord_vector (qmckl_context context, const double * cord_vector, const int64_t size_max); #+end_src #+NAME:pre2 @@ -718,7 +752,12 @@ return QMCKL_SUCCESS; #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -qmckl_exit_code qmckl_set_jastrow_ord_num(qmckl_context context, const int64_t aord_num, const int64_t bord_num, const int64_t cord_num) { +qmckl_exit_code +qmckl_set_jastrow_ord_num(qmckl_context context, + const int64_t aord_num, + const int64_t bord_num, + const int64_t cord_num) +{ <> if (aord_num <= 0) { @@ -750,7 +789,10 @@ qmckl_exit_code qmckl_set_jastrow_ord_num(qmckl_context context, const int64_t a <> } -qmckl_exit_code qmckl_set_jastrow_type_nucl_num(qmckl_context context, const int64_t type_nucl_num) { + +qmckl_exit_code +qmckl_set_jastrow_type_nucl_num(qmckl_context context, const int64_t type_nucl_num) +{ <> if (type_nucl_num <= 0) { @@ -766,7 +808,12 @@ qmckl_exit_code qmckl_set_jastrow_type_nucl_num(qmckl_context context, const int <> } -qmckl_exit_code qmckl_set_jastrow_type_nucl_vector(qmckl_context context, int64_t const * type_nucl_vector, const int64_t nucl_num) { + +qmckl_exit_code +qmckl_set_jastrow_type_nucl_vector(qmckl_context context, + int64_t const * type_nucl_vector, + const int64_t nucl_num) +{ <> int32_t mask = 1 << 2; @@ -816,7 +863,12 @@ qmckl_exit_code qmckl_set_jastrow_type_nucl_vector(qmckl_context context, int64_ <> } -qmckl_exit_code qmckl_set_jastrow_aord_vector(qmckl_context context, double const * aord_vector, int64_t size_max) { + +qmckl_exit_code +qmckl_set_jastrow_aord_vector(qmckl_context context, + double const * aord_vector, + const int64_t size_max) +{ <> int32_t mask = 1 << 3; @@ -849,13 +901,13 @@ qmckl_exit_code qmckl_set_jastrow_aord_vector(qmckl_context context, double cons return qmckl_failwith( context, rc, "qmckl_set_ord_vector", NULL); -} + } } qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = (aord_num + 1) * type_nucl_num * sizeof(double); - if (size_max < mem_info.size/sizeof(double)) { + if ((size_t) size_max < mem_info.size/sizeof(double)) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_set_jastrow_aord_vector", @@ -878,7 +930,12 @@ qmckl_exit_code qmckl_set_jastrow_aord_vector(qmckl_context context, double cons <> } -qmckl_exit_code qmckl_set_jastrow_bord_vector(qmckl_context context, double const * bord_vector, int64_t size_max) { + +qmckl_exit_code +qmckl_set_jastrow_bord_vector(qmckl_context context, + double const * bord_vector, + const int64_t size_max) +{ <> int32_t mask = 1 << 4; @@ -913,7 +970,7 @@ qmckl_exit_code qmckl_set_jastrow_bord_vector(qmckl_context context, double cons qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = (bord_num + 1) * sizeof(double); - if (size_max < mem_info.size/sizeof(double)) { + if ((size_t) size_max < mem_info.size/sizeof(double)) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_set_jastrow_bord_vector", @@ -936,7 +993,12 @@ qmckl_exit_code qmckl_set_jastrow_bord_vector(qmckl_context context, double cons <> } -qmckl_exit_code qmckl_set_jastrow_cord_vector(qmckl_context context, double const * cord_vector, int64_t size_max) { + +qmckl_exit_code +qmckl_set_jastrow_cord_vector(qmckl_context context, + double const * cord_vector, + const int64_t size_max) +{ <> int32_t mask = 1 << 5; @@ -978,7 +1040,7 @@ qmckl_exit_code qmckl_set_jastrow_cord_vector(qmckl_context context, double cons qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = dim_cord_vect * type_nucl_num * sizeof(double); - if (size_max < mem_info.size/sizeof(double)) { + if ((size_t) size_max < mem_info.size/sizeof(double)) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_set_jastrow_cord_vector", @@ -1069,6 +1131,7 @@ double* elec_coord = &(n2_elec_coord[0][0][0]); const double* nucl_charge = n2_charge; int64_t nucl_num = n2_nucl_num; double* nucl_coord = &(n2_nucl_coord[0][0]); +int64_t size_max; /* Provide Electron data */ @@ -1246,11 +1309,17 @@ assert(qmckl_nucleus_provided(context)); *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes -qmckl_exit_code qmckl_get_jastrow_asymp_jasb(qmckl_context context, double* const asymp_jasb, int64_t* size_max); +qmckl_exit_code +qmckl_get_jastrow_asymp_jasb(qmckl_context context, + double* const asymp_jasb, + const int64_t size_max); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -qmckl_exit_code qmckl_get_jastrow_asymp_jasb(qmckl_context context, double* const asymp_jasb, int64_t* size_max) +qmckl_exit_code +qmckl_get_jastrow_asymp_jasb(qmckl_context context, + double* const asymp_jasb, + const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; @@ -1264,9 +1333,14 @@ qmckl_exit_code qmckl_get_jastrow_asymp_jasb(qmckl_context context, double* cons qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - size_t sze = 2; + int64_t sze = 2; + if (size_max < sze) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_get_jastrow_asymp_jasb", + "Array too small. Expected 2"); + } memcpy(asymp_jasb, ctx->jastrow.asymp_jasb, sze * sizeof(double)); - (*size_max) = sze; return QMCKL_SUCCESS; } @@ -1409,7 +1483,7 @@ qmckl_exit_code qmckl_compute_asymp_jasb ( } asym_one = bord_vector[0] * kappa_inv / (1.0 + bord_vector[1] * kappa_inv); - asymp_jasb[0] = asym_one; + asymp_jasb[0] = asym_one; asymp_jasb[1] = 0.5 * asym_one; for (int i = 0 ; i <= 1; ++i) { @@ -1418,7 +1492,7 @@ qmckl_exit_code qmckl_compute_asymp_jasb ( x = x * kappa_inv; asymp_jasb[i] = asymp_jasb[i] + bord_vector[p + 1] * x; } - } + } return QMCKL_SUCCESS; } @@ -1433,7 +1507,7 @@ qmckl_exit_code qmckl_compute_asymp_jasb ( const int64_t bord_num, const double* bord_vector, const double rescale_factor_kappa_ee, - double* const asymp_jasb ); + double* const asymp_jasb ); #+end_src @@ -1457,7 +1531,7 @@ print("asym_one : ", asym_one) print("asymp_jasb[0] : ", asymp_jasb[0]) print("asymp_jasb[1] : ", asymp_jasb[1]) #+end_src - + #+RESULTS: asymp_jasb : asym_one : 0.43340325572525706 : asymp_jasb[0] : 0.5323750557252571 @@ -1500,8 +1574,7 @@ assert(rc == QMCKL_SUCCESS); assert(qmckl_jastrow_provided(context)); double asymp_jasb[2]; -int64_t size_max=0; -rc = qmckl_get_jastrow_asymp_jasb(context, asymp_jasb,&size_max); +rc = qmckl_get_jastrow_asymp_jasb(context, asymp_jasb,2); // calculate asymp_jasb assert(fabs(asymp_jasb[0]-0.5323750557252571) < 1.e-12); @@ -1521,11 +1594,17 @@ f_{ee} = \sum_{i,jelectron.walk_num; + if (size_max < sze) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_get_jastrow_factor_ee", + "Array too small. Expected walk_num"); + } memcpy(factor_ee, ctx->jastrow.factor_ee, sze*sizeof(double)); - (*size_max) = sze; return QMCKL_SUCCESS; } @@ -1828,8 +1912,7 @@ print("factor_ee :",factor_ee) assert(qmckl_jastrow_provided(context)); double factor_ee[walk_num]; -size_max=0; -rc = qmckl_get_jastrow_factor_ee(context, factor_ee, &size_max); +rc = qmckl_get_jastrow_factor_ee(context, factor_ee, walk_num); // calculate factor_ee assert(fabs(factor_ee[0]+4.282760865958113) < 1.e-12); @@ -1848,11 +1931,17 @@ assert(fabs(factor_ee[0]+4.282760865958113) < 1.e-12); *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes -qmckl_exit_code qmckl_get_jastrow_factor_ee_deriv_e(qmckl_context context, double* const factor_ee_deriv_e, int64_t* size_max); +qmckl_exit_code +qmckl_get_jastrow_factor_ee_deriv_e(qmckl_context context, + double* const factor_ee_deriv_e, + const int64_t size_max); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -qmckl_exit_code qmckl_get_jastrow_factor_ee_deriv_e(qmckl_context context, double* const factor_ee_deriv_e, int64_t* size_max) +qmckl_exit_code +qmckl_get_jastrow_factor_ee_deriv_e(qmckl_context context, + double* const factor_ee_deriv_e, + const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; @@ -1867,8 +1956,14 @@ qmckl_exit_code qmckl_get_jastrow_factor_ee_deriv_e(qmckl_context context, doubl assert (ctx != NULL); int64_t sze = ctx->electron.walk_num * 4 * ctx->electron.num; + if (size_max < sze) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_get_jastrow_factor_ee_deriv_e", + "Array too small. Expected 4*walk_num*elec_num"); + } + memcpy(factor_ee_deriv_e, ctx->jastrow.factor_ee_deriv_e, sze * sizeof(double)); - (*size_max) = sze; return QMCKL_SUCCESS; } @@ -2026,7 +2121,7 @@ integer function qmckl_compute_factor_ee_deriv_e_f(context, walk_num, elec_num, dx(1) = ee_distance_rescaled_deriv_e(1, i, j, nw) dx(2) = ee_distance_rescaled_deriv_e(2, i, j, nw) dx(3) = ee_distance_rescaled_deriv_e(3, i, j, nw) - dx(4) = ee_distance_rescaled_deriv_e(4, i, j, nw) + dx(4) = ee_distance_rescaled_deriv_e(4, i, j, nw) if((i .LE. up_num .AND. j .LE. up_num ) .OR. & (i .GT. up_num .AND. j .GT. up_num)) then @@ -2241,8 +2336,7 @@ assert(qmckl_jastrow_provided(context)); // calculate factor_ee_deriv_e double factor_ee_deriv_e[walk_num][4][elec_num]; -size_max=0; -rc = qmckl_get_jastrow_factor_ee_deriv_e(context, &(factor_ee_deriv_e[0][0][0]),&size_max); +rc = qmckl_get_jastrow_factor_ee_deriv_e(context, &(factor_ee_deriv_e[0][0][0]),walk_num*4*elec_num); // check factor_ee_deriv_e assert(fabs(factor_ee_deriv_e[0][0][0]-0.16364894652107934) < 1.e-12); @@ -2264,11 +2358,17 @@ f_{en} = \sum_{i,jelectron.walk_num; + if (size_max < sze) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_get_jastrow_factor_en", + "Array too small. Expected walk_num"); + } memcpy(factor_en, ctx->jastrow.factor_en, sze*sizeof(double)); - (*size_max)=sze; return QMCKL_SUCCESS; } @@ -2562,8 +2667,7 @@ print("factor_en :",factor_en) assert(qmckl_jastrow_provided(context)); double factor_en[walk_num]; -size_max=0; -rc = qmckl_get_jastrow_factor_en(context, factor_en,&size_max); +rc = qmckl_get_jastrow_factor_en(context, factor_en,walk_num); // calculate factor_en assert(fabs(factor_en[0]+5.865822569188727) < 1.e-12); @@ -2579,11 +2683,17 @@ assert(fabs(factor_en[0]+5.865822569188727) < 1.e-12); *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes -qmckl_exit_code qmckl_get_jastrow_factor_en_deriv_e(qmckl_context context, double* const factor_en_deriv_e, int64_t* size_max); +qmckl_exit_code +qmckl_get_jastrow_factor_en_deriv_e(qmckl_context context, + double* const factor_en_deriv_e, + const int64_t size_max); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -qmckl_exit_code qmckl_get_jastrow_factor_en_deriv_e(qmckl_context context, double* const factor_en_deriv_e, int64_t* size_max) +qmckl_exit_code +qmckl_get_jastrow_factor_en_deriv_e(qmckl_context context, + double* const factor_en_deriv_e, + const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; @@ -2598,8 +2708,13 @@ qmckl_exit_code qmckl_get_jastrow_factor_en_deriv_e(qmckl_context context, doubl assert (ctx != NULL); int64_t sze = ctx->electron.walk_num * 4 * ctx->electron.num; + if (size_max < sze) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_get_jastrow_factor_en_deriv_e", + "Array too small. Expected 4*walk_num*elec_num"); + } memcpy(factor_en_deriv_e, ctx->jastrow.factor_en_deriv_e, sze*sizeof(double)); - (*size_max) = sze; return QMCKL_SUCCESS; } @@ -2966,8 +3081,7 @@ assert(qmckl_jastrow_provided(context)); // calculate factor_en_deriv_e double factor_en_deriv_e[walk_num][4][elec_num]; -size_max=0; -rc = qmckl_get_jastrow_factor_en_deriv_e(context, &(factor_en_deriv_e[0][0][0]),&size_max); +rc = qmckl_get_jastrow_factor_en_deriv_e(context, &(factor_en_deriv_e[0][0][0]),walk_num*4*elec_num); // check factor_en_deriv_e assert(fabs(factor_en_deriv_e[0][0][0]-0.11609919541763383) < 1.e-12); @@ -2991,11 +3105,17 @@ assert(fabs(factor_en_deriv_e[0][3][0]+0.9667363412285741 ) < 1.e-12); *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes -qmckl_exit_code qmckl_get_jastrow_een_rescaled_e(qmckl_context context, double* const distance_rescaled, int64_t* size_max); +qmckl_exit_code +qmckl_get_jastrow_een_rescaled_e(qmckl_context context, + double* const distance_rescaled, + const int64_t size_max); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -qmckl_exit_code qmckl_get_jastrow_een_rescaled_e(qmckl_context context, double* const distance_rescaled, int64_t* size_max) +qmckl_exit_code +qmckl_get_jastrow_een_rescaled_e(qmckl_context context, + double* const distance_rescaled, + const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; @@ -3009,9 +3129,14 @@ qmckl_exit_code qmckl_get_jastrow_een_rescaled_e(qmckl_context context, double* qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - size_t sze = ctx->electron.num * ctx->electron.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1); + int64_t sze = ctx->electron.num * ctx->electron.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1); + if (size_max < sze) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_get_jastrow_factor_een_rescaled_e", + "Array too small. Expected ctx->electron.num * ctx->electron.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1)"); + } memcpy(distance_rescaled, ctx->jastrow.een_rescaled_e, sze * sizeof(double)); - (*size_max) = sze; return QMCKL_SUCCESS; } @@ -3223,7 +3348,7 @@ end function qmckl_compute_een_rescaled_e_f #+end_src *** Test - + #+begin_src python :results output :exports none :noweb yes import numpy as np @@ -3287,8 +3412,7 @@ assert(qmckl_electron_provided(context)); double een_rescaled_e[walk_num][(cord_num + 1)][elec_num][elec_num]; -size_max=0; -rc = qmckl_get_jastrow_een_rescaled_e(context, &(een_rescaled_e[0][0][0][0]),&size_max); +rc = qmckl_get_jastrow_een_rescaled_e(context, &(een_rescaled_e[0][0][0][0]),elec_num*elec_num*(cord_num+1)*walk_num); // value of (0,2,1) assert(fabs(een_rescaled_e[0][1][0][2]-0.08084493981483197) < 1.e-12); @@ -3313,11 +3437,17 @@ assert(fabs(een_rescaled_e[0][2][1][5]-0.3424402276009091) < 1.e-12); *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes -qmckl_exit_code qmckl_get_jastrow_een_rescaled_e_deriv_e(qmckl_context context, double* const distance_rescaled, int64_t* size_max); +qmckl_exit_code +qmckl_get_jastrow_een_rescaled_e_deriv_e(qmckl_context context, + double* const distance_rescaled, + const int64_t size_max); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -qmckl_exit_code qmckl_get_jastrow_een_rescaled_e_deriv_e(qmckl_context context, double* const distance_rescaled, int64_t* size_max) +qmckl_exit_code +qmckl_get_jastrow_een_rescaled_e_deriv_e(qmckl_context context, + double* const distance_rescaled, + const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; @@ -3331,9 +3461,14 @@ qmckl_exit_code qmckl_get_jastrow_een_rescaled_e_deriv_e(qmckl_context context, qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - size_t sze = ctx->electron.num * 4 * ctx->electron.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1); + int64_t sze = ctx->electron.num * 4 * ctx->electron.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1); + if (size_max < sze) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_get_jastrow_factor_een_deriv_e", + "Array too small. Expected ctx->electron.num * 4 * ctx->electron.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1)"); + } memcpy(distance_rescaled, ctx->jastrow.een_rescaled_e_deriv_e, sze * sizeof(double)); - (*size_max) = sze; return QMCKL_SUCCESS; } @@ -3654,8 +3789,9 @@ for l in range(0,cord_num+1): #+begin_src c :tangle (eval c_test) //assert(qmckl_electron_provided(context)); double een_rescaled_e_deriv_e[walk_num][(cord_num + 1)][elec_num][4][elec_num]; -size_max=0; -rc = qmckl_get_jastrow_een_rescaled_e_deriv_e(context, &(een_rescaled_e_deriv_e[0][0][0][0][0]),&size_max); +size_max=walk_num*(cord_num + 1)*elec_num*4*elec_num; +rc = qmckl_get_jastrow_een_rescaled_e_deriv_e(context, + &(een_rescaled_e_deriv_e[0][0][0][0][0]),size_max); // value of (0,0,0,2,1) assert(fabs(een_rescaled_e_deriv_e[0][1][0][0][2] + 0.05991352796887283 ) < 1.e-12); @@ -3680,11 +3816,17 @@ assert(fabs(een_rescaled_e_deriv_e[0][2][1][0][5] + 0.5880599146214673 ) < 1. *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes -qmckl_exit_code qmckl_get_jastrow_een_rescaled_n(qmckl_context context, double* const distance_rescaled, int64_t* size_max); +qmckl_exit_code +qmckl_get_jastrow_een_rescaled_n(qmckl_context context, + double* const distance_rescaled, + const int64_t size_max); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -qmckl_exit_code qmckl_get_jastrow_een_rescaled_n(qmckl_context context, double* const distance_rescaled, int64_t* size_max) +qmckl_exit_code +qmckl_get_jastrow_een_rescaled_n(qmckl_context context, + double* const distance_rescaled, + const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; @@ -3698,9 +3840,14 @@ qmckl_exit_code qmckl_get_jastrow_een_rescaled_n(qmckl_context context, double* qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - size_t sze = ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1); + int64_t sze = ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1); + if (size_max < sze) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_get_jastrow_factor_een_deriv_e", + "Array too small. Expected ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1)"); + } memcpy(distance_rescaled, ctx->jastrow.een_rescaled_n, sze * sizeof(double)); - (*size_max)=sze; return QMCKL_SUCCESS; } @@ -3982,8 +4129,8 @@ print(" een_rescaled_n[1, 5, 2] = ",een_rescaled_n[1, 5, 2]) assert(qmckl_electron_provided(context)); double een_rescaled_n[walk_num][(cord_num + 1)][nucl_num][elec_num]; -size_max=0; -rc = qmckl_get_jastrow_een_rescaled_n(context, &(een_rescaled_n[0][0][0][0]),&size_max); +size_max=walk_num*(cord_num + 1)*nucl_num*elec_num; +rc = qmckl_get_jastrow_een_rescaled_n(context, &(een_rescaled_n[0][0][0][0]),size_max); // value of (0,2,1) assert(fabs(een_rescaled_n[0][1][0][2]-0.10612983920006765) < 1.e-12); @@ -4003,11 +4150,17 @@ assert(fabs(een_rescaled_n[0][2][1][5]-0.01343938025140174) < 1.e-12); *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes -qmckl_exit_code qmckl_get_jastrow_een_rescaled_n_deriv_e(qmckl_context context, double* const distance_rescaled, int64_t* size_max); +qmckl_exit_code +qmckl_get_jastrow_een_rescaled_n_deriv_e(qmckl_context context, + double* const distance_rescaled, + const int64_t size_max); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -qmckl_exit_code qmckl_get_jastrow_een_rescaled_n_deriv_e(qmckl_context context, double* const distance_rescaled, int64_t* size_max) +qmckl_exit_code +qmckl_get_jastrow_een_rescaled_n_deriv_e(qmckl_context context, + double* const distance_rescaled, + const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; @@ -4021,9 +4174,14 @@ qmckl_exit_code qmckl_get_jastrow_een_rescaled_n_deriv_e(qmckl_context context, qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - size_t sze = ctx->electron.num * 4 * ctx->nucleus.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1); + int64_t sze = ctx->electron.num * 4 * ctx->nucleus.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1); + if (size_max < sze) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_get_jastrow_factor_een_deriv_e", + "Array too small. Expected ctx->electron.num * 4 * ctx->nucleus.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1)"); + } memcpy(distance_rescaled, ctx->jastrow.een_rescaled_n_deriv_e, sze * sizeof(double)); - (*size_max)=sze; return QMCKL_SUCCESS; } @@ -4282,9 +4440,9 @@ end function qmckl_compute_factor_een_rescaled_n_deriv_e_f end function qmckl_compute_factor_een_rescaled_n_deriv_e #+end_src - + *** Test - + #+begin_src python :results output :exports none :noweb yes import numpy as np @@ -4354,8 +4512,8 @@ print(" een_rescaled_n_deriv_e[2, 1, 6, 2] = ",een_rescaled_n_deriv_e[5, 0, 1, 2 assert(qmckl_electron_provided(context)); double een_rescaled_n_deriv_e[walk_num][(cord_num + 1)][nucl_num][4][elec_num]; -size_max=0; -rc = qmckl_get_jastrow_een_rescaled_n_deriv_e(context, &(een_rescaled_n_deriv_e[0][0][0][0][0]),&size_max); +size_max=walk_num*(cord_num + 1)*nucl_num*4*elec_num; +rc = qmckl_get_jastrow_een_rescaled_n_deriv_e(context, &(een_rescaled_n_deriv_e[0][0][0][0][0]),size_max); // value of (0,2,1) assert(fabs(een_rescaled_n_deriv_e[0][1][0][0][2]+0.07633444246999128 ) < 1.e-12); @@ -5397,11 +5555,17 @@ assert(fabs(dtmp_c[0][1][0][0][0][0] - 0.237440520852232) < 1e-12); *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes -qmckl_exit_code qmckl_get_jastrow_factor_een(qmckl_context context, double* const factor_een, int64_t* size_max); +qmckl_exit_code +qmckl_get_jastrow_factor_een(qmckl_context context, + double* const factor_een, + const int64_t size_max); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -qmckl_exit_code qmckl_get_jastrow_factor_een(qmckl_context context, double* const factor_een, int64_t* size_max) +qmckl_exit_code +qmckl_get_jastrow_factor_een(qmckl_context context, + double* const factor_een, + const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; @@ -5415,9 +5579,14 @@ qmckl_exit_code qmckl_get_jastrow_factor_een(qmckl_context context, double* cons qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - int64_t sze = ctx->electron.walk_num * ctx->electron.num; + int64_t sze = ctx->electron.walk_num; + if (size_max < sze) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_get_jastrow_factor_een", + "Array too small. Expected walk_num"); + } memcpy(factor_een, ctx->jastrow.factor_een, sze*sizeof(double)); - (*size_max)=sze; return QMCKL_SUCCESS; } @@ -5871,8 +6040,7 @@ print("factor_een:",factor_een) assert(qmckl_jastrow_provided(context)); double factor_een[walk_num]; -size_max=0; -rc = qmckl_get_jastrow_factor_een(context, &(factor_een[0]),&size_max); +rc = qmckl_get_jastrow_factor_een(context, &(factor_een[0]),walk_num); assert(fabs(factor_een[0] + 0.37407972141304213) < 1e-12); #+end_src @@ -5886,11 +6054,17 @@ assert(fabs(factor_een[0] + 0.37407972141304213) < 1e-12); *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes -qmckl_exit_code qmckl_get_jastrow_factor_een_deriv_e(qmckl_context context, double* const factor_een_deriv_e, int64_t* size_max); +qmckl_exit_code +qmckl_get_jastrow_factor_een_deriv_e(qmckl_context context, + double* const factor_een_deriv_e, + const int64_t size_max); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -qmckl_exit_code qmckl_get_jastrow_factor_een_deriv_e(qmckl_context context, double* const factor_een_deriv_e, int64_t* size_max) +qmckl_exit_code +qmckl_get_jastrow_factor_een_deriv_e(qmckl_context context, + double* const factor_een_deriv_e, + const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; @@ -5904,9 +6078,14 @@ qmckl_exit_code qmckl_get_jastrow_factor_een_deriv_e(qmckl_context context, doub qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - int64_t sze = ctx->electron.walk_num * ctx->electron.num; + int64_t sze = ctx->electron.walk_num * 4 * ctx->electron.num; + if (size_max < sze) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_get_jastrow_factor_een_deriv_e", + "Array too small. Expected 4*walk_num*elec_num"); + } memcpy(factor_een_deriv_e, ctx->jastrow.factor_een_deriv_e, sze*sizeof(double)); - (*size_max)=sze; return QMCKL_SUCCESS; } @@ -6394,7 +6573,7 @@ end function qmckl_compute_factor_een_deriv_e_f end function qmckl_compute_factor_een_deriv_e #+end_src - + *** Test #+begin_src python :results output :exports none :noweb yes import numpy as np @@ -6451,11 +6630,10 @@ print("factor_een:",factor_een) /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_provided(context)); -double factor_een_deriv_e[walk_num][elec_num]; -size_max=0; -rc = qmckl_get_jastrow_factor_een_deriv_e(context, &(factor_een_deriv_e[0][0]),&size_max); +double factor_een_deriv_e[4][walk_num][elec_num]; +rc = qmckl_get_jastrow_factor_een_deriv_e(context, &(factor_een_deriv_e[0][0][0]),4*walk_num*elec_num); -assert(fabs(factor_een_deriv_e[0][0] + 0.0005481671107226865) < 1e-12); +assert(fabs(factor_een_deriv_e[0][0][0] + 0.0005481671107226865) < 1e-12); #+end_src * End of files :noexport: diff --git a/org/qmckl_mo.org b/org/qmckl_mo.org index 709fb31..d4efb41 100644 --- a/org/qmckl_mo.org +++ b/org/qmckl_mo.org @@ -104,7 +104,7 @@ typedef struct qmckl_mo_basis_struct { double * coefficient; double * mo_vgl; - int64_t mo_vgl_date; + uint64_t mo_vgl_date; int32_t uninitialized; bool provided; diff --git a/org/qmckl_nucleus.org b/org/qmckl_nucleus.org index a562618..0191afe 100644 --- a/org/qmckl_nucleus.org +++ b/org/qmckl_nucleus.org @@ -333,7 +333,7 @@ qmckl_get_nucleus_coord (const qmckl_context context, rc = qmckl_transpose(context, ctx->nucleus.coord, At); if (rc != QMCKL_SUCCESS) return rc; rc = qmckl_double_of_matrix(context, At, coord, size_max); - qmckl_matrix_free(context, At); + qmckl_matrix_free(context, &At); } else { rc = qmckl_double_of_matrix(context, ctx->nucleus.coord, coord, size_max); } @@ -470,6 +470,12 @@ qmckl_set_nucleus_charge(qmckl_context context, ctx->nucleus.charge = qmckl_vector_alloc(context, num); rc = qmckl_vector_of_double(context, charge, num, &(ctx->nucleus.charge)); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, + QMCKL_FAILURE, + "qmckl_set_nucleus_charge", + "Error in vector->double* conversion"); + } <> } @@ -518,7 +524,7 @@ qmckl_set_nucleus_coord(qmckl_context context, const int64_t nucl_num = (int64_t) ctx->nucleus.num; if (ctx->nucleus.coord.data != NULL) { - rc = qmckl_matrix_free(context, ctx->nucleus.coord); + rc = qmckl_matrix_free(context, &(ctx->nucleus.coord)); if (rc != QMCKL_SUCCESS) return rc; } @@ -662,14 +668,14 @@ assert(!qmckl_nucleus_provided(context)); rc = qmckl_get_nucleus_coord (context, 'N', nucl_coord2, 3*nucl_num); assert(rc == QMCKL_SUCCESS); for (size_t k=0 ; k<3 ; ++k) { - for (size_t i=0 ; ipoint.coord, coord, size_max); } @@ -301,7 +301,7 @@ qmckl_set_point (qmckl_context context, if (ctx->point.num < num) { if (ctx->point.coord.data != NULL) { - rc = qmckl_matrix_free(context, ctx->point.coord); + rc = qmckl_matrix_free(context, &(ctx->point.coord)); assert (rc == QMCKL_SUCCESS); } @@ -318,8 +318,17 @@ qmckl_set_point (qmckl_context context, ctx->point.num = num; if (transp == 'T') { - memcpy(ctx->point.coord.data, coord, 3*num*sizeof(double)); + double *a = ctx->point.coord.data; +#ifdef HAVE_OPENMP + #pragma omp for +#endif + for (int64_t i=0 ; i<3*num ; ++i) { + a[i] = coord[i]; + } } else { +#ifdef HAVE_OPENMP + #pragma omp for +#endif for (int64_t i=0 ; ipoint.coord, i, 0) = coord[3*i ]; qmckl_mat(ctx->point.coord, i, 1) = coord[3*i+1]; @@ -328,8 +337,8 @@ qmckl_set_point (qmckl_context context, } /* Increment the date of the context */ - ctx->date += 1UL; - ctx->point.date = ctx->date; + rc = qmckl_context_touch(context); + assert (rc == QMCKL_SUCCESS); return QMCKL_SUCCESS; diff --git a/org/qmckl_trexio.org b/org/qmckl_trexio.org index e08f044..12b7153 100644 --- a/org/qmckl_trexio.org +++ b/org/qmckl_trexio.org @@ -429,6 +429,11 @@ qmckl_trexio_read_ao_X(qmckl_context context, trexio_t* const file) /* Reformat data */ rc = qmckl_set_ao_basis_nucleus_index(context, nucleus_index, nucleus_num); + if (rc != QMCKL_SUCCESS) { + qmckl_free(context, nucleus_index); + nucleus_index = NULL; + return rc; + } for (int i=shell_num-1 ; i>=0 ; --i) { const int k = tmp_array[i]; @@ -1048,21 +1053,44 @@ qmckl_trexio_read_mo_X(qmckl_context context, trexio_t* const file) * Read everything #+begin_src c :tangle (eval h_func) -qmckl_exit_code qmckl_trexio_read(const qmckl_context context, const char* file_name); +qmckl_exit_code +qmckl_trexio_read(const qmckl_context context, + const char* file_name, + const int64_t size_max); + #+end_src + + #+begin_src f90 :tangle (eval fh_func) :comments org :exports none + interface + integer(c_int32_t) function qmckl_trexio_read & + (context, file_name, size_max) & + bind(C) + use, intrinsic :: iso_c_binding + import + implicit none + + integer (c_int64_t) , intent(in) , value :: context + integer (c_int64_t) , intent(in) , value :: size_max + character(c_char ) , intent(in) :: file_name(size_max) + + end function qmckl_trexio_read + end interface #+end_src #+begin_src c :tangle (eval c) qmckl_exit_code -qmckl_trexio_read(const qmckl_context context, const char* file_name) +qmckl_trexio_read(const qmckl_context context, const char* file_name, const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return false; } qmckl_exit_code rc; + char file_name_new[size_max+1]; + strncpy(file_name_new, file_name, size_max+1); + file_name_new[size_max] = '\0'; #ifdef HAVE_TREXIO - trexio_t* file = qmckl_trexio_open_X(file_name, &rc); + trexio_t* file = qmckl_trexio_open_X(file_name_new, &rc); if (file == NULL) { trexio_close(file); return qmckl_failwith( context, @@ -1141,7 +1169,7 @@ strncat(fname, "/chbrclf", 255); printf("Test file: %s\n", fname); rc = qmckl_set_electron_walk_num(context, walk_num); -rc = qmckl_trexio_read(context, fname); +rc = qmckl_trexio_read(context, fname, 255); if (rc != QMCKL_SUCCESS) { printf("%s\n", qmckl_string_of_error(rc)); diff --git a/tools/build_doc.sh b/tools/build_doc.sh index 519b417..de7f21b 100755 --- a/tools/build_doc.sh +++ b/tools/build_doc.sh @@ -9,12 +9,7 @@ if [[ -z ${srcdir} ]] ; then exit 1 fi -if [[ -z ${top_builddir} ]] ; then - echo "Error: srcdir environment variable is not defined" - exit 1 -fi - -readonly DOCS=${top_builddir}/share/doc/qmckl/ +readonly DOCS=share/doc/qmckl/ readonly ORG=${srcdir}/org/ readonly HTMLIZE=${DOCS}/html/htmlize.el readonly CONFIG_DOC=${srcdir}/tools/config_doc.el @@ -47,7 +42,7 @@ function extract_doc() for i in $@ do exported=${i%.org}.exported - exported=${top_builddir}/src/$(basename $exported) + exported=${srcdir}/src/$(basename $exported) NOW=$(date +"%m%d%H%M.%S") extract_doc ${i} &> $exported rc=$? diff --git a/tools/build_makefile.py b/tools/build_makefile.py index f85bcdb..cb4b33e 100755 --- a/tools/build_makefile.py +++ b/tools/build_makefile.py @@ -49,8 +49,8 @@ def main(): f_test_o = "tests/test_"+i+"_f.$(OBJEXT)" c_test = "tests/test_"+i+".c" f_test = "tests/test_"+i+"_f.F90" - html = "share/doc/qmckl/html/"+i+".html" - text = "share/doc/qmckl/text/"+i+".txt" + html = "$(htmldir_loc)/"+i+".html" + text = "$(textdir_loc)/"+i+".txt" i="src/"+i @@ -267,7 +267,7 @@ def main(): for f in sorted(DEPS_DOC.keys()): output += [ DEPS_DOC[f] + ": $(srcdir)/" + f + " $(htmlize_el)", - "\t$(export_verbose)top_builddir=$(abs_top_builddir) srcdir=$(abs_srcdir) $(srcdir)/tools/missing bash $(srcdir)/tools/build_doc.sh $(srcdir)/"+f, + "\tsrcdir=$(abs_srcdir) $(srcdir)/tools/missing bash $(srcdir)/tools/build_doc.sh $(srcdir)/"+f, "" ] output += ["endif"]