TREX/qmckl
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-06-30 00:44:52 +02:00
Code Issues Projects Releases Wiki Activity

Merge branch 'master' of github.com:TREX-CoE/qmckl

Browse Source
This commit is contained in:
Anthony Scemama 2022-01-31 16:52:01 +01:00
commit 74e5da5c21
44 changed files with 131149 additions and 128433 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
.github/workflows/test-build.yml vendored
View File

@ -20,8 +20,9 @@ jobs:
- name: Install trexio
run: |
export VERSION=1.1.0
wget https://github.com/TREX-CoE/trexio/releases/download/v${VERSION}/trexio-${VERSION}.tar.gz
export TAG=v2.1.0
export VERSION=2.1.0
wget https://github.com/TREX-CoE/trexio/releases/download/${TAG}/trexio-${VERSION}.tar.gz
tar -zxf trexio-${VERSION}.tar.gz
cd trexio-${VERSION}
./configure --prefix=/usr
@ -31,7 +32,7 @@ jobs:
- name: Build
run: |
./autogen.sh
QMCKL_DEVEL=1 ./configure --enable-silent-rules --enable-maintainer-mode --enable-debug
./configure --enable-silent-rules --enable-debug
make -j 8
- name: Run test
@ -80,7 +81,7 @@ jobs:
# - name: Run test
# run: |
# ./autogen.sh
# QMCKL_DEVEL=1 ./configure --enable-silent-rules --enable-maintainer-mode --enable-debug
# ./configure --enable-silent-rules --enable-debug
# make -j 8
# make -j check
# make distcheck

11
.github/workflows/vfc_test_workflow.yml vendored
View File

@ -25,17 +25,18 @@ jobs:
run: |
ln -s /usr/bin/python3 /usr/bin/python
apt update
apt -y install emacs pkg-config wget libhdf5-dev
apt -y install emacs pkg-config wget libhdf5-dev libblas-dev liblapack-dev
- name: Install trexio
run: |
export VERSION=1.1.0
export VERSION=2.1.0
wget https://github.com/TREX-CoE/trexio/releases/download/v${VERSION}/trexio-${VERSION}.tar.gz
tar -zxf trexio-${VERSION}.tar.gz
cd trexio-${VERSION}
./configure --prefix=/usr
make -j 8
sudo make install
./configure --prefix=/usr CC="gcc-7" FC="gfortran-7"
# modify LDFLAGS to include -lhdf5_hl because autoconf sometime fails to detect the HL component
make LDFLAGS="-L/usr/lib/x86_64-linux-gnu/hdf5/serial -lhdf5_hl"
make install
- name: Run tests
run: vfc_ci test -g -r

58
Makefile.am
View File

@ -47,27 +47,27 @@ pkgconfig_DATA = pkgconfig/qmckl.pc
qmckl_h = include/qmckl.h
include_HEADERS = $(qmckl_h)
qmckl_f = share/qmckl/fortran/qmckl_f.f90
test_qmckl_f = tests/qmckl_f.f90
test_qmckl_fo = tests/qmckl_f.o
src_qmckl_f = src/qmckl_f.f90
src_qmckl_fo = src/qmckl_f.o
header_tests = tests/chbrclf.h tests/n2.h
fortrandir = $(datadir)/fortran
fortran_DATA = $(qmckl_f)
fortrandir = $(datadir)/qmckl/fortran
fortran_DATA = $(src_qmckl_f)
QMCKL_TEST_DIR = $(abs_srcdir)/share/qmckl/test_data/
AM_CPPFLAGS = -I$(srcdir)/src -I$(srcdir)/include
AM_CPPFLAGS += -DQMCKL_TEST_DIR="\"$(abs_srcdir)/share/qmckl/test_data/\""
AM_CPPFLAGS = -I$(top_builddir)/src -I$(top_builddir)/include
AM_CPPFLAGS += -I$(srcdir)/src -I$(srcdir)/include -I$(QMCKL_TEST_DIR)
AM_CPPFLAGS += -DQMCKL_TEST_DIR=\"$(QMCKL_TEST_DIR)\"
lib_LTLIBRARIES = src/libqmckl.la
src_libqmckl_la_SOURCES = $(qmckl_h) $(src_qmckl_f) $(C_FILES) $(F_FILES) $(H_PRIVATE_FUNC_FILES) $(H_PRIVATE_TYPE_FILES) $(header_tests)
src_libqmckl_la_SOURCES = $(qmckl_h) $(src_qmckl_f) $(C_FILES) $(F_FILES) $(H_PRIVATE_FUNC_FILES) $(H_PRIVATE_TYPE_FILES)
CLEANFILES+=$(test_qmckl_f) $(test_qmckl_fo) $(src_qmckl_f) $(src_qmckl_fo) $(test_qmckl_o) $(src_qmckl_o)
CLEANFILES+=$(test_qmckl_fo) $(src_qmckl_fo) $(test_qmckl_o) $(src_qmckl_o) $(FH_TYPE_FILES) $(FH_FUNC_FILES)
htmlize_el=share/doc/qmckl/html/htmlize.el
include generated.mk
@ -77,20 +77,18 @@ ln_s_verbose_ = $(ln_s_verbose_@AM_DEFAULT_V@)
ln_s_verbose_0 = @echo " LN_S $<";
htmldir = $(docdir)/html
dist_html_DATA = $(HTML_FILES) \
share/doc/qmckl/html/index.html \
share/doc/qmckl/html/qmckl.css
textdir = $(docdir)/text
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
dist_text_DATA = $(TEXT_FILES)
share/doc/qmckl/html/index.html: share/doc/qmckl/html/README.html
$(ln_s_verbose)cd share/doc/qmckl/html/ && \
rm -rf index.html && \
$(LN_S) README.html index.html
html: $(dist_html_DATA)
html-local: $(dist_html_DATA)
text: $(dist_text_DATA)
doc: html text
@ -98,16 +96,15 @@ doc: html text
if QMCKL_DEVEL
if VFC_CI
AM_LDFLAGS=-lvfc_probes -lvfc_probes_f
endif
dist_src_DATA = $(ORG_FILES) $(TANGLED_FILES) $(EXPORTED_FILES) $(qmckl_f) $(test_qmckl_f)
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) $(qmckl_f) $(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)
CLEANFILES += $(BUILT_SOURCES) $(C_TEST_FILES) $(F_TEST_FILES) $(TANGLED_FILES) $(C_TEST_FILES) $(F_TEST_FILES) $(qmckl_f) $(src_qmckl_f) $(test_qmckl_f) $(qmckl_h) $(HTML_FILES) $(TEXT_FILES) share/doc/qmckl/html/index.html $(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)
EXTRA_DIST += \
tools/build_doc.sh \
@ -147,25 +144,22 @@ SUFFIXES = .f90 .h .org .c _f.f90 _func.h _type.h _private_func.h _private_type.
$(test_qmckl_f): $(src_qmckl_f)
cp $(src_qmckl_f) $(test_qmckl_f)
$(qmckl_f): $(src_qmckl_f)
cp $(src_qmckl_f) $(qmckl_f)
$(qmckl_h): $(H_FUNC_FILES) $(H_TYPE_FILES)
$(cat_h_verbose)top_builddir=$(top_builddir) srcdir=$(srcdir) qmckl_h=$(qmckl_h) $(srcdir)/tools/build_qmckl_h.sh
$(cat_h_verbose)top_builddir=$(abs_top_builddir) srcdir=$(abs_srcdir) qmckl_h=$(qmckl_h) $(srcdir)/tools/build_qmckl_h.sh
$(src_qmckl_f): $(FH_FUNC_FILES) $(FH_TYPE_FILES)
$(cat_h_verbose)top_builddir=$(top_builddir) srcdir=$(srcdir) src_qmckl_f=$(src_qmckl_f) $(srcdir)/tools/build_qmckl_f.sh
echo $(FH_TYPE_FILES) $(FH_FUNC_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):
$(srcdir)/tools/install_htmlize.sh $(htmlize_el)
$(MKDIR_P) $(htmldir_loc)
$(MKDIR_P) $(textdir_loc)
abs_srcdir=$(abs_srcdir) $(srcdir)/tools/install_htmlize.sh $(htmlize_el)
tests/chbrclf.h: $(qmckl_h)
tests/n2.h: $(qmckl_h)
generated.mk: $(ORG_FILES)
top_builddir=$(top_builddir) srcdir=$(srcdir) qmckl_h=$(qmckl_h) src_qmckl_f=$(src_qmckl_f) $(PYTHON) $(srcdir)/tools/build_makefile.py
cppcheck: cppcheck.out
cppcheck.out: $(qmckl_h)
@ -175,7 +169,7 @@ cppcheck.out: $(qmckl_h)
--suppress="unusedFunction" \
--suppress="missingIncludeSystem" \
--language=c --std=c99 -rp --platform=unix64 \
-I../include *.c *.h 2>../$@
-I$(srcdir)/include -I$(top_builddir)/include *.c *.h 2>../$@
.PHONY: cppcheck

4
README.md
View File

@ -28,7 +28,7 @@ in org-mode files and produce the source code and the documentation from these f
```
./autogen.sh
QMCKL_DEVEL=1 ./configure --prefix=$PWD/_install --enable-silent-rules --enable-maintainer-mode
./configure --prefix=$PWD/_install
make
make check
@ -54,7 +54,7 @@ by the preprocessor otherwise). To enable vfc_ci support, the library should be
configured with the following command :
```
QMCKL_DEVEL=1 ./configure --prefix=$PWD/_install \ --enable-silent-rules --enable-maintainer-mode --enable-vfc_ci --host=x86_64 \ CC="verificarlo-f" FC="verificarlo-f"
./configure --prefix=$PWD/_install \ --enable-vfc_ci --host=x86_64 \ CC="verificarlo-f" FC="verificarlo-f"
```
where CC and FC are set to verificarlo-f, and support is explicitely enabled

8
autogen.sh
View File

@ -1,7 +1,5 @@
#!/bin/bash
python ./tools/build_makefile.py
autoreconf -i
echo "To configure in maintainer mode, use:
$ QMCKL_DEVEL=1 ./configure --enable-silent-rules --enable-maintainer-mode
"
export srcdir="."
python ${srcdir}/tools/build_makefile.py
autoreconf -i -Wall --no-recursive

116
configure.ac
View File

@ -35,11 +35,18 @@
AC_PREREQ([2.69])
AC_INIT([qmckl],[0.1.1],
[https://github.com/TREX-CoE/qmckl/issues], [],
[https://trex-coe.github.io/qmckl/index.html])
AC_CONFIG_AUX_DIR(tools)
AC_INIT([qmckl],[0.1.1],[https://github.com/TREX-CoE/qmckl/issues],[],[https://trex-coe.github.io/qmckl/index.html])
AC_CONFIG_AUX_DIR([tools])
AM_INIT_AUTOMAKE([subdir-objects color-tests parallel-tests silent-rules 1.11])
# Activate developer mode when the source is the git repository.
# Otherwise, it is the source distribution and the developer mode should not be activated.
AS_IF([test -d ${srcdir}/.git],
[enable_maintainer_mode="yes"],
[enable_maintainer_mode="no"]
)
AM_PROG_AR
AM_MAINTAINER_MODE()
LT_INIT
AC_CONFIG_SRCDIR([configure.ac])
@ -58,7 +65,7 @@ AC_SUBST([VERSION_PATCH])
#AM_SILENT_RULES(yes)
#AC_SUBST(SHARED_VERSION_INFO)
#AM_ENABLE_SHARED(no) dnl shared libs cause too many headaches to be default
#AC_ENABLE_SHARED([no]) dnl shared libs cause too many headaches to be default
AC_LANG(C)
@ -68,11 +75,12 @@ AC_PROG_CC
m4_version_prereq([2.70],[], [AC_PROG_CC_C99])
AS_IF([test "$ac_cv_prog_cc_c99" = "no"], [AC_MSG_ERROR([The compiler does not support C99])])
AC_PROG_CC_C_O
AC_PROG_F77
AC_PROG_FC
AC_PROG_FC_C_O
AC_FC_SRCEXT([f90])
AC_FC_FREEFORM
AC_PROG_LIBTOOL
LT_INIT
AC_PROG_INSTALL
AC_PROG_LN_S
PKG_PROG_PKG_CONFIG([])
@ -90,7 +98,7 @@ 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, [AC_HELP_STRING([--with-openmp],[enable OpenMP])], with_omp=$withval, with_omp=no)
#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([],
@ -98,9 +106,55 @@ AC_CHECK_LIB([pthread], [pthread_create])
# CFLAGS="${CFLAGS} ${OPENMP_CFLAGS}"
#fi
# CHAMELEON
AC_ARG_WITH(chameleon,
AS_HELP_STRING([--without-chameleon],
[Do not use Chameleon. Default: auto-detect]), [
case "$with_chameleon" in
no)
: ;;
yes)
PKG_CHECK_MODULES([LIBCHAMELEON],[chameleon >= 1.0.0],
[PKG_CFLAGS="$PKG_CFLAGS $LIBCHAMELEON_CFLAGS"
PKG_LIBS="$PKG_LIBS $LIBCHAMELEON_LIBS"]
,[
## something went wrong.
## try to find the package without pkg-config
## check that the library is actually new enough.
## by testing for a 1.0.0+ function which we use
AC_CHECK_LIB(chameleon,CHAMELEON_finalize,[LIBCHAMELEON_LIBS="-lchameleon"])
])
AH_TEMPLATE(HAVE_CHAMELEON,[Chameleon support is available])
;;
*)
if test ! -d "$withval" ; then
AC_MSG_ERROR([--with-chameleon path does not point to a directory])
fi
LIBCHAMELEON_LIBS="-L$with_chameleon/lib -lchameleon -lchameleon_starpu -lhqr -lcoreblas"
LIBCHAMELEON_CFLAGS="-I$with_chameleon/include $CFLAGS"
LIBCHAMELEON_CPPFLAGS="-I$with_chameleon/include $CPPFLAGS"
esac
])
if test "x$LIBCHAMELEON_LIBS" != "x" ; then
LIBS="$LIBS $LIBCHAMELEON_LIBS"
CFLAGS="$CFLAGS $LIBCHAMELEON_CFLAGS"
CPPFLAGS="$CPPFLAGS $LIBCHAMELEON_CPPFLAGS"
AC_CHECK_HEADERS([chameleon.h], [], [AC_MSG_ERROR("chamelon.h not found")])
AC_DEFINE_UNQUOTED([HAVE_CHAMELEON],1,[CHAMELEON support is available])
fi
#AS_IF([test "x$with_chameleon" != "xno"], [
#])
AC_MSG_NOTICE([CHAMELEON library support: ${with_CHAMELEON:=auto} ${LIBCHAMELEON_PATH} ${LIBCHAMELEON_LIBS}])
# TREXIO
AC_ARG_WITH(trexio, [AC_HELP_STRING([--without-trexio],[disable support for TREXIO])],
AC_ARG_WITH(trexio, [AS_HELP_STRING([--without-trexio],[disable support for TREXIO])],
with_trexio=$withval, with_trexio=yes)
AS_IF([test "x$with_trexio" != xno], [
@ -117,14 +171,30 @@ PKG_CFLAGS="$PKG_CFLAGS $TREXIO_CFLAGS"
PKG_LIBS="$PKG_LIBS $TREXIO_LIBS"
## BLAS
#AX_BLAS([], [AC_MSG_ERROR([BLAS was not found.])])
AX_BLAS([], [AC_MSG_ERROR([BLAS was not found.])])
## LAPACK
#AX_LAPACK([], [AC_MSG_ERROR([LAPACK was not found.])])
AX_LAPACK([], [AC_MSG_ERROR([LAPACK was not found.])])
# Specific options required with some compilers
case $FC in
ifort*)
FCFLAGS="$FCFLAGS -nofor-main"
;;
gfortran*)
# Order is important here
FCFLAGS="-cpp $FCFLAGS"
;;
esac
# Options.
AC_ARG_ENABLE(debug, [AC_HELP_STRING([--enable-debug],[compile for debugging])], ok=$enableval, ok=no)
AC_ARG_ENABLE(debug, [AS_HELP_STRING([--enable-debug],[compile for debugging])], ok=$enableval, ok=no)
if test "$ok" = "yes"; then
if test "$GCC" = "yes"; then
CFLAGS="$CFLAGS \
@ -146,20 +216,20 @@ if test "$ok" = "yes"; then
ARGS="${ARGS} debug"
fi
AC_ARG_ENABLE(malloc-trace, [AC_HELP_STRING([--enable-malloc-trace],[use debug malloc/free])], ok=$enableval, ok=no)
AC_ARG_ENABLE(malloc-trace, [AS_HELP_STRING([--enable-malloc-trace],[use debug malloc/free])], ok=$enableval, ok=no)
if test "$ok" = "yes"; then
AC_DEFINE(MALLOC_TRACE,"malloc_trace.dat",[Define to use debugging malloc/free])
ARGS="${ARGS} malloc-trace"
fi
AC_ARG_ENABLE(prof, [AC_HELP_STRING([--enable-prof],[compile for profiling])], ok=$enableval, ok=no)
AC_ARG_ENABLE(prof, [AS_HELP_STRING([--enable-prof],[compile for profiling])], ok=$enableval, ok=no)
if test "$ok" = "yes"; then
CFLAGS="${CFLAGS} -pg"
AC_DEFINE(ENABLE_PROF,1,[Define when using the profiler tool])
ARGS="${ARGS} prof"
fi
AC_ARG_WITH(efence, [AC_HELP_STRING([--with-efence],[use ElectricFence library])], ok=$withval, ok=no)
AC_ARG_WITH(efence, [AS_HELP_STRING([--with-efence],[use ElectricFence library])], ok=$withval, ok=no)
if test "$ok" = "yes"; then
AC_CHECK_LIB(efence, malloc)
ARGS="${ARGS} efence"
@ -189,14 +259,18 @@ AC_TYPE_UINT64_T
AC_CHECK_FUNCS([memset strerror])
# Development mode
if test "x$enable_maintainer_mode" == "xyes" ; then
QMCKL_DEVEL=" -- Developer mode"
else
QMCKL_DEVEL=""
fi
AM_CONDITIONAL([QMCKL_DEVEL],[test "x$QMCKL_DEVEL" != x])
if test "x${QMCKL_DEVEL}" != "x"; then
QMCKL_DEVEL=" -- Developer mode"
AC_PROG_AWK
AM_PATH_PYTHON
${PYTHON} ${srcdir}/tools/build_makefile.py
AC_CHECK_PROGS([EMACS],[emacs26 emacs],[no])
if test x${EMACS} == xno ; then
@ -224,13 +298,6 @@ AC_ARG_ENABLE([vfc_ci],
esac],[vfc_ci=false])
AM_CONDITIONAL([VFC_CI], [test x$vfc_ci = xtrue])
# Enable Fortran preprocessor
if test "$FC" = "gfortran"; then
AC_MSG_NOTICE(gfortran detected)
# Arguments order is important here
FCFLAGS="-cpp $FCFLAGS"
fi
if test "$FC" = "verificarlo-f"; then
AC_MSG_NOTICE(verificarlo-f detected)
# Arguments order is important here
@ -241,7 +308,7 @@ fi
#mkl-dynamic-lp64-seq
PKG_LIBS="$PKG_LIBS $LIBS"
LIBS="$LAPACK_LIBS $BLAS_LIBS $PKG_LIBS"
LIBS="$BLAS_LIBS $LAPACK_LIBS $BLAS_LIBS $PKG_LIBS"
CFLAGS="$CFLAGS $PKG_CFLAGS"
AC_SUBST([PKG_LIBS])
AC_SUBST([PKG_CFLAGS])
@ -269,6 +336,7 @@ FC..............: ${FC}
FCLAGS..........: ${FCFLAGS}
LDFLAGS:........: ${LDFLAGS}
LIBS............: ${LIBS}
USE CHAMELEON...: ${with_chameleon}
Package features:
${ARGS}

56
org/README.org
View File

@ -16,23 +16,25 @@ grep TITLE $(cat table_of_contents) | tr ':' ' '
#+end_src
#+RESULTS: toc
| qmckl.org | #+TITLE | Introduction | | |
| qmckl_error.org | #+TITLE | Error | handling | |
| qmckl_context.org | #+TITLE | Context | | |
| qmckl_memory.org | #+TITLE | Memory | management | |
| qmckl_numprec.org | #+TITLE | Numerical | precision | |
| qmckl_distance.org | #+TITLE | Inter-particle | distances | |
| qmckl_nucleus.org | #+TITLE | Nucleus | | |
| qmckl_electron.org | #+TITLE | Electrons | | |
| qmckl_ao.org | #+TITLE | Atomic | Orbitals | |
| qmckl_mo.org | #+TITLE | Molecular | Orbitals | |
| qmckl_jastrow.org | #+TITLE | Jastrow | Factor | |
| qmckl_sherman_morrison_woodbury.org | #+TITLE | Sherman-Morrison-Woodbury | | |
| qmckl_utils.org | #+TITLE | Utility | functions | |
| qmckl_blas.org | #+TITLE | BLAS | functions | |
| qmckl_trexio.org | #+TITLE | TREXIO | I/O | library |
| qmckl_verificarlo.org | #+TITLE | Verificarlo | CI | |
| qmckl_tests.org | #+TITLE | Data | for | Tests |
| qmckl.org | #+TITLE | Introduction | | |
| qmckl_ao.org | #+TITLE | Atomic | Orbitals | |
| qmckl_blas.org | #+TITLE | BLAS | functions | |
| qmckl_context.org | #+TITLE | Context | | |
| qmckl_determinant.org | #+TITLE | Slater | Determinant | |
| qmckl_distance.org | #+TITLE | Inter-particle | distances | |
| qmckl_electron.org | #+TITLE | Electrons | | |
| qmckl_error.org | #+TITLE | Error | handling | |
| qmckl_jastrow.org | #+TITLE | Jastrow | Factor | |
| qmckl_local_energy.org | #+TITLE | Local | Energy | |
| qmckl_memory.org | #+TITLE | Memory | management | |
| qmckl_mo.org | #+TITLE | Molecular | Orbitals | |
| qmckl_numprec.org | #+TITLE | Numerical | precision | |
| qmckl_nucleus.org | #+TITLE | Nucleus | | |
| qmckl_sherman_morrison_woodbury.org | #+TITLE | Sherman-Morrison-Woodbury | | |
| qmckl_utils.org | #+TITLE | Utility | functions | |
| qmckl_trexio.org | #+TITLE | TREXIO | I/O | library |
| qmckl_verificarlo.org | #+TITLE | Verificarlo | CI | |
| qmckl_tests.org | #+TITLE | Data | for | Tests |
#+begin_src python :var data=toc :exports results :results raw
result = []
@ -45,19 +47,21 @@ return '\n'.join(result)
#+RESULTS:
- [[./qmckl.html][Introduction]]
- [[./qmckl_error.html][Error handling]]
- [[./qmckl_context.html][Context]]
- [[./qmckl_memory.html][Memory management]]
- [[./qmckl_numprec.html][Numerical precision]]
- [[./qmckl_distance.html][Inter-particle distances]]
- [[./qmckl_nucleus.html][Nucleus]]
- [[./qmckl_electron.html][Electrons]]
- [[./qmckl_ao.html][Atomic Orbitals]]
- [[./qmckl_mo.html][Molecular Orbitals]]
- [[./qmckl_blas.html][BLAS functions]]
- [[./qmckl_context.html][Context]]
- [[./qmckl_determinant.html][Slater Determinant]]
- [[./qmckl_distance.html][Inter-particle distances]]
- [[./qmckl_electron.html][Electrons]]
- [[./qmckl_error.html][Error handling]]
- [[./qmckl_jastrow.html][Jastrow Factor]]
- [[./qmckl_local_energy.html][Local Energy]]
- [[./qmckl_memory.html][Memory management]]
- [[./qmckl_mo.html][Molecular Orbitals]]
- [[./qmckl_numprec.html][Numerical precision]]
- [[./qmckl_nucleus.html][Nucleus]]
- [[./qmckl_sherman_morrison_woodbury.html][Sherman-Morrison-Woodbury]]
- [[./qmckl_utils.html][Utility functions]]
- [[./qmckl_blas.html][BLAS functions]]
- [[./qmckl_trexio.html][TREXIO I/O library]]
- [[./qmckl_verificarlo.html][Verificarlo CI]]
- [[./qmckl_tests.html][Data for Tests]]

132
org/qmckl.org
View File

@ -3,6 +3,33 @@
#+SETUPFILE: ../tools/theme.setup
# -*- mode: org -*-
* Installing QMCkl
The latest version fo QMCkl can be downloaded
[[https://github.com/TREX-CoE/qmckl/releases/latest][here]], and the source code is accessible on the
[[https://github.com/TREX-CoE/qmckl][GitHub repository]].
** Installing from the released tarball (for end users)
QMCkl is built with GNU Autotools, so the usual
=configure ; make ; make check ; make install= scheme will be used.
As usual, the C compiler can be specified with the ~CC~ variable
and the Fortran compiler with the ~FC~ variable. The compiler
options are defined using ~CFLAGS~ and ~FCFLAGS~.
** Installing from the source repository (for developers)
To compile from the source repository, additional dependencies are
required to generated the source files:
- Emacs >= 26
- Autotools
- Python3
When the repository is downloaded, the Makefile is not yet
generated, as well as the configure script. =./autogen.sh= has
to be executed first.
* Using QMCkl
The =qmckl.h= header file installed in the =${prefix}/include= directory
@ -59,6 +86,9 @@ Both files are located in the =include/= directory.
Moreover, within the Emacs text editor the source code blocks can be executed
interactively, in the same spirit as Jupyter notebooks.
Note that Emacs is not needed for end users because the distributed
tarball contains the generated source code.
** Source code editing
For a tutorial on literate programming with org-mode, follow [[http://www.howardism.org/Technical/Emacs/literate-programming-tutorial.html][this link]].
@ -80,36 +110,50 @@ Both files are located in the =include/= directory.
** Choice of the programming language
Most of the codes of the [[https://trex-coe.eu][TREX CoE]] are written in Fortran with some scripts in
Bash and Python. Outside of the CoE, Fortran is also important (Casino, Amolqc),
and other important languages used by the community are C and C++ (QMCPack,
QWalk), and Julia is gaining in popularity. The library we design should be
compatible with all of these languages. The QMCkl API has to be compatible
with the C language since libraries with a C-compatible API can be used in
every other language.
Most of the codes of the [[https://trex-coe.eu][TREX CoE]] are written in Fortran with some
scripts in Bash and Python. Outside of the CoE, Fortran is also
important in QMC codes (Casino, Amolqc), and other important
languages used by the community are C and C++ (QMCPack, QWalk),
Julia and Rust are gaining in popularity. We want QMCkl to be
compatible with all of these languages, so the QMCkl API has to be
compatible with the C language since libraries with a C-compatible
API can be used in every other language.
High-performance versions of the QMCkl, with the same API, will be rewritten by
the experts in HPC. These optimized libraries will be tuned for specific
architectures, among which we can cite x86 based processors, and GPU
accelerators. Nowadays, the most efficient software tools to take advantage of
low-level features of the processor (intrinsics) and of GPUs are for C++
developers. It is highly probable that the optimized implementations will be
written in C++, and this is agreement with our choice to make the API
C-compatible.
High-performance versions of QMCkl, with the same API, can be
rewritten by HPC experts. These optimized libraries will be tuned
for specific architectures, among which we can cite x86 based
processors, and GPU accelerators. Nowadays, the most efficient
software tools to take advantage of low-level features
(intrinsics, prefetching, aligned or pinned memory allocation,
...) are for C++ developers. It is highly probable that optimized
implementations will be written in C++, but as the API is
C-compatible this doesn't pose any problem for linking the library
in other languages.
Fortran is one of the most common languages used by the community, and is simple
enough to make the algorithms readable both by experts in QMC, and experts in
HPC. Hence we propose in this pedagogical implementation of QMCkl to use Fortran
to express the QMC algorithms. As the main languages of the library is C, this
implies that the exposed C functions call the Fortran routine. However, for
internal functions related to system programming, the C language is more natural
than Fortran.
Fortran is one of the most common languages used by the community,
and is simple enough to make the algorithms readable both by
experts in QMC, and experts in HPC. Hence we propose in this
pedagogical implementation of QMCkl to use Fortran to express the
QMC algorithms. However, for internal functions related to system
programming, the C language is more natural than Fortran.
The Fortran source files should provide a C interface using the
~iso_c_binding~ module. The name of the Fortran source files should end with
=_f.f90= to be properly handled by the =Makefile=. The names of the functions
defined in Fortran should be the same as those exposed in the API suffixed by
=_f=.
As QMCkl appears like a C library, for each Fortran function there
is an ~iso_c_binding~ interface to make the Fortran function
callable from C. It is this C interface which is exposed to the
user. As a consequence, the Fortran users of the library never
call directly the Fortran routines, but call instead the C binding
function and an ~iso_c_binding~ is still required:
#+begin_example
ISO_C_BINDING ISO_C_BINDING
Fortran ---------------> C ---------------> Fortran
#+end_example
The name of the Fortran source files should end with =_f.f90= to
be properly handled by the =Makefile= and to avoid collision of
object files (=*.o=) with the compiled C source files. The names
of the functions defined in Fortran should be the same as those
exposed in the API suffixed by =_f=.
For more guidelines on using Fortran to generate a C interface, see
[[http://fortranwiki.org/fortran/show/Generating+C+Interfaces][this link]].
@ -123,6 +167,8 @@ Both files are located in the =include/= directory.
#+begin_src bash
cppcheck --addon=cert --enable=all *.c &> cppcheck.out
# or
make cppcheck ; cat cppcheck.out
#+end_src
** Design of the library
@ -142,8 +188,6 @@ cppcheck --addon=cert --enable=all *.c &> cppcheck.out
produced C files should be =xxx.c= and =xxx.h= and the name of the
produced Fortran file should be =xxx.f90=.
Arrays are in uppercase and scalars are in lowercase.
In the names of the variables and functions, only the singular
form is allowed.
@ -240,33 +284,25 @@ cppcheck --addon=cert --enable=all *.c &> cppcheck.out
conversions. These functions are also responsible for allocating
temporary storage, to simplify the use of accelerators.
The high-level functions should be pure, unless the introduction
of non-purity is justified. All the side effects should be made in
the =context= variable.
# TODO : We need an identifier for impure functions
# Suggestion (VJ): using *_unsafe_* for impure functions ?
** Numerical precision
The number of bits of precision required for a function should be
given as an input of low-level computational functions. This input
will be used to define the values of the different thresholds that
might be used to avoid computing unnecessary noise. High-level
functions will use the precision specified in the =context=
variable.
The minimal number of bits of precision required for a function
should be given as an input of low-level computational
functions. This input will be used to define the values of the
different thresholds that might be used to avoid computing
unnecessary noise. High-level functions will use the precision
specified in the =context= variable.
In order to automatize numerical accuracy tests, QMCkl uses
[[https://github.com/verificarlo/verificarlo][Verificarlo]] and
its CI functionality. You can read Verificarlo CI's documentation
at the [[https://github.com/verificarlo/verificarlo/blob/master/doc/06-Postprocessing.md#verificarlo-ci][following link]].
Reading it is advised to understand the remainder of this section.
[[https://github.com/verificarlo/verificarlo][Verificarlo]] and its CI functionality. You can read Verificarlo CI's
documentation at the [[https://github.com/verificarlo/verificarlo/blob/master/doc/06-Postprocessing.md#verificarlo-ci][following link]]. Reading it is advised to
understand the remainder of this section.
To enable support for Verificarlo CI tests when building the
library, use the following configure command :
#+begin_src bash
QMCKL_DEVEL=1 ./configure --prefix=$PWD/_install --enable-silent-rules --enable-maintainer-mode CC=verificarlo-f FC=verificarlo-f --host=x86_64 --enable-vfc_ci
./configure CC=verificarlo-f FC=verificarlo-f --host=x86_64 --enable-vfc_ci
#+end_src
Note that this does require an install of Verificarlo *with
@ -290,7 +326,7 @@ cppcheck --addon=cert --enable=all *.c &> cppcheck.out
- ~qmckl_probe_check_relative~ : place a probe with a relative check. If ~vfc_ci~ is disabled, this will return the result of a relative check (|val - ref| / ref < accuracy target?). If the check fails, true is returned (false otherwise).
If you need more details on these functions or their Fortran
If you need more detail on these functions or their Fortran
interfaces, have a look at the ~tools/qmckl_probes~ files.
Finally, if you need to add a QMCkl kernel to the CI tests

4515
org/qmckl_ao.org
View File

File diff suppressed because it is too large Load Diff

2059
org/qmckl_blas.org
View File

File diff suppressed because it is too large Load Diff

32
org/qmckl_context.org
View File

@ -28,15 +28,21 @@ int main() {
#include "qmckl_error_private_type.h"
#include "qmckl_memory_private_type.h"
#include "qmckl_numprec_private_type.h"
#include "qmckl_point_private_type.h"
#include "qmckl_nucleus_private_type.h"
#include "qmckl_electron_private_type.h"
#include "qmckl_ao_private_type.h"
#include "qmckl_mo_private_type.h"
#include "qmckl_jastrow_private_type.h"
#include "qmckl_determinant_private_type.h"
#include "qmckl_local_energy_private_type.h"
#include "qmckl_point_private_func.h"
#include "qmckl_nucleus_private_func.h"
#include "qmckl_electron_private_func.h"
#include "qmckl_ao_private_func.h"
#include "qmckl_mo_private_func.h"
#include "qmckl_determinant_private_func.h"
#include "qmckl_local_energy_private_func.h"
#+end_src
#+begin_src c :tangle (eval c)
@ -117,16 +123,19 @@ typedef struct qmckl_context_struct {
/* Current date */
uint64_t date;
/* Points */
qmckl_point_struct point;
/* -- Molecular system -- */
qmckl_nucleus_struct nucleus;
qmckl_electron_struct electron;
qmckl_ao_basis_struct ao_basis;
qmckl_mo_basis_struct mo_basis;
qmckl_jastrow_struct jastrow;
qmckl_nucleus_struct nucleus;
qmckl_electron_struct electron;
qmckl_ao_basis_struct ao_basis;
qmckl_mo_basis_struct mo_basis;
qmckl_jastrow_struct jastrow;
qmckl_determinant_struct det;
qmckl_local_energy_struct local_energy;
/* To be implemented:
qmckl_mo_struct mo;
qmckl_determinant_struct det;
,*/
} qmckl_context_struct;
@ -232,6 +241,9 @@ qmckl_context qmckl_context_create() {
ctx->numprec.precision = QMCKL_DEFAULT_PRECISION;
ctx->numprec.range = QMCKL_DEFAULT_RANGE;
rc = qmckl_init_point(context);
assert (rc == QMCKL_SUCCESS);
rc = qmckl_init_electron(context);
assert (rc == QMCKL_SUCCESS);
@ -240,6 +252,12 @@ qmckl_context qmckl_context_create() {
rc = qmckl_init_ao_basis(context);
assert (rc == QMCKL_SUCCESS);
rc = qmckl_init_mo_basis(context);
assert (rc == QMCKL_SUCCESS);
rc = qmckl_init_determinant(context);
assert (rc == QMCKL_SUCCESS);
}
/* Allocate qmckl_memory_struct */

2100
org/qmckl_determinant.org Normal file
View File

File diff suppressed because it is too large Load Diff

135
org/qmckl_distance.org
View File

@ -40,19 +40,21 @@ int main() {
\]
#+NAME: qmckl_distance_sq_args
| qmckl_context | context | in | Global state |
| char | transa | in | Array ~A~ is ~'N'~: Normal, ~'T'~: Transposed |
| char | transb | in | Array ~B~ is ~'N'~: Normal, ~'T'~: Transposed |
| int64_t | m | in | Number of points in the first set |
| int64_t | n | in | Number of points in the second set |
| double | A[][lda] | in | Array containing the $m \times 3$ matrix $A$ |
| int64_t | lda | in | Leading dimension of array ~A~ |
| double | B[][ldb] | in | Array containing the $n \times 3$ matrix $B$ |
| int64_t | ldb | in | Leading dimension of array ~B~ |
| double | C[n][ldc] | out | Array containing the $m \times n$ matrix $C$ |
| int64_t | ldc | in | Leading dimension of array ~C~ |
| Variable | Type | In/Out | Description |
|-----------+------------------+--------+-----------------------------------------------|
| ~context~ | ~qmckl_context~ | in | Global state |
| ~transa~ | ~char~ | in | Array ~A~ is ~'N'~: Normal, ~'T'~: Transposed |
| ~transb~ | ~char~ | in | Array ~B~ is ~'N'~: Normal, ~'T'~: Transposed |
| ~m~ | ~int64_t~ | in | Number of points in the first set |
| ~n~ | ~int64_t~ | in | Number of points in the second set |
| ~A~ | ~double[][lda]~ | in | Array containing the $m \times 3$ matrix $A$ |
| ~lda~ | ~int64_t~ | in | Leading dimension of array ~A~ |
| ~B~ | ~double[][ldb]~ | in | Array containing the $n \times 3$ matrix $B$ |
| ~ldb~ | ~int64_t~ | in | Leading dimension of array ~B~ |
| ~C~ | ~double[n][ldc]~ | out | Array containing the $m \times n$ matrix $C$ |
| ~ldc~ | ~int64_t~ | in | Leading dimension of array ~C~ |
*** Requirements
Requirements:
- ~context~ is not ~QMCKL_NULL_CONTEXT~
- ~m > 0~
@ -66,8 +68,6 @@ int main() {
- ~B~ is allocated with at least $3 \times n \times 8$ bytes
- ~C~ is allocated with at least $m \times n \times 8$ bytes
*** C header
#+CALL: generate_c_header(table=qmckl_distance_sq_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
#+RESULTS:
@ -83,10 +83,9 @@ int main() {
const double* B,
const int64_t ldb,
double* const C,
const int64_t ldc );
const int64_t ldc );
#+end_src
*** Source
#+begin_src f90 :tangle (eval f)
integer function qmckl_distance_sq_f(context, transa, transb, m, n, &
A, LDA, B, LDB, C, LDC) &
@ -222,8 +221,6 @@ end function qmckl_distance_sq_f
This function is more efficient when ~A~ and ~B~ are
transposed.
** C interface :noexport:
#+CALL: generate_c_interface(table=qmckl_distance_sq_args,rettyp=get_value("FRetType"),fname=get_value("Name"))
#+RESULTS:
@ -284,7 +281,6 @@ end function qmckl_distance_sq_f
*** Test :noexport:
#+begin_src f90 :tangle (eval f_test)
integer(qmckl_exit_code) function test_qmckl_distance_sq(context) bind(C)
use qmckl
@ -441,17 +437,19 @@ end function test_qmckl_distance_sq
the leading dimension: ~[n][3]~ in C and ~(3,n)~ in Fortran.
#+NAME: qmckl_distance_args
| qmckl_context | context | in | Global state |
| char | transa | in | Array ~A~ is ~'N'~: Normal, ~'T'~: Transposed |
| char | transb | in | Array ~B~ is ~'N'~: Normal, ~'T'~: Transposed |
| int64_t | m | in | Number of points in the first set |
| int64_t | n | in | Number of points in the second set |
| double | A[][lda] | in | Array containing the $m \times 3$ matrix $A$ |
| int64_t | lda | in | Leading dimension of array ~A~ |
| double | B[][ldb] | in | Array containing the $n \times 3$ matrix $B$ |
| int64_t | ldb | in | Leading dimension of array ~B~ |
| double | C[n][ldc] | out | Array containing the $m \times n$ matrix $C$ |
| int64_t | ldc | in | Leading dimension of array ~C~ |
| Variable | Type | In/Out | Description |
|-----------+------------------+--------+-----------------------------------------------|
| ~context~ | ~qmckl_context~ | in | Global state |
| ~transa~ | ~char~ | in | Array ~A~ is ~'N'~: Normal, ~'T'~: Transposed |
| ~transb~ | ~char~ | in | Array ~B~ is ~'N'~: Normal, ~'T'~: Transposed |
| ~m~ | ~int64_t~ | in | Number of points in the first set |
| ~n~ | ~int64_t~ | in | Number of points in the second set |
| ~A~ | ~double[][lda]~ | in | Array containing the $m \times 3$ matrix $A$ |
| ~lda~ | ~int64_t~ | in | Leading dimension of array ~A~ |
| ~B~ | ~double[][ldb]~ | in | Array containing the $n \times 3$ matrix $B$ |
| ~ldb~ | ~int64_t~ | in | Leading dimension of array ~B~ |
| ~C~ | ~double[n][ldc]~ | out | Array containing the $m \times n$ matrix $C$ |
| ~ldc~ | ~int64_t~ | in | Leading dimension of array ~C~ |
*** Requirements
@ -484,7 +482,7 @@ end function test_qmckl_distance_sq
const double* B,
const int64_t ldb,
double* const C,
const int64_t ldc );
const int64_t ldc );
#+end_src
*** Source
@ -881,18 +879,20 @@ end function test_qmckl_dist
the leading dimension: ~[n][3]~ in C and ~(3,n)~ in Fortran.
#+NAME: qmckl_distance_rescaled_args
| qmckl_context | context | in | Global state |
| char | transa | in | Array ~A~ is ~'N'~: Normal, ~'T'~: Transposed |
| char | transb | in | Array ~B~ is ~'N'~: Normal, ~'T'~: Transposed |
| int64_t | m | in | Number of points in the first set |
| int64_t | n | in | Number of points in the second set |
| double | A[][lda] | in | Array containing the $m \times 3$ matrix $A$ |
| int64_t | lda | in | Leading dimension of array ~A~ |
| double | B[][ldb] | in | Array containing the $n \times 3$ matrix $B$ |
| int64_t | ldb | in | Leading dimension of array ~B~ |
| double | C[n][ldc] | out | Array containing the $m \times n$ matrix $C$ |
| int64_t | ldc | in | Leading dimension of array ~C~ |
| double | rescale_factor_kappa | in | Factor for calculating rescaled distances |
| Variable | Type | In/Out | Description |
|------------------------+------------------+--------+-----------------------------------------------|
| ~context~ | ~qmckl_context~ | in | Global state |
| ~transa~ | ~char~ | in | Array ~A~ is ~'N'~: Normal, ~'T'~: Transposed |
| ~transb~ | ~char~ | in | Array ~B~ is ~'N'~: Normal, ~'T'~: Transposed |
| ~m~ | ~int64_t~ | in | Number of points in the first set |
| ~n~ | ~int64_t~ | in | Number of points in the second set |
| ~A~ | ~double[][lda]~ | in | Array containing the $m \times 3$ matrix $A$ |
| ~lda~ | ~int64_t~ | in | Leading dimension of array ~A~ |
| ~B~ | ~double[][ldb]~ | in | Array containing the $n \times 3$ matrix $B$ |
| ~ldb~ | ~int64_t~ | in | Leading dimension of array ~B~ |
| ~C~ | ~double[n][ldc]~ | out | Array containing the $m \times n$ matrix $C$ |
| ~ldc~ | ~int64_t~ | in | Leading dimension of array ~C~ |
| ~rescale_factor_kappa~ | ~double~ | in | Factor for calculating rescaled distances |
*** Requirements
@ -926,7 +926,7 @@ end function test_qmckl_dist
const int64_t ldb,
double* const C,
const int64_t ldc,
const double rescale_factor_kappa);
const double rescale_factor_kappa );
#+end_src
*** Source
@ -1215,20 +1215,22 @@ end function qmckl_distance_rescaled_f
the leading dimension: ~[n][3]~ in C and ~(3,n)~ in Fortran.
#+NAME: qmckl_distance_rescaled_deriv_e_args
| qmckl_context | context | in | Global state |
| char | transa | in | Array ~A~ is ~'N'~: Normal, ~'T'~: Transposed |
| char | transb | in | Array ~B~ is ~'N'~: Normal, ~'T'~: Transposed |
| int64_t | m | in | Number of points in the first set |
| int64_t | n | in | Number of points in the second set |
| double | A[][lda] | in | Array containing the $m \times 3$ matrix $A$ |
| int64_t | lda | in | Leading dimension of array ~A~ |
| double | B[][ldb] | in | Array containing the $n \times 3$ matrix $B$ |
| int64_t | ldb | in | Leading dimension of array ~B~ |
| double | C[4][n][ldc] | out | Array containing the $4 \times m \times n$ matrix $C$ |
| int64_t | ldc | in | Leading dimension of array ~C~ |
| double | rescale_factor_kappa | in | Factor for calculating rescaled distances derivatives |
| Variable | Type | In/Out | Description |
|------------------------+---------------------+--------+-------------------------------------------------------|
| ~context~ | ~qmckl_context~ | in | Global state |
| ~transa~ | ~char~ | in | Array ~A~ is ~'N'~: Normal, ~'T'~: Transposed |
| ~transb~ | ~char~ | in | Array ~B~ is ~'N'~: Normal, ~'T'~: Transposed |
| ~m~ | ~int64_t~ | in | Number of points in the first set |
| ~n~ | ~int64_t~ | in | Number of points in the second set |
| ~A~ | ~double[][lda]~ | in | Array containing the $m \times 3$ matrix $A$ |
| ~lda~ | ~int64_t~ | in | Leading dimension of array ~A~ |
| ~B~ | ~double[][ldb]~ | in | Array containing the $n \times 3$ matrix $B$ |
| ~ldb~ | ~int64_t~ | in | Leading dimension of array ~B~ |
| ~C~ | ~double[4][n][ldc]~ | out | Array containing the $4 \times m \times n$ matrix $C$ |
| ~ldc~ | ~int64_t~ | in | Leading dimension of array ~C~ |
| ~rescale_factor_kappa~ | ~double~ | in | Factor for calculating rescaled distances derivatives |
*** Requirements
Requirements:
- ~context~ is not ~QMCKL_NULL_CONTEXT~
- ~m > 0~
@ -1241,9 +1243,7 @@ end function qmckl_distance_rescaled_f
- ~A~ is allocated with at least $3 \times m \times 8$ bytes
- ~B~ is allocated with at least $3 \times n \times 8$ bytes
- ~C~ is allocated with at least $4 \times m \times n \times 8$ bytes
*** C header
#+CALL: generate_c_header(table=qmckl_distance_rescaled_deriv_e_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
#+RESULTS:
@ -1260,10 +1260,9 @@ end function qmckl_distance_rescaled_f
const int64_t ldb,
double* const C,
const int64_t ldc,
const double rescale_factor_kappa);
const double rescale_factor_kappa );
#+end_src
*** Source
#+begin_src f90 :tangle (eval f)
integer function qmckl_distance_rescaled_deriv_e_f(context, transa, transb, m, n, &
A, LDA, B, LDB, C, LDC, rescale_factor_kappa) &
@ -1450,13 +1449,9 @@ integer function qmckl_distance_rescaled_deriv_e_f(context, transa, transb, m, n
end function qmckl_distance_rescaled_deriv_e_f
#+end_src
*** Performance
This function is more efficient when ~A~ and ~B~ are transposed.
** C interface :noexport:
#+CALL: generate_c_interface(table=qmckl_distance_rescaled_deriv_e_args,rettyp=get_value("FRetType"),fname=get_value("Name"))
#+CALL: generate_c_interface(table=qmckl_distance_rescaled_deriv_e_args,fname=get_value("Name"))
#+RESULTS:
#+begin_src f90 :tangle (eval f) :comments org :exports none
@ -1476,7 +1471,7 @@ end function qmckl_distance_rescaled_deriv_e_f
integer (c_int64_t) , intent(in) , value :: lda
real (c_double ) , intent(in) :: B(ldb,*)
integer (c_int64_t) , intent(in) , value :: ldb
real (c_double ) , intent(out) :: C(4,ldc,n)
real (c_double ) , intent(out) :: C(ldc,n,4)
integer (c_int64_t) , intent(in) , value :: ldc
real (c_double ) , intent(in) , value :: rescale_factor_kappa
@ -1508,7 +1503,7 @@ end function qmckl_distance_rescaled_deriv_e_f
integer (c_int64_t) , intent(in) , value :: lda
real (c_double ) , intent(in) :: B(ldb,*)
integer (c_int64_t) , intent(in) , value :: ldb
real (c_double ) , intent(out) :: C(4,ldc,n)
real (c_double ) , intent(out) :: C(ldc,n,4)
integer (c_int64_t) , intent(in) , value :: ldc
real (c_double ) , intent(in) , value :: rescale_factor_kappa

1014
org/qmckl_electron.org
View File

File diff suppressed because it is too large Load Diff

1581
org/qmckl_jastrow.org
View File

File diff suppressed because it is too large Load Diff

1692
org/qmckl_local_energy.org Normal file
View File

File diff suppressed because it is too large Load Diff

2
org/qmckl_memory.org
View File

@ -194,7 +194,7 @@ b[2] = 3; assert(b[2] == 3);
case some important information has been stored related to memory
allocation and needs to be updated.
#+begin_src c :tangle (eval h_func)
#+begin_src c :tangle (eval h_private_func)
qmckl_exit_code qmckl_free(qmckl_context context,
void * const ptr);
#+end_src

238
org/qmckl_mo.org
View File

@ -2,7 +2,7 @@
#+SETUPFILE: ../tools/theme.setup
#+INCLUDE: ../tools/lib.org
The molecular orbitals (MOs) are defined in the basis of AOs along with a AO to MO
The molecular orbitals (MOs) are defined in the basis of AOs along with a AO to MO
coefficient matrix \[C\]. Using these coefficients (e.g. from Hartree Fock SCF method)
the MOs are defined as follows:
@ -83,12 +83,13 @@ int main() {
The following arrays are stored in the context:
|---------------+--------------------+----------------------|
| ~mo_num~ | | Number of MOs |
| ~coefficient~ | ~[mo_num][ao_num]~ | Orbital coefficients |
Computed data:
|---------------+-------------------------+----------------------------------------------------------------------------------------|
|---------------+-------------------------+----------------------------------------------------------------------------------------|
| ~mo_vgl~ | ~[5][elec_num][mo_num]~ | Value, gradients, Laplacian of the MOs at electron positions |
@ -99,7 +100,7 @@ int main() {
#+begin_src c :comments org :tangle (eval h_private_type)
typedef struct qmckl_mo_basis_struct {
int64_t mo_num;
int64_t mo_num;
double * coefficient;
double * mo_vgl;
@ -117,6 +118,26 @@ typedef struct qmckl_mo_basis_struct {
Some values are initialized by default, and are not concerned by
this mechanism.
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_exit_code qmckl_init_mo_basis(qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c)
qmckl_exit_code qmckl_init_mo_basis(qmckl_context context) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return false;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
ctx->mo_basis.uninitialized = (1 << 2) - 1;
return QMCKL_SUCCESS;
}
#+end_src
** Access functions
#+begin_src c :comments org :tangle (eval h_func) :exports none
@ -205,7 +226,7 @@ qmckl_get_mo_basis_coefficient (const qmckl_context context,
assert (ctx->mo_basis.coefficient != NULL);
memcpy(coefficient, ctx->mo_basis.coefficient,
ctx->ao_basis.ao_num * ctx->mo_basis.mo_num * sizeof(int64_t));
ctx->ao_basis.ao_num * ctx->mo_basis.mo_num * sizeof(double));
return QMCKL_SUCCESS;
}
@ -258,10 +279,9 @@ qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
ctx->mo_basis.uninitialized &= ~mask;
ctx->mo_basis.provided = (ctx->mo_basis.uninitialized == 0);
if (ctx->mo_basis.provided) {
qmckl_exit_code rc_ = qmckl_finalize_basis(context);
qmckl_exit_code rc_ = qmckl_finalize_mo_basis(context);
if (rc_ != QMCKL_SUCCESS) return rc_;
}
}
return QMCKL_SUCCESS;
#+end_src
@ -318,6 +338,28 @@ qmckl_exit_code qmckl_set_mo_basis_coefficient(qmckl_context context, const dou
When the basis set is completely entered, other data structures are
computed to accelerate the calculations.
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_finalize_mo_basis(qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_finalize_mo_basis(qmckl_context context) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT,
"qmckl_finalize_mo_basis",
NULL);
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
qmckl_exit_code rc = QMCKL_SUCCESS;
return rc;
}
#+end_src
* Computation
** Computation of MOs
@ -330,7 +372,7 @@ qmckl_exit_code qmckl_get_mo_basis_vgl(qmckl_context context, double* const mo_v
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_mo_basis_vgl(qmckl_context context, double* const mo_vgl) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
@ -377,6 +419,7 @@ qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context);
qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context)
{
qmckl_exit_code rc;
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
@ -391,13 +434,21 @@ qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context)
NULL);
}
rc = qmckl_provide_ao_vgl(context);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_ao_basis",
NULL);
}
if(!(ctx->electron.provided)) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_electron",
NULL);
}
if (!ctx->mo_basis.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
@ -424,7 +475,6 @@ qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context)
ctx->mo_basis.mo_vgl = mo_vgl;
}
qmckl_exit_code rc;
rc = qmckl_compute_mo_basis_vgl(context,
ctx->ao_basis.ao_num,
ctx->mo_basis.mo_num,
@ -442,7 +492,7 @@ qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context)
return QMCKL_SUCCESS;
}
#+end_src
*** Compute
:PROPERTIES:
:Name: qmckl_compute_mo_basis_vgl
@ -450,7 +500,7 @@ qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context)
:FRetType: qmckl_exit_code
:END:
#+NAME: qmckl_mo_basis_gaussian_vgl_args
#+NAME: qmckl_mo_basis_vgl_args
| ~qmckl_context~ | ~context~ | in | Global state |
| ~int64_t~ | ~ao_num~ | in | Number of AOs |
| ~int64_t~ | ~mo_num~ | in | Number of MOs |
@ -458,7 +508,8 @@ qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context)
| ~double~ | ~coef_normalized[mo_num][ao_num]~ | in | AO to MO transformation matrix |
| ~double~ | ~ao_vgl[5][elec_num][ao_num]~ | in | Value, gradients and Laplacian of the AOs |
| ~double~ | ~mo_vgl[5][elec_num][mo_num]~ | out | Value, gradients and Laplacian of the MOs |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_mo_basis_vgl_f(context, &
ao_num, mo_num, elec_num, &
@ -472,86 +523,71 @@ integer function qmckl_compute_mo_basis_vgl_f(context, &
double precision , intent(in) :: ao_vgl(ao_num,elec_num,5)
double precision , intent(in) :: coef_normalized(ao_num,mo_num)
double precision , intent(out) :: mo_vgl(mo_num,elec_num,5)
logical*8 :: TransA, TransB
character :: TransA, TransB
double precision,dimension(:,:),allocatable :: mo_vgl_big
double precision,dimension(:,:),allocatable :: ao_vgl_big
!double precision,dimension(:,:),allocatable :: coef_trans
!double precision,dimension(:),allocatable :: coef_all
double precision :: alpha, beta
integer :: info_qmckl_dgemm_value
integer :: info_qmckl_dgemm_Gx
integer :: info_qmckl_dgemm_Gy
integer :: info_qmckl_dgemm_Gz
integer :: info_qmckl_dgemm_lap
integer*8 :: M, N, K, LDA, LDB, LDC, i,j
integer*8 :: M, N, K, LDA, LDB, LDC, i,j, idx
integer*8 :: inucl, iprim, iwalk, ielec, ishell
double precision :: x, y, z, two_a, ar2, r2, v, cutoff
TransA = .False.
TransB = .False.
allocate(mo_vgl_big(mo_num,elec_num*5))
allocate(ao_vgl_big(ao_num,elec_num*5))
!allocate(coef_all(mo_num*ao_num))
!allocate(coef_trans(mo_num,ao_num))
TransA = 'T'
TransB = 'N'
alpha = 1.0d0
beta = 0.0d0
info = QMCKL_SUCCESS
info_qmckl_dgemm_value = QMCKL_SUCCESS
info_qmckl_dgemm_Gx = QMCKL_SUCCESS
info_qmckl_dgemm_Gy = QMCKL_SUCCESS
info_qmckl_dgemm_Gz = QMCKL_SUCCESS
info_qmckl_dgemm_lap = QMCKL_SUCCESS
! Don't compute exponentials when the result will be almost zero.
! TODO : Use numerical precision here
cutoff = -dlog(1.d-15)
M = 1_8
N = mo_num * 1_8
M = mo_num
N = elec_num*5
K = ao_num * 1_8
LDA = M
LDB = K
LDC = M
LDA = size(coef_normalized,1)
idx = 0
!do j = 1,ao_num
!do i = 1,mo_num
! idx = idx + 1
! coef_all(idx) = coef_normalized(i,j)
!end do
!end do
!idx = 0
!do j = 1,mo_num
!do i = 1,ao_num
! idx = idx + 1
! coef_trans(j,i) = coef_all(idx)
!end do
!end do
do ielec = 1, elec_num
! Value
info_qmckl_dgemm_value = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, &
ao_vgl(:, ielec, 1), LDA, &
coef_normalized(1:ao_num,1:mo_num),size(coef_normalized,1) * 1_8, &
beta, &
mo_vgl(:,ielec,1),LDC)
! Grad_x
info_qmckl_dgemm_Gx = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, &
ao_vgl(:, ielec, 2), LDA, &
coef_normalized(1:ao_num,1:mo_num),size(coef_normalized,1) * 1_8, &
beta, &
mo_vgl(:,ielec,2),LDC)
! Grad_y
info_qmckl_dgemm_Gy = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, &
ao_vgl(:, ielec, 3), LDA, &
coef_normalized(1:ao_num,1:mo_num),size(coef_normalized,1) * 1_8, &
beta, &
mo_vgl(:,ielec,3),LDC)
! Grad_z
info_qmckl_dgemm_Gz = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, &
ao_vgl(:, ielec, 4), LDA, &
coef_normalized(1:ao_num,1:mo_num),size(coef_normalized,1) * 1_8, &
beta, &
mo_vgl(:,ielec,4),LDC)
! Lapl_z
info_qmckl_dgemm_lap = qmckl_dgemm(context, TransA, TransB, M, N, K, alpha, &
ao_vgl(:, ielec, 5), LDA, &
coef_normalized(1:ao_num,1:mo_num),size(coef_normalized,1) * 1_8, &
beta, &
mo_vgl(:,ielec,5),LDC)
end do
ao_vgl_big = reshape(ao_vgl(:, :, :),(/ao_num, elec_num*5_8/))
LDB = size(ao_vgl_big,1)
LDC = size(mo_vgl_big,1)
info = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, &
coef_normalized,size(coef_normalized,1)*1_8, &
ao_vgl_big, size(ao_vgl_big,1)*1_8, &
beta, &
mo_vgl_big,LDC)
mo_vgl = reshape(mo_vgl_big,(/mo_num,elec_num,5_8/))
deallocate(mo_vgl_big)
deallocate(ao_vgl_big)
if(info_qmckl_dgemm_value .eq. QMCKL_SUCCESS .and. &
info_qmckl_dgemm_Gx .eq. QMCKL_SUCCESS .and. &
info_qmckl_dgemm_Gy .eq. QMCKL_SUCCESS .and. &
info_qmckl_dgemm_Gz .eq. QMCKL_SUCCESS .and. &
info_qmckl_dgemm_lap .eq. QMCKL_SUCCESS ) then
info = QMCKL_SUCCESS
else
info = QMCKL_FAILURE
end if
end function qmckl_compute_mo_basis_vgl_f
#+end_src
#+CALL: generate_c_header(table=qmckl_mo_basis_gaussian_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_vgl"))
#+CALL: generate_c_header(table=qmckl_mo_basis_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_vgl"))
#+RESULTS:
#+begin_src c :tangle (eval h_func) :comments org
@ -562,11 +598,11 @@ end function qmckl_compute_mo_basis_vgl_f
const int64_t elec_num,
const double* coef_normalized,
const double* ao_vgl,
double* const mo_vgl );
double* const mo_vgl );
#+end_src
#+CALL: generate_c_interface(table=qmckl_mo_basis_gaussian_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_vgl"))
#+CALL: generate_c_interface(table=qmckl_mo_basis_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_vgl"))
#+RESULTS:
#+begin_src f90 :tangle (eval f) :comments org :exports none
@ -598,7 +634,7 @@ end function qmckl_compute_mo_basis_vgl_f
import numpy as np
def f(a,x,y):
return np.sum( [c * np.exp( -b*(np.linalg.norm(x-y))**2) for b,c in a] )
return np.sum( [c * np.exp( -b*(np.linalg.norm(x-y))**2) for b,c in a] )
def df(a,x,y,n):
h0 = 1.e-6
@ -670,7 +706,7 @@ const int64_t nucl_num = chbrclf_nucl_num;
const double* nucl_charge = chbrclf_charge;
const double* nucl_coord = &(chbrclf_nucl_coord[0][0]);
rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num);
rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num);
assert (rc == QMCKL_SUCCESS);
rc = qmckl_set_electron_walk_num (context, walk_num);
@ -678,16 +714,16 @@ assert (rc == QMCKL_SUCCESS);
assert(qmckl_electron_provided(context));
rc = qmckl_set_electron_coord (context, 'N', elec_coord);
rc = qmckl_set_electron_coord (context, 'N', elec_coord, walk_num*elec_num*3);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_num (context, nucl_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]));
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), nucl_num*3);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_charge(context, nucl_charge);
rc = qmckl_set_nucleus_charge(context, nucl_charge, nucl_num);
assert(rc == QMCKL_SUCCESS);
assert(qmckl_nucleus_provided(context));
@ -719,45 +755,45 @@ rc = qmckl_set_ao_basis_prim_num (context, chbrclf_prim_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_nucleus_index (context, nucleus_index);
rc = qmckl_set_ao_basis_nucleus_index (context, nucleus_index, nucl_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_nucleus_shell_num (context, nucleus_shell_num);
rc = qmckl_set_ao_basis_nucleus_shell_num (context, nucleus_shell_num, nucl_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_shell_ang_mom (context, shell_ang_mom);
rc = qmckl_set_ao_basis_shell_ang_mom (context, shell_ang_mom, chbrclf_shell_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_shell_factor (context, shell_factor);
rc = qmckl_set_ao_basis_shell_factor (context, shell_factor, chbrclf_shell_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_shell_prim_num (context, shell_prim_num);
rc = qmckl_set_ao_basis_shell_prim_num (context, shell_prim_num, chbrclf_shell_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_shell_prim_index (context, shell_prim_index);
rc = qmckl_set_ao_basis_shell_prim_index (context, shell_prim_index, chbrclf_shell_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_exponent (context, exponent);
rc = qmckl_set_ao_basis_exponent (context, exponent, chbrclf_prim_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_coefficient (context, coefficient);
rc = qmckl_set_ao_basis_coefficient (context, coefficient, chbrclf_prim_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_prim_factor (context, prim_factor);
rc = qmckl_set_ao_basis_prim_factor (context, prim_factor, chbrclf_prim_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_ao_basis_ao_num(context, chbrclf_ao_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_ao_basis_ao_factor (context, ao_factor);
rc = qmckl_set_ao_basis_ao_factor (context, ao_factor, chbrclf_ao_num);
assert(rc == QMCKL_SUCCESS);
assert(qmckl_ao_basis_provided(context));
@ -765,7 +801,8 @@ assert(qmckl_ao_basis_provided(context));
double ao_vgl[5][walk_num][elec_num][chbrclf_ao_num];
rc = qmckl_get_ao_vgl(context, &(ao_vgl[0][0][0][0]));
rc = qmckl_get_ao_basis_ao_vgl(context, &(ao_vgl[0][0][0][0]),
(int64_t) 5*walk_num*elec_num*chbrclf_ao_num);
assert (rc == QMCKL_SUCCESS);
/* Set up MO data */
@ -785,10 +822,10 @@ rc = qmckl_get_mo_basis_vgl(context, &(mo_vgl[0][0][0]));
assert (rc == QMCKL_SUCCESS);
// Test overlap of MO
//double point_x[100];
//double point_y[100];
//double point_z[100];
//int32_t npoints=100;
//double point_x[10];
//double point_y[10];
//double point_z[10];
//int32_t npoints=10;
//// obtain points
//double dr = 20./(npoints-1);
//double dr3 = dr*dr*dr;
@ -802,15 +839,16 @@ assert (rc == QMCKL_SUCCESS);
//double ovlmo1 = 0.0;
//// Calculate overlap
//for (int i=0;i<npoints;++i) {
// printf(".");
// fflush(stdout);
// for (int j=0;j<npoints;++j) {
// printf(" .. ");
// for (int k=0;k<npoints;++k) {
// printf(" . ");
// // Set point
// elec_coord[0] = point_x[i];
// elec_coord[1] = point_y[j];
// elec_coord[2] = point_z[k];
// rc = qmckl_set_electron_coord (context, 'N', elec_coord);
// rc = qmckl_set_electron_coord (context, 'N', elec_coord);
// assert(rc == QMCKL_SUCCESS);
//
// // Calculate value of MO (1st electron)
@ -838,8 +876,8 @@ printf(" mo_vgl mo_vgl[1][26][223] %25.15e\n", mo_vgl[1][2][3]);
printf(" mo_vgl mo_vgl[0][26][224] %25.15e\n", mo_vgl[0][2][3]);
printf(" mo_vgl mo_vgl[1][26][224] %25.15e\n", mo_vgl[1][2][3]);
printf("\n");
}
}
#+end_src
* End of files :noexport:

540
org/qmckl_nucleus.org
View File

@ -67,35 +67,42 @@ int main() {
The following data stored in the context:
| ~uninitialized~ | int32_t | Keeps bit set for uninitialized data |
| ~num~ | int64_t | Total number of nuclei |
| ~provided~ | bool | If true, ~nucleus~ is valid |
| ~charge~ | double[num] | Nuclear charges |
| ~coord~ | double[3][num] | Nuclear coordinates, in transposed format |
| ~nn_distance~ | double[num][num] | Nucleus-nucleus distances |
| ~nn_distance_date~ | int64_t | Date when Nucleus-nucleus distances were computed |
| ~nn_distance_rescaled~ | double[num][num] | Nucleus-nucleus rescaled distances |
| ~nn_distance_rescaled_date~ | int64_t | Date when Nucleus-nucleus rescaled distances were computed |
| ~repulsion~ | double | Nuclear repulsion energy |
| ~repulsion_date~ | int64_t | Date when the nuclear repulsion energy was computed |
| ~rescale_factor_kappa~ | double | The distance scaling factor |
|------------------------+--------------+-------------------------------------------|
| ~uninitialized~ | int32_t | Keeps bit set for uninitialized data |
| ~num~ | int64_t | Total number of nuclei |
| ~provided~ | bool | If true, ~nucleus~ is valid |
| ~charge~ | qmckl_vector | Nuclear charges |
| ~coord~ | qmckl_matrix | Nuclear coordinates, in transposed format |
| ~coord_date~ | int64_t | Nuclear coordinates, date if modified |
| ~rescale_factor_kappa~ | double | The distance scaling factor |
Computed data:
|-----------------------------+--------------+------------------------------------------------------------|
| ~nn_distance~ | qmckl_matrix | Nucleus-nucleus distances |
| ~nn_distance_date~ | int64_t | Date when Nucleus-nucleus distances were computed |
| ~nn_distance_rescaled~ | qmckl_matrix | Nucleus-nucleus rescaled distances |
| ~nn_distance_rescaled_date~ | int64_t | Date when Nucleus-nucleus rescaled distances were computed |
| ~repulsion~ | double | Nuclear repulsion energy |
| ~repulsion_date~ | int64_t | Date when the nuclear repulsion energy was computed |
** Data structure
#+begin_src c :comments org :tangle (eval h_private_type)
typedef struct qmckl_nucleus_struct {
int64_t num;
int64_t repulsion_date;
int64_t nn_distance_date;
int64_t nn_distance_rescaled_date;
double* coord;
double* charge;
double* nn_distance;
double* nn_distance_rescaled;
double repulsion;
double rescale_factor_kappa;
int32_t uninitialized;
bool provided;
int64_t num;
int64_t repulsion_date;
int64_t nn_distance_date;
int64_t nn_distance_rescaled_date;
int64_t coord_date;
qmckl_vector charge;
qmckl_matrix coord;
qmckl_matrix nn_distance;
qmckl_matrix nn_distance_rescaled;
double repulsion;
double rescale_factor_kappa;
int32_t uninitialized;
bool provided;
} qmckl_nucleus_struct;
#+end_src
@ -130,18 +137,7 @@ qmckl_exit_code qmckl_init_nucleus(qmckl_context context) {
}
#+end_src
** Access functions
#+begin_src c :comments org :tangle (eval h_func) :exports none
qmckl_exit_code qmckl_get_nucleus_num (const qmckl_context context, int64_t* const num);
qmckl_exit_code qmckl_get_nucleus_charge (const qmckl_context context, double* const charge);
qmckl_exit_code qmckl_get_nucleus_coord (const qmckl_context context, const char transp, double* const coord);
qmckl_exit_code qmckl_get_nucleus_rescale_factor (const qmckl_context context, double* const rescale_factor_kappa);
#+end_src
#+NAME:post
#+begin_src c :exports none
if ( (ctx->nucleus.uninitialized & mask) != 0) {
@ -149,6 +145,13 @@ if ( (ctx->nucleus.uninitialized & mask) != 0) {
}
#+end_src
#+begin_src c :comments org :tangle (eval h_func) :exports none
qmckl_exit_code
qmckl_get_nucleus_num(const qmckl_context context,
int64_t* const num);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_nucleus_num (const qmckl_context context, int64_t* const num) {
@ -182,10 +185,22 @@ qmckl_get_nucleus_num (const qmckl_context context, int64_t* const num) {
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src c :comments org :tangle (eval h_func) :exports none
qmckl_exit_code
qmckl_get_nucleus_charge (const qmckl_context context, double* const charge) {
qmckl_get_nucleus_charge(const qmckl_context context,
double* const charge,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_nucleus_charge (const qmckl_context context,
double* const charge,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_INVALID_CONTEXT;
@ -210,16 +225,31 @@ qmckl_get_nucleus_charge (const qmckl_context context, double* const charge) {
"nucleus data is not provided");
}
assert (ctx->nucleus.charge != NULL);
assert (ctx->nucleus.charge.data != NULL);
int64_t nucl_num = ctx->nucleus.num;
if (ctx->nucleus.num > size_max) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_nucleus_charge",
"Array too small");
}
memcpy(charge, ctx->nucleus.charge, nucl_num*sizeof(double));
qmckl_exit_code rc;
rc = qmckl_double_of_vector(context, ctx->nucleus.charge, charge, size_max);
return QMCKL_SUCCESS;
return rc;
}
#+end_src
#+begin_src c :comments org :tangle (eval h_func) :exports none
qmckl_exit_code
qmckl_get_nucleus_rescale_factor(const qmckl_context context,
double* const rescale_factor_kappa);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_nucleus_rescale_factor (const qmckl_context context,
double* const rescale_factor_kappa)
@ -245,10 +275,24 @@ qmckl_get_nucleus_rescale_factor (const qmckl_context context,
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src c :comments org :tangle (eval h_func) :exports none
qmckl_exit_code
qmckl_get_nucleus_coord (const qmckl_context context, const char transp, double* const coord) {
qmckl_get_nucleus_coord(const qmckl_context context,
const char transp,
double* const coord,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_nucleus_coord (const qmckl_context context,
const char transp,
double* const coord,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_INVALID_CONTEXT;
@ -280,25 +324,21 @@ qmckl_get_nucleus_coord (const qmckl_context context, const char transp, double*
"nucleus data is not provided");
}
int64_t nucl_num = ctx->nucleus.num;
assert (ctx->nucleus.coord.data != NULL);
assert (ctx->nucleus.coord != NULL);
qmckl_exit_code rc;
if (transp == 'N') {
qmckl_exit_code rc;
rc = qmckl_transpose(context, nucl_num, 3,
ctx->nucleus.coord, nucl_num,
coord, 3);
qmckl_matrix At = qmckl_matrix_alloc(context, 3, ctx->nucleus.coord.size[0]);
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);
} else {
memcpy(coord, ctx->nucleus.coord, 3*nucl_num*sizeof(double));
rc = qmckl_double_of_matrix(context, ctx->nucleus.coord, coord, size_max);
}
return QMCKL_SUCCESS;
return rc;
}
#+end_src
@ -325,51 +365,6 @@ bool qmckl_nucleus_provided(const qmckl_context context) {
** Initialization functions
To set the data relative to the nuclei in the context, the
following functions need to be called.
#+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code qmckl_set_nucleus_num (qmckl_context context, const int64_t num);
qmckl_exit_code qmckl_set_nucleus_charge (qmckl_context context, const double* charge);
qmckl_exit_code qmckl_set_nucleus_coord (qmckl_context context, const char transp, const double* coord);
qmckl_exit_code qmckl_set_nucleus_rescale_factor (qmckl_context context, const double rescale_factor_kappa);
#+end_src
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_set_nucleus_num(context, num) &
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 :: num
end function
end interface
interface
integer(c_int32_t) function qmckl_set_nucleus_charge(context, charge) &
bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
real (c_double) , intent(in) :: charge(*)
end function
end interface
interface
integer(c_int32_t) function qmckl_set_nucleus_coord(context, transp, coord) &
bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
character(c_char) , intent(in) , value :: transp
real (c_double) , intent(in) :: coord(*)
end function
end interface
#+end_src
#+NAME:pre2
#+begin_src c :exports none
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
@ -387,12 +382,22 @@ ctx->nucleus.provided = (ctx->nucleus.uninitialized == 0);
return QMCKL_SUCCESS;
#+end_src
To set the data relative to the nuclei in the context, the
following functions need to be called.
To set the number of nuclei, use
#+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_set_nucleus_num(qmckl_context context,
const int64_t num);
#+end_src
Sets the number of nuclei.
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_set_nucleus_num(qmckl_context context, const int64_t num) {
qmckl_set_nucleus_num(qmckl_context context,
const int64_t num)
{
<<pre2>>
if (num <= 0) {
@ -410,11 +415,35 @@ qmckl_set_nucleus_num(qmckl_context context, const int64_t num) {
}
#+end_src
The following function sets the nuclear charges of all the atoms.
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_set_nucleus_num(context, num) &
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 :: num
end function
end interface
#+end_src
#+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_set_nucleus_charge(qmckl_context context,
const double* charge,
const int64_t size_max);
#+end_src
Sets the nuclear charges of all the atoms.
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_set_nucleus_charge(qmckl_context context, const double* charge) {
qmckl_set_nucleus_charge(qmckl_context context,
const double* charge,
const int64_t size_max)
{
<<pre2>>
if (charge == NULL) {
@ -432,34 +461,126 @@ qmckl_set_nucleus_charge(qmckl_context context, const double* charge) {
rc = qmckl_get_nucleus_num(context, &num);
if (rc != QMCKL_SUCCESS) return rc;
if (ctx->nucleus.charge != NULL) {
qmckl_free(context, ctx->nucleus.charge);
ctx->nucleus.charge= NULL;
}
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = num*sizeof(double);
assert (ctx->nucleus.charge == NULL);
ctx->nucleus.charge = (double*) qmckl_malloc(context, mem_info);
if (ctx->nucleus.charge == NULL) {
if (num > size_max) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
QMCKL_INVALID_ARG_3,
"qmckl_set_nucleus_charge",
NULL);
"Array too small");
}
ctx->nucleus.charge= memcpy(ctx->nucleus.charge, charge, num*sizeof(double));
assert (ctx->nucleus.charge != NULL);
ctx->nucleus.charge = qmckl_vector_alloc(context, num);
rc = qmckl_vector_of_double(context, charge, num, &(ctx->nucleus.charge));
<<post2>>
}
#+end_src
The following function sets the rescale parameter for the nuclear distances.
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_set_nucleus_charge(context, charge, size_max) &
bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
real (c_double) , intent(in) :: charge(*)
integer (c_int64_t) , intent(in) , value :: size_max
end function
end interface
#+end_src
#+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_set_nucleus_coord(qmckl_context context,
const char transp,
const double* coord,
const int64_t size_max);
#+end_src
Sets the nuclear coordinates of all the atoms. The coordinates
are be given in atomic units.
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_set_nucleus_rescale_factor(qmckl_context context, const double rescale_factor_kappa) {
qmckl_set_nucleus_coord(qmckl_context context,
const char transp,
const double* coord,
const int64_t size_max)
{
<<pre2>>
qmckl_exit_code rc;
int32_t mask = 1 << 2;
const int64_t nucl_num = (int64_t) ctx->nucleus.num;
if (ctx->nucleus.coord.data != NULL) {
rc = qmckl_matrix_free(context, ctx->nucleus.coord);
if (rc != QMCKL_SUCCESS) return rc;
}
ctx->nucleus.coord = qmckl_matrix_alloc(context, nucl_num, 3);
if (ctx->nucleus.coord.data == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_set_nucleus_coord",
NULL);
}
if (size_max < 3*nucl_num) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_4,
"qmckl_set_nucleus_coord",
"Array too small");
}
if (transp == 'N') {
qmckl_matrix At;
At = qmckl_matrix_alloc(context, 3, nucl_num);
rc = qmckl_matrix_of_double(context, coord, 3*nucl_num, &At);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_transpose(context, At, ctx->nucleus.coord);
} else {
rc = qmckl_matrix_of_double(context, coord, nucl_num*3, &(ctx->nucleus.coord));
}
if (rc != QMCKL_SUCCESS) return rc;
<<post2>>
}
#+end_src
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_set_nucleus_coord(context, transp, coord, size_max) &
bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
character(c_char) , intent(in) , value :: transp
real (c_double) , intent(in) :: coord(*)
integer (c_int64_t) , intent(in) , value :: size_max
end function
end interface
#+end_src
#+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_set_nucleus_rescale_factor(qmckl_context context,
const double kappa);
#+end_src
Sets the rescale parameter for the nuclear distances.
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_set_nucleus_rescale_factor(qmckl_context context,
const double rescale_factor_kappa)
{
<<pre2>>
if (rescale_factor_kappa <= 0.0) {
@ -475,50 +596,17 @@ qmckl_set_nucleus_rescale_factor(qmckl_context context, const double rescale_fac
}
#+end_src
The following function sets the nuclear coordinates of all the
atoms. The coordinates should be given in atomic units.
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_set_nucleus_coord(qmckl_context context, const char transp, const double* coord) {
<<pre2>>
int64_t nucl_num = (int64_t) 0;
qmckl_exit_code rc;
int32_t mask = 1 << 2;
rc = qmckl_get_nucleus_num(context, &nucl_num);
if (rc != QMCKL_SUCCESS) return rc;
if (ctx->nucleus.coord != NULL) {
qmckl_free(context, ctx->nucleus.coord);
ctx->nucleus.coord = NULL;
}
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = 3*nucl_num*sizeof(double);
assert(ctx->nucleus.coord == NULL);
ctx->nucleus.coord = (double*) qmckl_malloc(context, mem_info);
if (coord == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_set_nucleus_coord",
NULL);
}
if (transp == 'N') {
rc = qmckl_transpose(context, 3, nucl_num,
coord, 3,
ctx->nucleus.coord, nucl_num);
if (rc != QMCKL_SUCCESS) return rc;
} else {
memcpy(ctx->nucleus.coord, coord, 3*nucl_num*sizeof(double));
}
<<post2>>
}
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_set_rescale_factor(context, kappa) &
bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
real (c_double) , intent(in) , value :: kappa
end function
end interface
#+end_src
** Test
@ -563,15 +651,15 @@ assert(k == nucl_rescale_factor_kappa);
double nucl_coord2[3*nucl_num];
rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2);
rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2, 3*nucl_num);
assert(rc == QMCKL_NOT_PROVIDED);
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]));
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), 3*nucl_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_nucleus_provided(context));
rc = qmckl_get_nucleus_coord (context, 'N', nucl_coord2);
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 ; i<nucl_num ; ++i) {
@ -579,7 +667,7 @@ for (size_t k=0 ; k<3 ; ++k) {
}
}
rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2);
rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2, 3*nucl_num);
assert(rc == QMCKL_SUCCESS);
for (size_t i=0 ; i<3*nucl_num ; ++i) {
assert( nucl_coord[i] == nucl_coord2[i] );
@ -587,13 +675,13 @@ for (size_t i=0 ; i<3*nucl_num ; ++i) {
double nucl_charge2[nucl_num];
rc = qmckl_get_nucleus_charge(context, nucl_charge2);
rc = qmckl_get_nucleus_charge(context, nucl_charge2, nucl_num);
assert(rc == QMCKL_NOT_PROVIDED);
rc = qmckl_set_nucleus_charge(context, nucl_charge);
rc = qmckl_set_nucleus_charge(context, nucl_charge, nucl_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_get_nucleus_charge(context, nucl_charge2);
rc = qmckl_get_nucleus_charge(context, nucl_charge2, nucl_num);
assert(rc == QMCKL_SUCCESS);
for (size_t i=0 ; i<nucl_num ; ++i) {
assert( nucl_charge[i] == nucl_charge2[i] );
@ -616,11 +704,17 @@ assert(qmckl_nucleus_provided(context));
*** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_nucleus_nn_distance(qmckl_context context, double* distance);
qmckl_exit_code
qmckl_get_nucleus_nn_distance(qmckl_context context,
double* distance,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_nucleus_nn_distance(qmckl_context context, double* distance)
qmckl_exit_code
qmckl_get_nucleus_nn_distance(qmckl_context context,
double* distance,
const int64_t size_max)
{
/* Check input parameters */
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
@ -633,22 +727,29 @@ qmckl_exit_code qmckl_get_nucleus_nn_distance(qmckl_context context, double* dis
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
size_t sze = ctx->nucleus.num * ctx->nucleus.num;
memcpy(distance, ctx->nucleus.nn_distance, sze * sizeof(double));
const int64_t sze = ctx->nucleus.num * ctx->nucleus.num;
if (sze > size_max) {
return qmckl_failwith(context,
QMCKL_INVALID_ARG_3,
"qmckl_get_nucleus_nn_distance",
"Array too small");
}
rc = qmckl_double_of_matrix(context, ctx->nucleus.nn_distance, distance, size_max);
return QMCKL_SUCCESS;
return rc;
}
#+end_src
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_get_nucleus_nn_distance(context, distance) &
integer(c_int32_t) function qmckl_get_nucleus_nn_distance(context, distance, size_max) &
bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
real (c_double ) , intent(out) :: distance(*)
integer (c_int64_t) , intent(in) , value :: size_max
end function
end interface
#+end_src
@ -673,26 +774,23 @@ qmckl_exit_code qmckl_provide_nn_distance(qmckl_context context)
if (!ctx->nucleus.provided) return QMCKL_NOT_PROVIDED;
/* Allocate array */
if (ctx->nucleus.nn_distance == NULL) {
if (ctx->nucleus.nn_distance.data == NULL) {
ctx->nucleus.nn_distance =
qmckl_matrix_alloc(context, ctx->nucleus.num, ctx->nucleus.num);
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->nucleus.num * ctx->nucleus.num * sizeof(double);
double* nn_distance = (double*) qmckl_malloc(context, mem_info);
if (nn_distance == NULL) {
if (ctx->nucleus.nn_distance.data == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_nn_distance",
NULL);
}
ctx->nucleus.nn_distance = nn_distance;
}
qmckl_exit_code rc =
qmckl_compute_nn_distance(context,
ctx->nucleus.num,
ctx->nucleus.coord,
ctx->nucleus.nn_distance);
ctx->nucleus.coord.data,
ctx->nucleus.nn_distance.data);
if (rc != QMCKL_SUCCESS) {
return rc;
}
@ -783,24 +881,29 @@ qmckl_exit_code qmckl_compute_nn_distance (
assert(qmckl_nucleus_provided(context));
double distance[nucl_num*nucl_num];
rc = qmckl_get_nucleus_nn_distance(context, distance);
rc = qmckl_get_nucleus_nn_distance(context, distance, nucl_num*nucl_num);
assert(distance[0] == 0.);
assert(distance[1] == distance[nucl_num]);
assert(fabs(distance[1]-2.070304721365169) < 1.e-12);
#+end_src
** Nucleus-nucleus rescaled distances
*** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_nucleus_nn_distance_rescaled(qmckl_context context, double* distance_rescaled);
qmckl_exit_code
qmckl_get_nucleus_nn_distance_rescaled(qmckl_context context,
double* 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_nucleus_nn_distance_rescaled(qmckl_context context, double* distance_rescaled)
qmckl_exit_code
qmckl_get_nucleus_nn_distance_rescaled(qmckl_context context,
double* distance_rescaled,
const int64_t size_max)
{
/* Check input parameters */
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
@ -813,13 +916,35 @@ qmckl_exit_code qmckl_get_nucleus_nn_distance_rescaled(qmckl_context context, do
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
size_t sze = ctx->nucleus.num * ctx->nucleus.num;
memcpy(distance_rescaled, ctx->nucleus.nn_distance_rescaled, sze * sizeof(double));
return QMCKL_SUCCESS;
const int64_t sze = ctx->nucleus.num * ctx->nucleus.num;
if (sze > size_max) {
return qmckl_failwith(context,
QMCKL_INVALID_ARG_3,
"qmckl_get_nucleus_nn_distance_rescaled",
"Array too small");
}
rc = qmckl_double_of_matrix(context,
ctx->nucleus.nn_distance_rescaled,
distance_rescaled,
size_max);
return rc;
}
#+end_src
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_get_nucleus_nn_distance_rescaled(context, distance_rescaled, size_max) &
bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
real (c_double ) , intent(out) :: distance_rescaled(*)
integer (c_int64_t) , intent(in) , value :: size_max
end function
end interface
#+end_src
*** Provide :noexport:
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
@ -840,27 +965,24 @@ qmckl_exit_code qmckl_provide_nn_distance_rescaled(qmckl_context context)
if (!ctx->nucleus.provided) return QMCKL_NOT_PROVIDED;
/* Allocate array */
if (ctx->nucleus.nn_distance_rescaled == NULL) {
if (ctx->nucleus.nn_distance_rescaled.data == NULL) {
ctx->nucleus.nn_distance_rescaled =
qmckl_matrix_alloc(context, ctx->nucleus.num, ctx->nucleus.num);
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->nucleus.num * ctx->nucleus.num * sizeof(double);
double* nn_distance_rescaled = (double*) qmckl_malloc(context, mem_info);
if (nn_distance_rescaled == NULL) {
if (ctx->nucleus.nn_distance_rescaled.data == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_nn_distance_rescaled",
NULL);
}
ctx->nucleus.nn_distance_rescaled = nn_distance_rescaled;
}
qmckl_exit_code rc =
qmckl_compute_nn_distance_rescaled(context,
ctx->nucleus.num,
ctx->nucleus.rescale_factor_kappa,
ctx->nucleus.coord,
ctx->nucleus.nn_distance_rescaled);
ctx->nucleus.coord.data,
ctx->nucleus.nn_distance_rescaled.data);
if (rc != QMCKL_SUCCESS) {
return rc;
}
@ -955,7 +1077,7 @@ qmckl_exit_code qmckl_compute_nn_distance_rescaled (
//assert(qmckl_nucleus_provided(context));
//
//double distance[nucl_num*nucl_num];
//rc = qmckl_get_nucleus_nn_distance(context, distance);
//rc = qmckl_get_nucleus_nn_distance(context, distance, nucl_num*nucl_num);
//assert(distance[0] == 0.);
//assert(distance[1] == distance[nucl_num]);
//assert(fabs(distance[1]-2.070304721365169) < 1.e-12);
@ -971,11 +1093,11 @@ qmckl_exit_code qmckl_compute_nn_distance_rescaled (
*** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_nucleus_repulsion(qmckl_context context, double* energy);
qmckl_exit_code qmckl_get_nucleus_repulsion(qmckl_context context, double* const energy);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_nucleus_repulsion(qmckl_context context, double* energy)
qmckl_exit_code qmckl_get_nucleus_repulsion(qmckl_context context, double* const energy)
{
/* Check input parameters */
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
@ -1033,8 +1155,8 @@ qmckl_exit_code qmckl_provide_nucleus_repulsion(qmckl_context context)
rc = qmckl_compute_nucleus_repulsion(context,
ctx->nucleus.num,
ctx->nucleus.charge,
ctx->nucleus.nn_distance,
ctx->nucleus.charge.data,
ctx->nucleus.nn_distance.data,
&(ctx->nucleus.repulsion));
if (rc != QMCKL_SUCCESS) {
return rc;
@ -1083,7 +1205,9 @@ integer function qmckl_compute_nucleus_repulsion_f(context, nucl_num, charge, nn
energy = 0.d0
do j=2, nucl_num
do i=1, j-1
energy = energy + charge(i) * charge(j) / nn_distance(i,j)
if (dabs(nn_distance(i,j)) > 1e-5) then
energy = energy + charge(i) * charge(j) / nn_distance(i,j)
endif
end do
end do

457
org/qmckl_point.org Normal file
View File

@ -0,0 +1,457 @@
#+TITLE: Point
#+SETUPFILE: ../tools/theme.setup
#+INCLUDE: ../tools/lib.org
This data structure contains cartesian coordinates where the functions
will be evaluated. For DFT codes these may be the integration grid
points. For QMC codes, these are the electron coordinates of all the
walkers.
* Headers :noexport:
#+begin_src elisp :noexport :results none
(org-babel-lob-ingest "../tools/lib.org")
#+end_src
#+begin_src c :tangle (eval h_private_type)
#ifndef QMCKL_POINT_HPT
#define QMCKL_POINT_HPT
#include <stdbool.h>
#include "qmckl_blas_private_type.h"
#+end_src
#+begin_src c :tangle (eval h_private_func)
#ifndef QMCKL_POINT_HPF
#define QMCKL_POINT_HPF
#include "qmckl_blas_private_type.h"
#include "qmckl_blas_private_func.h"
#+end_src
#+begin_src c :tangle (eval c_test) :noweb yes
#include "qmckl.h"
#include <assert.h>
#include <math.h>
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include "chbrclf.h"
#include "qmckl_blas_private_type.h"
#include "qmckl_blas_private_func.h"
int main() {
qmckl_context context;
context = qmckl_context_create();
#+end_src
#+begin_src c :tangle (eval c)
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#ifdef HAVE_STDINT_H
#include <stdint.h>
#elif HAVE_INTTYPES_H
#include <inttypes.h>
#endif
#include <stdlib.h>
#include <string.h>
#include <stdbool.h>
#include <assert.h>
#include <math.h>
#include <stdio.h>
#include "qmckl.h"
#include "qmckl_context_private_type.h"
#include "qmckl_memory_private_type.h"
#include "qmckl_blas_private_type.h"
#include "qmckl_point_private_func.h"
#include "qmckl_memory_private_func.h"
#include "qmckl_blas_private_func.h"
#+end_src
* Context
The following data stored in the context:
| Variable | Type | Description |
|----------+----------------+-------------------------------------------|
| ~num~ | ~int64_t~ | Total number of points |
| ~date~ | ~uint64_t~ | Last modification date of the coordinates |
| ~coord~ | ~qmckl_matrix~ | ~num~ \times 3 matrix3 |
We consider that the matrix is stored 'transposed' and 'normal'
corresponds to the 3 \times ~num~ matrix.
** Data structure
#+begin_src c :comments org :tangle (eval h_private_type)
typedef struct qmckl_point_struct {
int64_t num;
uint64_t date;
qmckl_matrix coord;
} qmckl_point_struct;
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_exit_code qmckl_init_point(qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c)
qmckl_exit_code qmckl_init_point(qmckl_context context) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return false;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
memset(&(ctx->point), 0, sizeof(qmckl_point_struct));
return QMCKL_SUCCESS;
}
#+end_src
** Access functions
Access functions return ~QMCKL_SUCCESS~ when the data has been
successfully retrieved. They return ~QMCKL_INVALID_CONTEXT~ when
the context is not a valid context. If the function returns
successfully, the variable pointed by the pointer given in argument
contains the requested data. Otherwise, this variable is untouched.
*** Number of points
#+begin_src c :comments org :tangle (eval h_func) :exports none
qmckl_exit_code qmckl_get_point_num (const qmckl_context context, int64_t* const num);
#+end_src
Returns the number of points stored in the context.
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_point_num (const qmckl_context context, int64_t* const num) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_INVALID_CONTEXT;
}
if (num == NULL) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
"qmckl_get_point_num",
"num is a null pointer");
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
assert (ctx->point.num >= (int64_t) 0);
,*num = ctx->point.num;
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src f90 :comments org :tangle (eval fh_func) :noweb yes
interface
integer(c_int32_t) function qmckl_get_point_num(context, num) bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(out) :: num
end function
end interface
#+end_src
*** Point coordinates
#+begin_src c :comments org :tangle (eval h_func) :exports none
qmckl_exit_code qmckl_get_point(const qmckl_context context,
const char transp,
double* const coord,
const int64_t size_max);
#+end_src
Returns the point coordinates as sequences of (x,y,z).
The pointer is assumed to point on a memory block of size
~size_max~ \ge ~3 * point_num~.
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_point(const qmckl_context context,
const char transp,
double* const coord,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_INVALID_CONTEXT;
}
if (coord == NULL) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
"qmckl_get_point_coord",
"coord is a null pointer");
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
int64_t point_num = ctx->point.num;
if (point_num == 0) return QMCKL_NOT_PROVIDED;
assert (ctx->point.coord.data != NULL);
if (size_max < 3*point_num) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_point_coord",
"size_max too small");
}
qmckl_exit_code rc;
if (transp == 'N') {
qmckl_matrix At = qmckl_matrix_alloc( context, 3, point_num);
rc = qmckl_transpose( context, ctx->point.coord, At );
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_double_of_matrix( context, At, coord, size_max);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_matrix_free(context, At);
} else {
rc = qmckl_double_of_matrix( context, ctx->point.coord, coord, size_max);
}
if (rc != QMCKL_SUCCESS) return rc;
return rc;
}
#+end_src
#+begin_src f90 :comments org :tangle (eval fh_func) :noweb yes
interface
integer(c_int32_t) function qmckl_get_point(context, transp, coord, size_max) bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
character(c_char) , intent(in) , value :: transp
real (c_double ) , intent(out) :: coord(*)
integer (c_int64_t) , intent(in) :: size_max
end function
end interface
#+end_src
** Initialization functions
When the data is set in the context, if the arrays are large
enough, we overwrite the data contained in them.
To set the data relative to the points in the context, one of the
following functions need to be called.
#+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code qmckl_set_point (qmckl_context context,
const char transp,
const double* coord,
const int64_t num);
#+end_src
Copy a sequence of ~num~ points $(x,y,z)$ into the context.
#+begin_src c :comments org :tangle (eval c) :noweb yes
qmckl_exit_code
qmckl_set_point (qmckl_context context,
const char transp,
const double* coord,
const int64_t num)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
if (transp != 'N' && transp != 'T') {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
"qmckl_set_point",
"transp should be 'N' or 'T'");
}
if (coord == NULL) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_set_point",
"coord is a NULL pointer");
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
qmckl_exit_code rc;
if (ctx->point.num < num) {
if (ctx->point.coord.data != NULL) {
rc = qmckl_matrix_free(context, ctx->point.coord);
assert (rc == QMCKL_SUCCESS);
}
ctx->point.coord = qmckl_matrix_alloc(context, num, 3);
if (ctx->point.coord.data == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_set_point",
NULL);
}
};
ctx->point.num = num;
if (transp == 'T') {
memcpy(ctx->point.coord.data, coord, 3*num*sizeof(double));
} else {
for (int64_t i=0 ; i<num ; ++i) {
qmckl_mat(ctx->point.coord, i, 0) = coord[3*i ];
qmckl_mat(ctx->point.coord, i, 1) = coord[3*i+1];
qmckl_mat(ctx->point.coord, i, 2) = coord[3*i+2];
}
}
/* Increment the date of the context */
ctx->date += 1UL;
ctx->point.date = ctx->date;
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src f90 :comments org :tangle (eval fh_func) :noweb yes
interface
integer(c_int32_t) function qmckl_set_point(context, &
transp, coord, size_max) bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
character(c_char) , intent(in) , value :: transp
real (c_double ) , intent(in) :: coord(*)
integer (c_int64_t) , intent(in) , value :: size_max
end function
end interface
#+end_src
** Test
#+begin_src c :tangle (eval c_test)
/* Reference input data */
int64_t point_num = chbrclf_elec_num;
const double* coord = &(chbrclf_elec_coord[0][0][0]);
/* --- */
qmckl_exit_code rc;
double coord2[point_num*3];
double coord3[point_num*3];
rc = qmckl_get_point (context, 'N', coord2, (point_num*3));
assert(rc == QMCKL_NOT_PROVIDED);
rc = qmckl_set_point (context, 'N', coord, point_num);
assert(rc == QMCKL_SUCCESS);
int64_t n;
rc = qmckl_get_point_num (context, &n);
assert(rc == QMCKL_SUCCESS);
assert(n == point_num);
rc = qmckl_get_point (context, 'N', coord2, (point_num*3));
assert(rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i<3*point_num ; ++i) {
assert( coord[i] == coord2[i] );
}
rc = qmckl_get_point (context, 'T', coord2, (point_num*3));
assert(rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i<point_num ; ++i) {
assert( coord[3*i+0] == coord2[i] );
assert( coord[3*i+1] == coord2[i+point_num] );
assert( coord[3*i+2] == coord2[i+point_num*2] );
}
rc = qmckl_set_point (context, 'T', coord2, point_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_get_point (context, 'N', coord3, (point_num*3));
assert(rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i<3*point_num ; ++i) {
assert( coord[i] == coord3[i] );
}
#+end_src
* End of files :noexport:
#+begin_src c :tangle (eval h_private_type)
#endif
#+end_src
#+begin_src c :tangle (eval h_private_func)
#endif
#+end_src
*** Test
#+begin_src c :tangle (eval c_test)
if (qmckl_context_destroy(context) != QMCKL_SUCCESS)
return QMCKL_FAILURE;
return 0;
}
#+end_src
*** Compute file names
#+begin_src emacs-lisp
; The following is required to compute the file names
(setq pwd (file-name-directory buffer-file-name))
(setq name (file-name-nondirectory (substring buffer-file-name 0 -4)))
(setq f (concat pwd name "_f.f90"))
(setq fh (concat pwd name "_fh.f90"))
(setq c (concat pwd name ".c"))
(setq h (concat name ".h"))
(setq h_private (concat name "_private.h"))
(setq c_test (concat pwd "test_" name ".c"))
(setq f_test (concat pwd "test_" name "_f.f90"))
; Minted
(require 'ox-latex)
(setq org-latex-listings 'minted)
(add-to-list 'org-latex-packages-alist '("" "listings"))
(add-to-list 'org-latex-packages-alist '("" "color"))
#+end_src
#+RESULTS:
| | color |
| | listings |
# -*- mode: org -*-
# vim: syntax=c

166
org/qmckl_sherman_morrison_woodbury.org
View File

File diff suppressed because one or more lines are too long

117839
org/qmckl_tests.org
View File

File diff suppressed because it is too large Load Diff

311
org/qmckl_trexio.org
View File

@ -2,7 +2,7 @@
#+SETUPFILE: ../tools/theme.setup
#+INCLUDE: ../tools/lib.org
The [[https://github.com/trex-coe/trexio][TREXIO library]] enables easy and efficient input/output of wave
The [[https://trex-coe.github.io/trexio/trex.html][TREXIO library]] enables easy and efficient input/output of wave
function parameters. In this section we provide high-level functions
to prepare the context by reading the required data from a TREXIO file.
@ -200,7 +200,7 @@ qmckl_trexio_read_nucleus_X(qmckl_context context, trexio_t* const file)
assert (nucl_charge != NULL);
rcio = trexio_read_nucleus_charge_64(file, nucl_charge);
rcio = trexio_read_safe_nucleus_charge_64(file, nucl_charge, nucleus_num);
if (rcio != TREXIO_SUCCESS) {
return qmckl_failwith( context,
QMCKL_FAILURE,
@ -208,7 +208,7 @@ qmckl_trexio_read_nucleus_X(qmckl_context context, trexio_t* const file)
trexio_string_of_error(rcio));
}
rc = qmckl_set_nucleus_charge(context, nucl_charge);
rc = qmckl_set_nucleus_charge(context, nucl_charge, nucleus_num);
qmckl_free(context, nucl_charge);
nucl_charge = NULL;
@ -240,7 +240,7 @@ qmckl_trexio_read_nucleus_X(qmckl_context context, trexio_t* const file)
assert (nucl_coord != NULL);
rcio = trexio_read_nucleus_coord_64(file, nucl_coord);
rcio = trexio_read_safe_nucleus_coord_64(file, nucl_coord, 3*nucleus_num);
if (rcio != TREXIO_SUCCESS) {
return qmckl_failwith( context,
QMCKL_FAILURE,
@ -248,7 +248,7 @@ qmckl_trexio_read_nucleus_X(qmckl_context context, trexio_t* const file)
trexio_string_of_error(rcio));
}
rc = qmckl_set_nucleus_coord(context, 'N', nucl_coord);
rc = qmckl_set_nucleus_coord(context, 'N', nucl_coord, 3*nucleus_num);
qmckl_free(context, nucl_coord);
nucl_coord = NULL;
@ -322,7 +322,7 @@ qmckl_trexio_read_ao_X(qmckl_context context, trexio_t* const file)
#+begin_src c :tangle (eval c)
int64_t shell_num = 0L;
rcio = trexio_read_basis_num_64(file, &shell_num);
rcio = trexio_read_basis_shell_num_64(file, &shell_num);
if (rcio != TREXIO_SUCCESS) {
return qmckl_failwith( context,
QMCKL_FAILURE,
@ -385,8 +385,9 @@ qmckl_trexio_read_ao_X(qmckl_context context, trexio_t* const file)
#+begin_src c :tangle (eval c)
{
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = nucleus_num * sizeof(int64_t);
/* Allocate array for data */
mem_info.size = nucleus_num * sizeof(int64_t);
int64_t* nucleus_index = (int64_t*) qmckl_malloc(context, mem_info);
if (nucleus_index == NULL) {
@ -398,15 +399,57 @@ qmckl_trexio_read_ao_X(qmckl_context context, trexio_t* const file)
assert (nucleus_index != NULL);
rcio = trexio_read_basis_nucleus_index_64(file, nucleus_index);
/* Allocate temporary array */
mem_info.size = shell_num * sizeof(int64_t);
int64_t* tmp_array = (int64_t*) qmckl_malloc(context, mem_info);
if (tmp_array == NULL) {
qmckl_free(context, nucleus_index);
nucleus_index = NULL;
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_trexio_read_basis_nucleus_index_X",
NULL);
}
assert (tmp_array != NULL);
/* Read in the temporary array */
rcio = trexio_read_safe_basis_nucleus_index_64(file, tmp_array, shell_num);
if (rcio != TREXIO_SUCCESS) {
qmckl_free(context, tmp_array);
tmp_array = NULL;
qmckl_free(context, nucleus_index);
nucleus_index = NULL;
return qmckl_failwith( context,
QMCKL_FAILURE,
"trexio_read_basis_nucleus_index",
trexio_string_of_error(rcio));
}
rc = qmckl_set_ao_basis_nucleus_index(context, nucleus_index);
/* Reformat data */
rc = qmckl_set_ao_basis_nucleus_index(context, nucleus_index, nucleus_num);
for (int i=shell_num-1 ; i>=0 ; --i) {
const int k = tmp_array[i];
if (k < 0 || k >= nucleus_num) {
qmckl_free(context, tmp_array);
tmp_array = NULL;
qmckl_free(context, nucleus_index);
nucleus_index = NULL;
return qmckl_failwith( context,
QMCKL_FAILURE,
"trexio_read_basis_nucleus_index",
"Irrelevant data in TREXIO file");
}
nucleus_index[k] = i;
}
qmckl_free(context, tmp_array);
tmp_array = NULL;
/* Store data */
rc = qmckl_set_ao_basis_nucleus_index(context, nucleus_index, shell_num);
qmckl_free(context, nucleus_index);
nucleus_index = NULL;
@ -421,8 +464,9 @@ qmckl_trexio_read_ao_X(qmckl_context context, trexio_t* const file)
#+begin_src c :tangle (eval c)
{
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = nucleus_num * sizeof(int64_t);
/* Allocate array for data */
mem_info.size = nucleus_num * sizeof(int64_t);
int64_t* nucleus_shell_num = (int64_t*) qmckl_malloc(context, mem_info);
if (nucleus_shell_num == NULL) {
@ -434,15 +478,60 @@ qmckl_trexio_read_ao_X(qmckl_context context, trexio_t* const file)
assert (nucleus_shell_num != NULL);
rcio = trexio_read_basis_nucleus_shell_num_64(file, nucleus_shell_num);
/* Allocate temporary array */
mem_info.size = shell_num * sizeof(int64_t);
int64_t* tmp_array = (int64_t*) qmckl_malloc(context, mem_info);
if (tmp_array == NULL) {
qmckl_free(context, nucleus_shell_num);
nucleus_shell_num = NULL;
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_trexio_read_basis_nucleus_shell_num_X",
NULL);
}
assert (tmp_array != NULL);
/* Read in the temporary array */
rcio = trexio_read_safe_basis_nucleus_index_64(file, tmp_array, shell_num);
if (rcio != TREXIO_SUCCESS) {
qmckl_free(context, tmp_array);
tmp_array = NULL;
qmckl_free(context, nucleus_shell_num);
nucleus_shell_num = NULL;
return qmckl_failwith( context,
QMCKL_FAILURE,
"trexio_read_basis_nucleus_shell_num",
trexio_string_of_error(rcio));
}
rc = qmckl_set_ao_basis_nucleus_shell_num(context, nucleus_shell_num);
/* Reformat data */
for (int i=0 ; i<nucleus_num ; ++i) {
nucleus_shell_num[i] = 0;
}
for (int i=0 ; i<shell_num ; ++i) {
const int k = tmp_array[i];
if (k < 0 || k >= nucleus_num) {
qmckl_free(context, tmp_array);
tmp_array = NULL;
qmckl_free(context, nucleus_shell_num);
nucleus_shell_num = NULL;
return qmckl_failwith( context,
QMCKL_FAILURE,
"trexio_read_basis_nucleus_shell_num",
"Irrelevant data in TREXIO file");
}
nucleus_shell_num[k] += 1;
}
qmckl_free(context, tmp_array);
tmp_array = NULL;
/* Store data */
rc = qmckl_set_ao_basis_nucleus_shell_num(context, nucleus_shell_num, shell_num);
qmckl_free(context, nucleus_shell_num);
nucleus_shell_num = NULL;
@ -457,6 +546,8 @@ qmckl_trexio_read_ao_X(qmckl_context context, trexio_t* const file)
#+begin_src c :tangle (eval c)
{
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
/* Allocate array for data */
mem_info.size = shell_num * sizeof(int32_t);
int32_t* shell_ang_mom = (int32_t*) qmckl_malloc(context, mem_info);
@ -470,15 +561,19 @@ qmckl_trexio_read_ao_X(qmckl_context context, trexio_t* const file)
assert (shell_ang_mom != NULL);
rcio = trexio_read_basis_shell_ang_mom_32(file, shell_ang_mom);
/* Read data */
rcio = trexio_read_safe_basis_shell_ang_mom_32(file, shell_ang_mom, shell_num);
if (rcio != TREXIO_SUCCESS) {
qmckl_free(context, shell_ang_mom);
shell_ang_mom = NULL;
return qmckl_failwith( context,
QMCKL_FAILURE,
"trexio_read_basis_shell_ang_mom",
trexio_string_of_error(rcio));
}
rc = qmckl_set_ao_basis_shell_ang_mom(context, shell_ang_mom);
/* Store data */
rc = qmckl_set_ao_basis_shell_ang_mom(context, shell_ang_mom, shell_num);
qmckl_free(context, shell_ang_mom);
shell_ang_mom = NULL;
@ -493,6 +588,8 @@ qmckl_trexio_read_ao_X(qmckl_context context, trexio_t* const file)
#+begin_src c :tangle (eval c)
{
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
/* Allocate array for data */
mem_info.size = shell_num * sizeof(int64_t);
int64_t* shell_prim_num = (int64_t*) qmckl_malloc(context, mem_info);
@ -506,15 +603,58 @@ qmckl_trexio_read_ao_X(qmckl_context context, trexio_t* const file)
assert (shell_prim_num != NULL);
rcio = trexio_read_basis_shell_prim_num_64(file, shell_prim_num);
/* Allocate temporary array */
mem_info.size = prim_num * sizeof(int64_t);
int64_t* tmp_array = (int64_t*) qmckl_malloc(context, mem_info);
if (tmp_array == NULL) {
qmckl_free(context, shell_prim_num);
shell_prim_num = NULL;
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_trexio_read_basis_shell_prim_num_X",
NULL);
}
assert (tmp_array != NULL);
/* Read data */
rcio = trexio_read_safe_basis_shell_index_64 (file, tmp_array, prim_num);
if (rcio != TREXIO_SUCCESS) {
qmckl_free(context, shell_prim_num);
shell_prim_num = NULL;
qmckl_free(context, tmp_array);
tmp_array = NULL;
return qmckl_failwith( context,
QMCKL_FAILURE,
"trexio_read_basis_shell_prim_num",
trexio_string_of_error(rcio));
}
rc = qmckl_set_ao_basis_shell_prim_num(context, shell_prim_num);
/* Reformat data */
for (int i=0 ; i<shell_num ; ++i) {
shell_prim_num[i] = 0;
}
for (int i=0 ; i<prim_num ; ++i) {
const int k = tmp_array[i];
if (k < 0 || k >= shell_num) {
qmckl_free(context, tmp_array);
qmckl_free(context, shell_prim_num);
return qmckl_failwith( context,
QMCKL_FAILURE,
"trexio_read_basis_shell_prim_num",
"Irrelevant data in TREXIO file");
}
shell_prim_num[k] += 1;
}
qmckl_free(context, tmp_array);
tmp_array = NULL;
/* Store data */
rc = qmckl_set_ao_basis_shell_prim_num(context, shell_prim_num, shell_num);
qmckl_free(context, shell_prim_num);
shell_prim_num = NULL;
@ -529,6 +669,8 @@ qmckl_trexio_read_ao_X(qmckl_context context, trexio_t* const file)
#+begin_src c :tangle (eval c)
{
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
/* Allocate array for data */
mem_info.size = shell_num * sizeof(int64_t);
int64_t* shell_prim_index = (int64_t*) qmckl_malloc(context, mem_info);
@ -542,15 +684,54 @@ qmckl_trexio_read_ao_X(qmckl_context context, trexio_t* const file)
assert (shell_prim_index != NULL);
rcio = trexio_read_basis_shell_prim_index_64(file, shell_prim_index);
/* Allocate temporary array */
mem_info.size = prim_num * sizeof(int64_t);
int64_t* tmp_array = (int64_t*) qmckl_malloc(context, mem_info);
if (tmp_array == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_trexio_read_basis_shell_prim_index_X",
NULL);
}
assert (tmp_array != NULL);
/* Read data */
rcio = trexio_read_safe_basis_shell_index_64(file, tmp_array, prim_num);
if (rcio != TREXIO_SUCCESS) {
qmckl_free(context, shell_prim_index);
shell_prim_index = NULL;
qmckl_free(context, tmp_array);
tmp_array = NULL;
return qmckl_failwith( context,
QMCKL_FAILURE,
"trexio_read_basis_shell_prim_index",
trexio_string_of_error(rcio));
}
rc = qmckl_set_ao_basis_shell_prim_index(context, shell_prim_index);
/* Reformat data */
for (int i=prim_num-1 ; i>=0 ; --i) {
const int k = tmp_array[i];
if (k < 0 || k >= shell_num) {
qmckl_free(context, tmp_array);
tmp_array = NULL;
qmckl_free(context, shell_prim_index);
shell_prim_index = NULL;
return qmckl_failwith( context,
QMCKL_FAILURE,
"trexio_read_basis_shell_prim_index",
"Irrelevant data in TREXIO file");
}
shell_prim_index[k] = i;
}
qmckl_free(context, tmp_array);
tmp_array = NULL;
/* Store data */
rc = qmckl_set_ao_basis_shell_prim_index(context, shell_prim_index, shell_num);
qmckl_free(context, shell_prim_index);
shell_prim_index = NULL;
@ -565,6 +746,8 @@ qmckl_trexio_read_ao_X(qmckl_context context, trexio_t* const file)
#+begin_src c :tangle (eval c)
{
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
/* Allocate array for data */
mem_info.size = shell_num * sizeof(double);
double* shell_factor = (double*) qmckl_malloc(context, mem_info);
@ -578,15 +761,19 @@ qmckl_trexio_read_ao_X(qmckl_context context, trexio_t* const file)
assert (shell_factor != NULL);
rcio = trexio_read_basis_shell_factor_64(file, shell_factor);
/* Read data */
rcio = trexio_read_safe_basis_shell_factor_64(file, shell_factor, shell_num);
if (rcio != TREXIO_SUCCESS) {
qmckl_free(context, shell_factor);
shell_factor = NULL;
return qmckl_failwith( context,
QMCKL_FAILURE,
"trexio_read_basis_shell_factor",
trexio_string_of_error(rcio));
}
rc = qmckl_set_ao_basis_shell_factor(context, shell_factor);
/* Store data */
rc = qmckl_set_ao_basis_shell_factor(context, shell_factor, shell_num);
qmckl_free(context, shell_factor);
shell_factor = NULL;
@ -601,6 +788,8 @@ qmckl_trexio_read_ao_X(qmckl_context context, trexio_t* const file)
#+begin_src c :tangle (eval c)
{
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
/* Allocate array for data */
mem_info.size = prim_num * sizeof(double);
double* exponent = (double*) qmckl_malloc(context, mem_info);
@ -614,15 +803,19 @@ qmckl_trexio_read_ao_X(qmckl_context context, trexio_t* const file)
assert (exponent != NULL);
rcio = trexio_read_basis_exponent_64(file, exponent);
/* Read data */
rcio = trexio_read_safe_basis_exponent_64(file, exponent, prim_num);
if (rcio != TREXIO_SUCCESS) {
qmckl_free(context, exponent);
exponent = NULL;
return qmckl_failwith( context,
QMCKL_FAILURE,
"trexio_read_basis_exponent",
trexio_string_of_error(rcio));
}
rc = qmckl_set_ao_basis_exponent(context, exponent);
/* Store data */
rc = qmckl_set_ao_basis_exponent(context, exponent, prim_num);
qmckl_free(context, exponent);
exponent = NULL;
@ -637,6 +830,8 @@ qmckl_trexio_read_ao_X(qmckl_context context, trexio_t* const file)
#+begin_src c :tangle (eval c)
{
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
/* Allocate array for data */
mem_info.size = prim_num * sizeof(double);
double* coefficient = (double*) qmckl_malloc(context, mem_info);
@ -650,15 +845,19 @@ qmckl_trexio_read_ao_X(qmckl_context context, trexio_t* const file)
assert (coefficient != NULL);
rcio = trexio_read_basis_coefficient_64(file, coefficient);
/* Read data */
rcio = trexio_read_safe_basis_coefficient_64(file, coefficient, prim_num);
if (rcio != TREXIO_SUCCESS) {
qmckl_free(context, coefficient);
coefficient = NULL;
return qmckl_failwith( context,
QMCKL_FAILURE,
"trexio_read_basis_coefficient",
trexio_string_of_error(rcio));
}
rc = qmckl_set_ao_basis_coefficient(context, coefficient);
/* Store data */
rc = qmckl_set_ao_basis_coefficient(context, coefficient, prim_num);
qmckl_free(context, coefficient);
coefficient = NULL;
@ -673,6 +872,8 @@ qmckl_trexio_read_ao_X(qmckl_context context, trexio_t* const file)
#+begin_src c :tangle (eval c)
{
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
/* Allocate array for data */
mem_info.size = prim_num * sizeof(double);
double* prim_factor = (double*) qmckl_malloc(context, mem_info);
@ -686,15 +887,19 @@ qmckl_trexio_read_ao_X(qmckl_context context, trexio_t* const file)
assert (prim_factor != NULL);
rcio = trexio_read_basis_prim_factor_64(file, prim_factor);
/* Read data */
rcio = trexio_read_safe_basis_prim_factor_64(file, prim_factor, prim_num);
if (rcio != TREXIO_SUCCESS) {
qmckl_free(context, prim_factor);
prim_factor = NULL;
return qmckl_failwith( context,
QMCKL_FAILURE,
"trexio_read_basis_prim_factor",
trexio_string_of_error(rcio));
}
rc = qmckl_set_ao_basis_prim_factor(context, prim_factor);
/* Read data */
rc = qmckl_set_ao_basis_prim_factor(context, prim_factor, prim_num);
qmckl_free(context, prim_factor);
prim_factor = NULL;
@ -704,6 +909,48 @@ qmckl_trexio_read_ao_X(qmckl_context context, trexio_t* const file)
}
#+end_src
*** AO Normalization
#+begin_src c :tangle (eval c)
{
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
/* Allocate array for data */
mem_info.size = ao_num * sizeof(double);
double* ao_normalization = (double*) qmckl_malloc(context, mem_info);
if (ao_normalization == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_trexio_read_ao_normalization_X",
NULL);
}
assert (ao_normalization != NULL);
/* Read data */
rcio = trexio_read_safe_ao_normalization_64(file, ao_normalization, ao_num);
if (rcio != TREXIO_SUCCESS) {
qmckl_free(context, ao_normalization);
ao_normalization = NULL;
return qmckl_failwith( context,
QMCKL_FAILURE,
"trexio_read_ao_normalization",
trexio_string_of_error(rcio));
}
/* Store data */
rc = qmckl_set_ao_basis_ao_factor(context, ao_normalization, ao_num);
qmckl_free(context, ao_normalization);
ao_normalization = NULL;
if (rc != QMCKL_SUCCESS)
return rc;
}
#+end_src
#+begin_src c :tangle (eval c)
@ -711,7 +958,6 @@ qmckl_trexio_read_ao_X(qmckl_context context, trexio_t* const file)
}
#endif
#+end_src
** Molecular orbitals
In this section we read the MO coefficients.
@ -880,6 +1126,7 @@ qmckl_trexio_read(const qmckl_context context, const char* file_name)
#+begin_src c :tangle (eval c_test)
#ifdef HAVE_TREXIO
#define walk_num 2
qmckl_exit_code rc;
char fname[256];
@ -892,6 +1139,8 @@ char message[256];
strncpy(fname, QMCKL_TEST_DIR,255);
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);
if (rc != QMCKL_SUCCESS) {
@ -932,7 +1181,7 @@ assert (nucl_num == chbrclf_nucl_num);
printf("Nuclear charges\n");
double * charge = (double*) malloc (nucl_num * sizeof(double));
rc = qmckl_get_nucleus_charge(context, charge);
rc = qmckl_get_nucleus_charge(context, charge, nucl_num);
assert (rc == QMCKL_SUCCESS);
for (int i=0 ; i<nucl_num ; i++) {
assert (charge[i] == chbrclf_charge[i]);
@ -942,7 +1191,7 @@ charge = NULL;
printf("Nuclear coordinates\n");
double * coord = (double*) malloc (nucl_num * 3 * sizeof(double));
rc = qmckl_get_nucleus_coord(context, 'T', coord);
rc = qmckl_get_nucleus_coord(context, 'T', coord, 3*nucl_num);
assert (rc == QMCKL_SUCCESS);
int k=0;
for (int j=0 ; j<3 ; ++j) {
@ -1077,7 +1326,9 @@ double * mo_coef = (double*) malloc (ao_num * mo_num * sizeof(double));
rc = qmckl_get_mo_basis_coefficient(context, mo_coef, mo_num*ao_num);
assert (rc == QMCKL_SUCCESS);
for (int i=0 ; i<ao_num * mo_num ; i++) {
assert (mo_coef[i] == chbrclf_mo_coef[i]);
printf("%d %e %e %e\n", i, mo_coef[i], chbrclf_mo_coef[i],
( fabs(mo_coef[i] - chbrclf_mo_coef[i])/fabs(mo_coef[i])) );
assert ( fabs(mo_coef[i] - chbrclf_mo_coef[i])/fabs(mo_coef[i]) < 1.e-12 );
}
free(mo_coef);
charge = NULL;

229
org/qmckl_utils.org
View File

@ -1,229 +0,0 @@
#+TITLE: Utility functions
#+SETUPFILE: ../tools/theme.setup
#+INCLUDE: ../tools/lib.org
* Headers :noexport:
#+begin_src elisp :noexport :results none
(org-babel-lob-ingest "../tools/lib.org")
#+end_src
#+begin_src c :comments link :tangle (eval c_test) :noweb yes
#include "qmckl.h"
#include "assert.h"
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
int main() {
qmckl_context context;
context = qmckl_context_create();
#+end_src
* Matrix operations
** ~qmckl_transpose~
Transposes a matrix: $B_{ji} = A_{ij}$
#+NAME: qmckl_transpose_args
| qmckl_context | context | in | Global state |
| int64_t | m | in | Number of rows of the input matrix |
| int64_t | n | in | Number of columns of the input matrix |
| double | A[][lda] | in | Array containing the $m \times n$ matrix $A$ |
| int64_t | lda | in | Leading dimension of array ~A~ |
| double | B[][ldb] | out | Array containing the $n \times m$ matrix $B$ |
| int64_t | ldb | in | Leading dimension of array ~B~ |
*** Requirements
- ~context~ is not ~QMCKL_NULL_CONTEXT~
- ~m > 0~
- ~n > 0~
- ~lda >= m~
- ~ldb >= n~
- ~A~ is allocated with at least $m \times n \times 8$ bytes
- ~B~ is allocated with at least $n \times m \times 8$ bytes
*** C header
#+CALL: generate_c_header(table=qmckl_transpose_args,rettyp="qmckl_exit_code",fname="qmckl_transpose")
#+RESULTS:
#+begin_src c :tangle (eval h_func) :comments org
qmckl_exit_code qmckl_transpose (
const qmckl_context context,
const int64_t m,
const int64_t n,
const double* A,
const int64_t lda,
double* const B,
const int64_t ldb );
#+end_src
*** Source
#+begin_src f90 :tangle (eval f)
integer function qmckl_transpose_f(context, m, n, A, LDA, B, LDB) &
result(info)
use qmckl
implicit none
integer(qmckl_context) , intent(in) :: context
integer*8 , intent(in) :: m, n
integer*8 , intent(in) :: lda
real*8 , intent(in) :: A(lda,*)
integer*8 , intent(in) :: ldb
real*8 , intent(out) :: B(ldb,*)
integer*8 :: i,j
info = QMCKL_SUCCESS
if (context == QMCKL_NULL_CONTEXT) then
info = QMCKL_INVALID_CONTEXT
return
endif
if (m <= 0_8) then
info = QMCKL_INVALID_ARG_2
return
endif
if (n <= 0_8) then
info = QMCKL_INVALID_ARG_3
return
endif
if (LDA < m) then
info = QMCKL_INVALID_ARG_5
return
endif
if (LDB < n) then
info = QMCKL_INVALID_ARG_7
return
endif
do j=1,m
do i=1,n
B(i,j) = A(j,i)
end do
end do
end function qmckl_transpose_f
#+end_src
*** C interface :noexport:
#+CALL: generate_c_interface(table=qmckl_transpose_args,rettyp="qmckl_exit_code",fname="qmckl_transpose")
#+RESULTS:
#+begin_src f90 :tangle (eval f) :comments org :exports none
integer(c_int32_t) function qmckl_transpose &
(context, m, n, A, lda, B, ldb) &
bind(C) result(info)
use, intrinsic :: iso_c_binding
implicit none
integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(in) , value :: m
integer (c_int64_t) , intent(in) , value :: n
real (c_double ) , intent(in) :: A(lda,*)
integer (c_int64_t) , intent(in) , value :: lda
real (c_double ) , intent(out) :: B(ldb,*)
integer (c_int64_t) , intent(in) , value :: ldb
integer(c_int32_t), external :: qmckl_transpose_f
info = qmckl_transpose_f &
(context, m, n, A, lda, B, ldb)
end function qmckl_transpose
#+end_src
#+CALL: generate_f_interface(table=qmckl_transpose_args,rettyp="qmckl_exit_code",fname="qmckl_transpose")
#+RESULTS:
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_transpose &
(context, m, n, A, lda, B, ldb) &
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 :: m
integer (c_int64_t) , intent(in) , value :: n
real (c_double ) , intent(in) :: A(lda,*)
integer (c_int64_t) , intent(in) , value :: lda
real (c_double ) , intent(out) :: B(ldb,*)
integer (c_int64_t) , intent(in) , value :: ldb
end function qmckl_transpose
end interface
#+end_src
*** Test :noexport:
#+begin_src f90 :tangle (eval f_test)
integer(qmckl_exit_code) function test_qmckl_transpose(context) bind(C)
use qmckl
implicit none
integer(qmckl_context), intent(in), value :: context
double precision, allocatable :: A(:,:), B(:,:)
integer*8 :: m, n, LDA, LDB
integer*8 :: i,j
double precision :: x
m = 5
n = 6
LDA = m+3
LDB = n+1
allocate( A(LDA,n), B(LDB,m) )
A = 0.d0
B = 0.d0
do j=1,n
do i=1,m
A(i,j) = -10.d0 + dble(i+j)
end do
end do
test_qmckl_transpose = qmckl_transpose(context, m, n, A, LDA, B, LDB)
if (test_qmckl_transpose /= QMCKL_SUCCESS) return
test_qmckl_transpose = QMCKL_FAILURE
x = 0.d0
do j=1,n
do i=1,m
x = x + (A(i,j)-B(j,i))**2
end do
end do
if (dabs(x) <= 1.d-15) then
test_qmckl_transpose = QMCKL_SUCCESS
endif
deallocate(A,B)
end function test_qmckl_transpose
#+end_src
#+begin_src c :comments link :tangle (eval c_test)
qmckl_exit_code test_qmckl_transpose(qmckl_context context);
assert(QMCKL_SUCCESS == test_qmckl_transpose(context));
#+end_src
* End of files :noexport:
#+begin_src c :comments link :tangle (eval c_test)
assert (qmckl_context_destroy(context) == QMCKL_SUCCESS);
return 0;
}
#+end_src
# -*- mode: org -*-
# vim: syntax=c

13
org/table_of_contents
View File

@ -1,17 +1,20 @@
qmckl.org
qmckl_error.org
qmckl_context.org
qmckl_error.org
qmckl_blas.org
qmckl_memory.org
qmckl_numprec.org
qmckl_distance.org
qmckl_point.org
qmckl_nucleus.org
qmckl_electron.org
qmckl_distance.org
qmckl_ao.org
qmckl_mo.org
qmckl_jastrow.org
qmckl_determinant.org
qmckl_sherman_morrison_woodbury.org
qmckl_jastrow.org
qmckl_local_energy.org
qmckl_utils.org
qmckl_blas.org
qmckl_trexio.org
qmckl_verificarlo.org
qmckl_tests.org
qmckl_verificarlo.org

1
share/doc/qmckl/html/index.html Normal file
View File

@ -0,0 +1 @@
<meta http-equiv="refresh" content="0;URL=README.html">

3620
share/qmckl/test_data/be2/ao_1e_int.txt
View File

File diff suppressed because it is too large Load Diff

102
share/qmckl/test_data/be2/basis.txt
View File

@ -1,34 +1,37 @@
rank_basis_nucleus_index 1
dims_basis_nucleus_index 0 2
rank_basis_nucleus_shell_num 1
dims_basis_nucleus_shell_num 0 2
dims_basis_nucleus_index 0 12
rank_basis_shell_ang_mom 1
dims_basis_shell_ang_mom 0 12
rank_basis_shell_prim_num 1
dims_basis_shell_prim_num 0 12
rank_basis_shell_factor 1
dims_basis_shell_factor 0 12
rank_basis_shell_prim_index 1
dims_basis_shell_prim_index 0 12
rank_basis_shell_index 1
dims_basis_shell_index 0 50
rank_basis_exponent 1
dims_basis_exponent 0 50
rank_basis_coefficient 1
dims_basis_coefficient 0 50
rank_basis_prim_factor 1
dims_basis_prim_factor 0 50
basis_num_isSet 1
basis_num 12
basis_prim_num_isSet 1
basis_prim_num 50
basis_shell_num_isSet 1
basis_shell_num 12
len_basis_type 9
basis_type
Gaussian
basis_nucleus_index
0
6
basis_nucleus_shell_num
6
6
0
0
0
0
0
1
1
1
1
1
1
basis_shell_ang_mom
0
0
@ -42,19 +45,6 @@ basis_shell_ang_mom
1
1
2
basis_shell_prim_num
9
9
1
4
1
1
9
9
1
4
1
1
basis_shell_factor
9.9999971897081508e-01
9.9999963111699008e-01
@ -68,19 +58,57 @@ basis_shell_factor
1.0000002163655846e+00
1.0000000000000000e+00
1.0000000000000002e+00
basis_shell_prim_index
basis_shell_index
0
0
0
0
0
0
0
0
0
1
1
1
1
1
1
1
1
1
2
3
3
3
3
4
5
6
6
6
6
6
6
6
6
6
7
7
7
7
7
7
7
7
7
8
9
18
19
23
24
25
34
43
44
48
49
9
9
9
10
11
basis_exponent
2.9400000000000000e+03
4.4119999999999999e+02

2
share/qmckl/test_data/be2/metadata.txt
View File

@ -4,7 +4,7 @@ metadata_code_num_isSet 0
metadata_author_num_isSet 0
len_metadata_package_version 6
metadata_package_version
1.1.0
2.0.0
len_metadata_description 0
metadata_description
metadata_code

1680
share/qmckl/test_data/be2/mo.txt
View File

File diff suppressed because it is too large Load Diff

2370
share/qmckl/test_data/be2/mo_1e_int.txt
View File

File diff suppressed because it is too large Load Diff

2
share/qmckl/test_data/be2/nucleus.txt
View File

@ -7,6 +7,8 @@ rank_nucleus_label 1
dims_nucleus_label 0 2
nucleus_num_isSet 1
nucleus_num 2
nucleus_repulsion_isSet 1
nucleus_repulsion 3.4507806369169232e+00
len_nucleus_point_group 0
nucleus_point_group
nucleus_charge

483
share/qmckl/test_data/chbrclf/basis.txt
View File

@ -1,40 +1,97 @@
rank_basis_nucleus_index 1
dims_basis_nucleus_index 0 5
rank_basis_nucleus_shell_num 1
dims_basis_nucleus_shell_num 0 5
dims_basis_nucleus_index 0 72
rank_basis_shell_ang_mom 1
dims_basis_shell_ang_mom 0 72
rank_basis_shell_prim_num 1
dims_basis_shell_prim_num 0 72
rank_basis_shell_factor 1
dims_basis_shell_factor 0 72
rank_basis_shell_prim_index 1
dims_basis_shell_prim_index 0 72
rank_basis_shell_index 1
dims_basis_shell_index 0 297
rank_basis_exponent 1
dims_basis_exponent 0 297
rank_basis_coefficient 1
dims_basis_coefficient 0 297
rank_basis_prim_factor 1
dims_basis_prim_factor 0 297
basis_num_isSet 1
basis_num 72
basis_prim_num_isSet 1
basis_prim_num 297
basis_shell_num_isSet 1
basis_shell_num 72
len_basis_type 9
basis_type
Gaussian
basis_nucleus_index
0
14
23
37
53
basis_nucleus_shell_num
14
9
14
16
19
0
0
0
0
0
0
0
0
0
0
0
0
0
1
1
1
1
1
1
1
1
1
2
2
2
2
2
2
2
2
2
2
2
2
2
2
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
4
4
4
4
4
4
4
4
4
4
4
4
4
4
4
4
4
4
4
basis_shell_ang_mom
0
0
@ -108,79 +165,6 @@ basis_shell_ang_mom
2
3
3
basis_shell_prim_num
10
10
1
1
1
5
1
1
1
1
1
1
1
1
5
1
1
1
1
1
1
1
1
10
10
1
1
1
5
1
1
1
1
1
1
1
1
15
15
15
1
1
1
9
9
1
1
1
1
1
1
1
1
20
20
20
20
1
1
1
13
13
13
1
1
1
8
1
1
1
1
1
basis_shell_factor
9.9999973955019195e-01
9.9999948164842034e-01
@ -254,13 +238,83 @@ basis_shell_factor
1.0000000000000002e+00
1.0000000000000002e+00
1.0000000000000002e+00
basis_shell_prim_index
basis_shell_index
0
0
0
0
0
0
0
0
0
0
1
1
1
1
1
1
1
1
1
1
2
3
4
5
5
5
5
5
6
7
8
9
10
11
12
13
14
14
14
14
14
15
16
17
18
19
20
21
22
23
23
23
23
23
23
23
23
23
23
24
24
24
24
24
24
24
24
24
24
25
26
27
28
28
28
28
28
29
30
@ -270,63 +324,218 @@ basis_shell_prim_index
34
35
36
37
37
37
37
37
37
37
37
37
37
37
37
37
37
37
38
38
38
38
38
38
38
38
38
38
38
38
38
38
38
39
39
39
39
39
39
39
39
39
39
39
39
39
39
39
40
41
42
43
43
43
43
43
43
43
43
43
44
44
44
44
44
44
44
44
44
45
46
47
48
49
50
51
52
53
53
53
53
53
53
53
53
53
53
53
53
53
53
53
53
53
53
53
53
54
54
54
54
54
54
54
54
54
54
54
54
54
54
54
54
54
54
54
54
55
55
55
55
55
55
55
55
55
55
55
55
55
55
55
55
55
55
55
55
56
56
56
56
56
56
56
56
56
56
56
56
56
56
56
56
56
56
56
56
57
58
59
60
60
60
60
60
60
60
60
60
60
60
60
60
61
61
61
61
61
61
61
61
61
61
61
61
61
62
62
62
62
62
62
62
62
62
62
62
62
62
63
64
65
66
66
66
66
66
66
66
66
67
68
69
70
71
72
77
78
79
80
81
82
83
84
85
100
115
130
131
132
133
142
151
152
153
154
155
156
157
158
159
179
199
219
239
240
241
242
255
268
281
282
283
284
292
293
294
295
296
basis_exponent
8.2360000000000000e+03
1.2350000000000000e+03

3
share/qmckl/test_data/chbrclf/metadata.txt
View File

@ -4,8 +4,9 @@ metadata_code_num_isSet 0
metadata_author_num_isSet 0
len_metadata_package_version 6
metadata_package_version
0.3.0
2.0.0
len_metadata_description 0
metadata_description
metadata_code
metadata_author
metadata_unsafe_isSet 0

117822
share/qmckl/test_data/chbrclf/mo.txt
View File

File diff suppressed because it is too large Load Diff

2
share/qmckl/test_data/chbrclf/nucleus.txt
View File

@ -7,6 +7,8 @@ rank_nucleus_label 1
dims_nucleus_label 0 5
nucleus_num_isSet 1
nucleus_num 5
nucleus_repulsion_isSet 1
nucleus_repulsion 3.1823098794361579e+02
len_nucleus_point_group 0
nucleus_point_group
nucleus_charge

37
share/qmckl/test_data/sm_test.h Normal file
View File

File diff suppressed because one or more lines are too long

4
tools/build_doc.sh
View File

@ -32,7 +32,7 @@ function extract_doc()
--load="${HTMLIZE}" \
--load="${CONFIG_DOC}" \
${org} \
--load="../${CONFIG_TANGLE}" \
--load="${CONFIG_TANGLE}" \
-f org-html-export-to-html \
-f org-ascii-export-to-ascii
@ -47,7 +47,7 @@ function extract_doc()
for i in $@
do
exported=${i%.org}.exported
exported=$(dirname $exported)/.$(basename $exported)
exported=${top_builddir}/src/$(basename $exported)
NOW=$(date +"%m%d%H%M.%S")
extract_doc ${i} &> $exported
rc=$?

40
tools/build_makefile.py
View File

@ -4,8 +4,17 @@
from __future__ import print_function
from glob import glob
import os
import subprocess
def main():
wd = os.getcwd()
try:
srcdir = os.environ["srcdir"]
os.chdir(srcdir)
except KeyError:
pass
C_FILES = []
C_O_FILES = []
F_FILES = []
@ -33,8 +42,8 @@ def main():
for org in glob("org/*.org"):
i = org.split('/')[-1].rsplit(".",1)[0]
tangled = "org/."+i+".tangled"
exported = "org/."+i+".exported"
tangled = "src/"+i+".tangled"
exported = "src/"+i+".exported"
c_test_x = "tests/test_"+i
c_test_o = "tests/test_"+i+".$(OBJEXT)"
f_test_o = "tests/test_"+i+"_f.$(OBJEXT)"
@ -181,7 +190,7 @@ def main():
output = ["",
"## Source files",
"",
"ORG_FILES="+" ".join(ORG_FILES),
"ORG_FILES="+" ".join([ "$(srcdir)/"+ x for x in ORG_FILES]),
"TANGLED_FILES="+" ".join(TANGLED_FILES),
"EXPORTED_FILES="+" ".join(EXPORTED_FILES),
"C_FILES="+" ".join(C_FILES),
@ -195,7 +204,7 @@ def main():
"H_PRIVATE_TYPE_FILES="+" ".join(H_PRIVATE_TYPE_FILES),
"C_TEST_FILES="+" ".join(C_TEST_FILES),
"F_TEST_FILES="+" ".join(F_TEST_FILES),
"TESTS="+" ".join(TESTS.keys()).replace("$(srcdir)/",""),
"TESTS="+" ".join(TESTS.keys()),
"HTML_FILES="+" ".join(HTML.values()),
"TEXT_FILES="+" ".join(TEXT.values()),
"" ]
@ -205,8 +214,8 @@ def main():
"",
"if QMCKL_DEVEL" ]
for f in DEPS_ORG.keys():
output += [ DEPS_ORG[f] + ": "+f,
"\t$(tangle_verbose)top_builddir=$(top_builddir) srcdir=$(srcdir) $(srcdir)/tools/tangle.sh "+f,
output += [ DEPS_ORG[f] + ": $(srcdir)/"+f,
"\t$(tangle_verbose)top_builddir=$(abs_top_builddir) srcdir=$(abs_srcdir) $(srcdir)/tools/missing bash $(srcdir)/tools/tangle.sh $(srcdir)/"+f,
"" ]
output += [ "endif",
"" ]
@ -227,21 +236,22 @@ def main():
output+= ["",
"## Test files",
"",
"$(header_tests): $(TANGLED_FILES)",
"$(test_qmckl_fo): $(test_qmckl_f)"]
output += ["",
"check_PROGRAMS = $(TESTS)" ]
for f in sorted(TESTS.keys()):
prefix = "tests_" + f.rsplit("/",1)[-1]
output += [ prefix + "_SOURCES = " + \
" ".join(TESTS[f]).replace("$(srcdir)",""),
" ".join(TESTS[f]) + " $(header_tests)",
prefix + "_LDADD = src/libqmckl.la",
prefix + "_LDFLAGS = -no-install",
"" ]
tmp = "EXTRA_DIST += "
for dir in glob("share/qmckl/test_data/*"):
for f in glob("%s/*"%(dir)):
tmp += " \\\n "+f
r = subprocess.check_output("git ls-tree --name-only -r HEAD".split())
for line in r.splitlines():
if b"share/qmckl/test_data/" in line:
tmp += " \\\n " + line.decode('utf8')
tmp += "\n"
output += tmp.split("\n")
@ -252,16 +262,16 @@ def main():
for f in sorted(ORG_FILES):
output += [ HTML[f] + ": " + DEPS_DOC[f],
TEXT[f] + ": " + DEPS_DOC[f],
TEXT[f] + ": " + DEPS_DOC[f],
"" ]
for f in sorted(DEPS_DOC.keys()):
output += [ DEPS_DOC[f] + ": " + f + " $(htmlize_el)",
"\t$(export_verbose)top_builddir=$(top_builddir) srcdir=$(srcdir) $(srcdir)/tools/build_doc.sh "+f,
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,
"" ]
output += ["endif"]
f = open("generated.mk","w")
f = open(srcdir+"/generated.mk","w")
f.write("\n".join(output))

2
tools/build_qmckl_h.sh
View File

@ -99,7 +99,7 @@ EOF
for i in ${HEADERS}
do
header=${srcdir}/src/$i
header=${top_builddir}/src/$i
if [ -f $header ] ; then
echo "/* $header */" >> ${qmckl_h}
cat $header >> ${qmckl_h}

32
tools/config_tangle.el
View File

@ -33,19 +33,21 @@
; The following is required to compute the file names
(setq pwd (file-name-directory buffer-file-name))
(setq wd (concat pwd "/../src/"))
(setq td (concat pwd "/../tests/"))
(setq name (file-name-nondirectory (substring buffer-file-name 0 -4)))
(setq f (concat wd name "_f.f90"))
(setq fh_func (concat wd name "_fh_func.f90"))
(setq fh_type (concat wd name "_fh_type.f90"))
(setq c (concat wd name ".c"))
(setq h_func (concat wd name "_func.h"))
(setq h_type (concat wd name "_type.h"))
(setq h_private_type (concat wd name "_private_type.h"))
(setq h_private_func (concat wd name "_private_func.h"))
(setq c_test (concat td "test_" name ".c"))
(setq f_test (concat td "test_" name "_f.f90"))
(org-babel-lob-ingest "../tools/lib.org")
(setq top_builddir (or (getenv "top_builddir") "."))
(setq srcdir (or (getenv "srcdir") "."))
(setq src (concat top_builddir "/src/"))
(setq tests (concat top_builddir "/tests/"))
(setq name (file-name-nondirectory (substring buffer-file-name 0 -4)))
(setq f (concat src name "_f.f90"))
(setq fh_func (concat src name "_fh_func.f90"))
(setq fh_type (concat src name "_fh_type.f90"))
(setq c (concat src name ".c"))
(setq h_func (concat src name "_func.h"))
(setq h_type (concat src name "_type.h"))
(setq h_private_type (concat src name "_private_type.h"))
(setq h_private_func (concat src name "_private_func.h"))
(setq c_test (concat tests "test_" name ".c"))
(setq f_test (concat tests "test_" name "_f.f90"))
(org-babel-lob-ingest (concat srcdir "/tools/lib.org"))

2
tools/install_htmlize.sh
View File

@ -2,7 +2,7 @@
#
# Installs the htmlize Emacs plugin
./tools/missing git clone "https://github.com/hniksic/emacs-htmlize"
${abs_srcdir}/tools/missing git clone "https://github.com/TREX-CoE/emacs-htmlize"
mv emacs-htmlize/htmlize.el $1
rm -rf emacs-htmlize

53
tools/lib.org
View File

@ -19,17 +19,19 @@
** Table of function arguments
#+NAME: test
| ~qmckl_context~ | ~context~ | in | Global state |
| ~char~ | ~transa~ | in | Array ~A~ is ~'N'~: Normal, ~'T'~: Transposed |
| ~char~ | ~transb~ | in | Array ~B~ is ~'N'~: Normal, ~'T'~: Transposed |
| ~int64_t~ | ~m~ | in | Number of points in the first set |
| ~int64_t~ | ~n~ | in | Number of points in the second set |
| ~double~ | ~A[][lda]~ | in | Array containing the $m \times 3$ matrix $A$ |
| ~int64_t~ | ~lda~ | in | Leading dimension of array ~A~ |
| ~double~ | ~B[][ldb]~ | in | Array containing the $n \times 3$ matrix $B$ |
| ~int64_t~ | ~ldb~ | in | Leading dimension of array ~B~ |
| ~double~ | ~C[n][ldc]~ | out | Array containing the $m \times n$ matrix $C$ |
| ~int64_t~ | ~ldc~ | in | Leading dimension of array ~C~ |
| Variable | Type | In/Out | Description |
|-----------+------------------+--------+-----------------------------------------------|
| ~context~ | ~qmckl_context~ | in | Global state |
| ~transa~ | ~char~ | in | Array ~A~ is ~'N'~: Normal, ~'T'~: Transposed |
| ~transb~ | ~char~ | in | Array ~B~ is ~'N'~: Normal, ~'T'~: Transposed |
| ~m~ | ~int64_t~ | in | Number of points in the first set |
| ~n~ | ~int64_t~ | in | Number of points in the second set |
| ~A~ | ~double[][lda]~ | in | Array containing the $m \times 3$ matrix $A$ |
| ~lda~ | ~int64_t~ | in | Leading dimension of array ~A~ |
| ~B~ | ~double[][ldb]~ | in | Array containing the $n \times 3$ matrix $B$ |
| ~ldb~ | ~int64_t~ | in | Leading dimension of array ~B~ |
| ~C~ | ~double[n][ldc]~ | out | Array containing the $m \times n$ matrix $C$ |
| ~ldc~ | ~int64_t~ | in | Leading dimension of array ~C~ |
*** Fortran-C type conversions
@ -57,8 +59,6 @@ ctypeid_d = { '' : ''
, 'qmckl_exit_code' : 'integer(c_int32_t)'
, 'integer' : 'integer(c_int32_t)'
, 'integer*8' : 'integer(c_int64_t)'
, 'integer' : 'integer(c_uint32_t)'
, 'integer*8' : 'integer(c_uint64_t)'
, 'real' : 'real(c_float)'
, 'real*8' : 'real(c_double)'
, 'character' : 'character(c_char)'
@ -74,9 +74,9 @@ def parse_table(table):
result = []
for line in [ [x.replace('~','') for x in y] for y in table]:
d = { "c_type" : line[0],
d = { "name" : line[0],
"c_type" : line[1],
"inout" : line[2].lower(),
"name" : line[1],
"comment" : line[3] }
# Handle inout
@ -88,12 +88,12 @@ def parse_table(table):
d["inout"] == "inout"
# Find dimensions (replace [] by [*] to get * in Fortran dimensions)
dims = d["name"].replace("[]","[*]").split('[')
dims = d["c_type"].replace("[]","[*]").split('[')
d["rank"] = len(dims) - 1
if d["rank"] == 0:
d["dims"] = []
else:
d["name"] = d["name"].split('[')[0].strip()
d["c_type"] = d["c_type"].split('[')[0].strip()
d["dims"] = [ x.replace(']','').strip() for x in dims[1:] ]
result.append(d)
@ -104,7 +104,7 @@ def parse_table(table):
*** Generates a C header
#+NAME: generate_c_header
#+BEGIN_SRC python :var table=test :var rettyp=[] :var fname=[] :results drawer :noweb yes :wrap "src c :tangle (eval h_func) :comments org"
#+BEGIN_SRC python :var table=test :var rettyp="qmckl_exit_code" :var fname=[] :results drawer :noweb yes :wrap "src c :tangle (eval h_func) :comments org"
<<parse_table>>
results = []
@ -133,6 +133,23 @@ return template
#+END_SRC
#+RESULTS: generate_c_header
#+begin_src c :tangle (eval h_func) :comments org
qmckl_exit_code [] (
const qmckl_context context,
const char transa,
const char transb,
const int64_t m,
const int64_t n,
const double* A,
const int64_t lda,
const double* B,
const int64_t ldb,
double* const C,
const int64_t ldc );
#+end_src
*** Generates a C interface to the Fortran function
#+NAME: generate_c_interface

8
tools/tangle.sh
View File

@ -14,7 +14,7 @@ if [[ -z ${srcdir} ]] ; then
fi
if [[ -z ${top_builddir} ]] ; then
echo "Error: srcdir environment variable is not defined"
echo "Error: top_builddir environment variable is not defined"
exit 1
fi
@ -31,14 +31,14 @@ function tangle()
fi
${srcdir}/tools/missing \
emacs --batch ${org_file} \
--load=${PWD}/tools/config_tangle.el \
--load=${srcdir}/tools/config_tangle.el \
-f org-babel-tangle
}
for i in $@
do
tangled=${i%.org}.tangled
tangled=$(dirname $tangled)/.$(basename $tangled)
tangled=${top_builddir}/src/$(basename $tangled)
NOW=$(date +"%m%d%H%M.%S")
tangle ${i} &> $tangled
rc=$?
@ -47,6 +47,6 @@ do
# Fail if tangling failed
if [[ $rc -ne 0 ]] ; then
cat $tangled
exit rc
exit $rc
fi
done