TREX/qmckl
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2025-01-03 10:06:09 +01:00
Code Issues Projects Releases Wiki Activity

Fix examples in doc

Browse Source
This commit is contained in:
Anthony Scemama 2022-06-15 19:17:05 +02:00
commit 28fe475425
34 changed files with 56071 additions and 1644 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -2,9 +2,7 @@ name: test-build
on:
push:
branches: [ master ]
pull_request:
branches: [ master ]
jobs:
x86_ubuntu:
@ -44,6 +42,7 @@ jobs:
cd _build
../configure --enable-silent-rules --enable-debug
make -j 4
sudo make install
- name: Run test
run: make -j 4 check
@ -60,60 +59,77 @@ jobs:
run: make distcheck
working-directory: _build
x86_macos:
runs-on: macos-latest
name: x86 MacOS latest
steps:
- uses: actions/checkout@v2
- name: install dependencies
run: brew install emacs hdf5 automake pkg-config
- name: Symlink gfortran (macOS)
if: runner.os == 'macOS'
- name: Setup the virtual environment
run: |
# make sure gfortran is available
# https://github.com/actions/virtual-environments/issues/2524
# https://github.com/cbg-ethz/dce/blob/master/.github/workflows/pkgdown.yaml
sudo ln -s /usr/local/bin/gfortran-10 /usr/local/bin/gfortran
sudo mkdir /usr/local/gfortran
sudo ln -s /usr/local/Cellar/gcc@10/*/lib/gcc/10 /usr/local/gfortran/lib
gfortran --version
python3 -m venv --clear pyqmckl
source pyqmckl/bin/activate
- name: Install the latest TREXIO from the GitHub clone
run: |
git clone https://github.com/TREX-CoE/trexio.git
cd trexio
./autogen.sh
./configure --prefix=${PWD}/_install --enable-silent-rules
make -j 4
make install
- name: Install the Python requirements
run: pip install -r requirements.txt
working-directory: python
- name: Test TREXIO
run: make -j 4 check
working-directory: trexio
- name: Install the Python API
run: make python-install
working-directory: _build
- name: Archive TREXIO test log file
if: failure()
uses: actions/upload-artifact@v2
with:
name: test-report-trexio-macos
path: trexio/test-suite.log
- name: Test the Python API
run: make python-test
working-directory: _build
- name: Build QMCkl
run: |
export PKG_CONFIG_PATH=${PWD}/trexio/_install/lib/pkgconfig:$PKG_CONFIG_PATH
./autogen.sh
./configure CC=gcc-10 FC=gfortran-10 --enable-silent-rules --enable-debug
make -j 4
- name: Run test
run: make -j 4 check
- name: Archive test log file
if: failure()
uses: actions/upload-artifact@v2
with:
name: test-report-macos
path: test-suite.log
# x86_macos:
#
# runs-on: macos-latest
# name: x86 MacOS latest
#
# steps:
# - uses: actions/checkout@v2
# - name: install dependencies
# run: brew install emacs hdf5 automake pkg-config
#
# - name: Symlink gfortran (macOS)
# if: runner.os == 'macOS'
# run: |
# # make sure gfortran is available
# # https://github.com/actions/virtual-environments/issues/2524
# # https://github.com/cbg-ethz/dce/blob/master/.github/workflows/pkgdown.yaml
# sudo ln -s /usr/local/bin/gfortran-10 /usr/local/bin/gfortran
# sudo mkdir /usr/local/gfortran
# sudo ln -s /usr/local/Cellar/gcc@10/*/lib/gcc/10 /usr/local/gfortran/lib
# gfortran --version
#
# - name: Install the latest TREXIO from the GitHub clone
# run: |
# git clone https://github.com/TREX-CoE/trexio.git
# cd trexio
# ./autogen.sh
# ./configure --prefix=${PWD}/_install --enable-silent-rules
# make -j 4
# make install
#
# - name: Test TREXIO
# run: make -j 4 check
# working-directory: trexio
#
# - name: Archive TREXIO test log file
# if: failure()
# uses: actions/upload-artifact@v2
# with:
# name: test-report-trexio-macos
# path: trexio/test-suite.log
#
# - name: Build QMCkl
# run: |
# export PKG_CONFIG_PATH=${PWD}/trexio/_install/lib/pkgconfig:$PKG_CONFIG_PATH
# ./autogen.sh
# ./configure CC=gcc-10 FC=gfortran-10 --enable-silent-rules
# make -j 4
#
# - name: Run test
# run: make -j 4 check
#
# - name: Archive test log file
# if: failure()
# uses: actions/upload-artifact@v2
# with:
# name: test-report-macos
# path: test-suite.log

41
Makefile.am
View File

@ -65,6 +65,7 @@ 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)
src_libqmckl_la_LDFLAGS = $(LDFLAGS)
CLEANFILES+=$(test_qmckl_fo) $(src_qmckl_fo) $(test_qmckl_o) $(src_qmckl_o) $(FH_TYPE_FILES) $(FH_FUNC_FILES)
@ -169,6 +170,44 @@ cppcheck.out: $(qmckl_h)
--language=c --std=c99 -rp --platform=unix64 \
-I$(srcdir)/include -I$(top_builddir)/include *.c *.h 2>../$@
.PHONY: cppcheck
setup_py = $(srcdir)/python/setup.py
process_header_py = $(srcdir)/python/src/process_header.py
test_py = $(srcdir)/python/test/test_api.py
qmckl_i = $(srcdir)/python/src/qmckl.i
numpy_i = $(srcdir)/python/src/numpy.i
qmckl_wrap_c = python/src/qmckl_wrap.c
qmckl_include_i = python/src/qmckl_include.i
qmckl_py = python/qmckl/qmckl.py
$(qmckl_include_i): $(qmckl_h) $(process_header_py)
$(MKDIR_P) python/src
python $(process_header_py) $(qmckl_h)
mv qmckl_include.i $(qmckl_include_i)
$(qmckl_py): $(qmckl_i) $(qmckl_include_i)
swig -Iinclude -Ipython/src -python -py3 -builtin -o $(qmckl_wrap_c) $(qmckl_i)
$(qmckl_wrap_c): $(qmckl_py)
python-install: $(qmckl_h) $(qmckl_i) $(setup_py) $(qmckl_py) $(qmckl_wrap_c)
$(MKDIR_P) python/src
cd python ; \
[[ ! -f pyproject.toml ]] && \
cp $(abs_srcdir)/python/{pyproject.toml,requirements.txt,README.md,setup.py} . ; \
cp src/qmckl.py . ; \
export QMCKL_INCLUDEDIR="$(prefix)/include" ; \
export QMCKL_LIBDIR="$(prefix)/lib" ; \
pip install .
python-test: $(test_py)
cd $(abs_srcdir)/python/test/ && \
python test_api.py
CLEANFILES += $(qmckl_wrap_c) \
$(qmckl_include_i) \
$(qmckl_py)
.PHONY: cppcheck python-test python-install
endif

2
autogen.sh
View File

@ -1,5 +1,5 @@
#!/bin/bash
export srcdir="."
python ${srcdir}/tools/build_makefile.py
python3 ${srcdir}/tools/build_makefile.py
autoreconf -i -Wall --no-recursive

173
configure.ac
View File

@ -93,6 +93,7 @@ AC_PROG_F77
m4_version_prereq([2.70],[], [AC_PROG_CC_C99])
AS_IF([test "$ac_cv_prog_cc_c99" = "no"], [AC_MSG_ERROR([The compiler does not support C99])])
AC_PROG_CC_C_O
AM_PROG_CC_C_O
AC_PROG_FC
AC_PROG_FC_C_O
AC_FC_PP_DEFINE
@ -137,10 +138,10 @@ case "$with_chameleon" in
[PKG_CFLAGS="$PKG_CFLAGS $LIBCHAMELEON_CFLAGS"
PKG_LIBS="$PKG_LIBS $LIBCHAMELEON_LIBS"]
,[
## something went wrong.
## try to find the package without pkg-config
## check that the library is actually new enough.
## by testing for a 1.0.0+ function which we use
AC_CHECK_LIB(chameleon,CHAMELEON_finalize,[LIBCHAMELEON_LIBS="-lchameleon"])
@ -212,6 +213,9 @@ esac
case $CC in
*gcc*)
CFLAGS="$CFLAGS -fPIC"
;;
*nvc*)
CFLAGS="$CFLAGS -fPIC"
;;
@ -224,6 +228,131 @@ AS_IF([test "$HAVE_HPC" = "yes"], [
AC_DEFINE([HAVE_HPC], [1], [If defined, activate HPC routines])
])
# Enable Verificarlo tests
AC_ARG_ENABLE([vfc_ci],
[ --enable-vfc_ci Build the library with vfc_ci support],
[case "${enableval}" in
yes) vfc_ci=true && FCFLAGS="-D VFC_CI $FCFLAGS" && CFLAGS="-D VFC_CI $CFLAGS";;
no) vfc_ci=false ;;
*) AC_MSG_ERROR([bad value ${enableval} for --enable_vfc_ci]) ;;
esac],[vfc_ci=false])
AM_CONDITIONAL([VFC_CI], [test x$vfc_ci = xtrue])
if test "$FC" = "verificarlo-f"; then
AC_MSG_NOTICE(verificarlo-f detected)
# Arguments order is important here
FCFLAGS="-Mpreprocess $FCFLAGS"
fi
## Enable GPU offloading
# GPU offloading
AC_ARG_ENABLE(gpu, [AS_HELP_STRING([--enable-gpu],[openmp|openacc : Use GPU-offloaded functions])], enable_gpu=$enableval, enable_gpu=no)
AS_IF([test "$enable_gpu" = "yes"], [enable_gpu="openmp"])
# OpenMP offloading
HAVE_OPENMP_OFFLOAD="no"
AS_IF([test "$enable_gpu" = "openmp"], [
AC_DEFINE([HAVE_OPENMP_OFFLOAD], [1], [If defined, activate OpenMP-offloaded routines])
HAVE_OPENMP_OFFLOAD="yes"
case $CC in
*gcc*)
CFLAGS="$CFLAGS -fopenmp"
;;
*nvc*)
CFLAGS="$CFLAGS -mp=gpu"
;;
esac
case $FC in
*gfortran*)
FCFLAGS="$FCFLAGS -fopenmp"
;;
*nvfortran*)
FCFLAGS="$FCFLAGS -mp=gpu"
;;
esac]
)
# OpenMP offloading
HAVE_OPENACC_OFFLOAD="no"
AS_IF([test "$enable_gpu" = "openacc"], [
AC_DEFINE([HAVE_OPENACC_OFFLOAD], [1], [If defined, activate OpenACC-offloaded routines])
HAVE_OPENACC_OFFLOAD="yes"
case $CC in
*gcc*)
CFLAGS="$CFLAGS -fopenacc"
;;
*nvc*)
CFLAGS="$CFLAGS -acc=gpu"
;;
esac
case $FC in
*gfortran*)
FCFLAGS="$FCFLAGS -fopenacc"
;;
*nvfortran*)
FCFLAGS="$FCFLAGS -acc=gpu"
;;
esac
])
# cuBLAS offloading
AC_ARG_WITH(cublas, [AS_HELP_STRING([--with-cublas],[Use cuBLAS-offloaded functions])], HAVE_CUBLAS_OFFLOAD=$withval, HAVE_CUBLAS_OFFLOAD=no)
AS_IF([test "$HAVE_CUBLAS_OFFLOAD" = "yes"], [
AC_DEFINE([HAVE_CUBLAS_OFFLOAD], [1], [If defined, activate cuBLAS-offloaded routines])
HAVE_OPENACC_OFFLOAD="yes"
case $CC in
*gcc*)
CFLAGS="$CFLAGS -fopenmp"
LDFLAGS="-lcublas"
;;
*nvc*)
CFLAGS="$CFLAGS -mp=gpu -cudalib=cublas"
;;
esac
case $FC in
*gfortran*)
FCFLAGS="$FCFLAGS -fopenmp"
;;
*nvfortran*)
FCFLAGS="$FCFLAGS -mp=gpu -cudalib=cublas"
;;
esac
])
AC_ARG_ENABLE(malloc-trace, [AS_HELP_STRING([--enable-malloc-trace],[use debug malloc/free])], ok=$enableval, ok=no)
if test "$ok" = "yes"; then
AC_DEFINE(MALLOC_TRACE,"malloc_trace.dat",[Define to use debugging malloc/free])
ARGS="${ARGS} malloc-trace"
fi
AC_ARG_ENABLE(prof, [AS_HELP_STRING([--enable-prof],[compile for profiling])], ok=$enableval, ok=no)
if test "$ok" = "yes"; then
CFLAGS="${CFLAGS} -pg"
AC_DEFINE(ENABLE_PROF,1,[Define when using the profiler tool])
ARGS="${ARGS} prof"
fi
AC_ARG_WITH(efence, [AS_HELP_STRING([--with-efence],[use ElectricFence library])], ok=$withval, ok=no)
if test "$ok" = "yes"; then
AC_CHECK_LIB([efence], [malloc])
ARGS="${ARGS} efence"
fi
##
AC_ARG_ENABLE(debug, [AS_HELP_STRING([--enable-debug],[compile for debugging])], ok=$enableval, ok=no)
if test "$ok" = "yes"; then
if test "$GCC" = "yes"; then
@ -247,26 +376,6 @@ if test "$ok" = "yes"; then
ARGS="${ARGS} debug"
fi
AC_ARG_ENABLE(malloc-trace, [AS_HELP_STRING([--enable-malloc-trace],[use debug malloc/free])], ok=$enableval, ok=no)
if test "$ok" = "yes"; then
AC_DEFINE(MALLOC_TRACE,"malloc_trace.dat",[Define to use debugging malloc/free])
ARGS="${ARGS} malloc-trace"
fi
AC_ARG_ENABLE(prof, [AS_HELP_STRING([--enable-prof],[compile for profiling])], ok=$enableval, ok=no)
if test "$ok" = "yes"; then
CFLAGS="${CFLAGS} -pg"
AC_DEFINE(ENABLE_PROF,1,[Define when using the profiler tool])
ARGS="${ARGS} prof"
fi
AC_ARG_WITH(efence, [AS_HELP_STRING([--with-efence],[use ElectricFence library])], ok=$withval, ok=no)
if test "$ok" = "yes"; then
AC_CHECK_LIB(efence, malloc)
ARGS="${ARGS} efence"
fi
# Checks for header files.
## qmckl
@ -317,23 +426,10 @@ if test "x${QMCKL_DEVEL}" != "x"; then
HAS_CPPCHECK=1
fi
AX_PKG_SWIG(4.0.0, [], AC_MSG_WARN([SWIG is required to build Python API.]) )
fi
# Enable Verificarlo tests
AC_ARG_ENABLE([vfc_ci],
[ --enable-vfc_ci Build the library with vfc_ci support],
[case "${enableval}" in
yes) vfc_ci=true && FCFLAGS="-D VFC_CI $FCFLAGS" && CFLAGS="-D VFC_CI $CFLAGS";;
no) vfc_ci=false ;;
*) AC_MSG_ERROR([bad value ${enableval} for --enable_vfc_ci]) ;;
esac],[vfc_ci=false])
AM_CONDITIONAL([VFC_CI], [test x$vfc_ci = xtrue])
if test "$FC" = "verificarlo-f"; then
AC_MSG_NOTICE(verificarlo-f detected)
# Arguments order is important here
FCFLAGS="-Mpreprocess $FCFLAGS"
fi
#PKG-CONFIG
#mkl-dynamic-lp64-seq
@ -369,6 +465,9 @@ LDFLAGS:........: ${LDFLAGS}
LIBS............: ${LIBS}
USE CHAMELEON...: ${with_chameleon}
HPC version.....: ${HAVE_HPC}
OpenMP offload..: ${HAVE_OPENMP_OFFLOAD}
OpenACC offload.: ${HAVE_OPENACC_OFFLOAD}
cuBLAS offload..: ${HAVE_CUBLAS_OFFLOAD}
Package features:
${ARGS}

139
m4/ax_pkg_swig.m4 Normal file
View File

@ -0,0 +1,139 @@
# ===========================================================================
# https://www.gnu.org/software/autoconf-archive/ax_pkg_swig.html
# ===========================================================================
#
# SYNOPSIS
#
# AX_PKG_SWIG([major.minor.micro], [action-if-found], [action-if-not-found])
#
# DESCRIPTION
#
# This macro searches for a SWIG installation on your system. If found,
# then SWIG is AC_SUBST'd; if not found, then $SWIG is empty. If SWIG is
# found, then SWIG_LIB is set to the SWIG library path, and AC_SUBST'd.
#
# You can use the optional first argument to check if the version of the
# available SWIG is greater than or equal to the value of the argument. It
# should have the format: N[.N[.N]] (N is a number between 0 and 999. Only
# the first N is mandatory.) If the version argument is given (e.g.
# 1.3.17), AX_PKG_SWIG checks that the swig package is this version number
# or higher.
#
# As usual, action-if-found is executed if SWIG is found, otherwise
# action-if-not-found is executed.
#
# In configure.in, use as:
#
# AX_PKG_SWIG(1.3.17, [], [ AC_MSG_ERROR([SWIG is required to build..]) ])
# AX_SWIG_ENABLE_CXX
# AX_SWIG_MULTI_MODULE_SUPPORT
# AX_SWIG_PYTHON
#
# LICENSE
#
# Copyright (c) 2008 Sebastian Huber <sebastian-huber@web.de>
# Copyright (c) 2008 Alan W. Irwin
# Copyright (c) 2008 Rafael Laboissiere <rafael@laboissiere.net>
# Copyright (c) 2008 Andrew Collier
# Copyright (c) 2011 Murray Cumming <murrayc@openismus.com>
# Copyright (c) 2021 Vincent Danjean <Vincent.Danjean@ens-lyon.org>
#
# This program is free software; you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by the
# Free Software Foundation; either version 2 of the License, or (at your
# option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
# Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program. If not, see <https://www.gnu.org/licenses/>.
#
# As a special exception, the respective Autoconf Macro's copyright owner
# gives unlimited permission to copy, distribute and modify the configure
# scripts that are the output of Autoconf when processing the Macro. You
# need not follow the terms of the GNU General Public License when using
# or distributing such scripts, even though portions of the text of the
# Macro appear in them. The GNU General Public License (GPL) does govern
# all other use of the material that constitutes the Autoconf Macro.
#
# This special exception to the GPL applies to versions of the Autoconf
# Macro released by the Autoconf Archive. When you make and distribute a
# modified version of the Autoconf Macro, you may extend this special
# exception to the GPL to apply to your modified version as well.
#serial 14
AC_DEFUN([AX_PKG_SWIG],[
# Find path to the "swig" executable.
AC_PATH_PROGS([SWIG],[swig swig3.0 swig2.0])
if test -z "$SWIG" ; then
m4_ifval([$3],[$3],[:])
elif test -z "$1" ; then
m4_ifval([$2],[$2],[:])
else
AC_MSG_CHECKING([SWIG version])
[swig_version=`$SWIG -version 2>&1 | grep 'SWIG Version' | sed 's/.*\([0-9][0-9]*\.[0-9][0-9]*\.[0-9][0-9]*\).*/\1/g'`]
AC_MSG_RESULT([$swig_version])
if test -n "$swig_version" ; then
# Calculate the required version number components
[required=$1]
[required_major=`echo $required | sed 's/[^0-9].*//'`]
if test -z "$required_major" ; then
[required_major=0]
fi
[required=`echo $required. | sed 's/[0-9]*[^0-9]//'`]
[required_minor=`echo $required | sed 's/[^0-9].*//'`]
if test -z "$required_minor" ; then
[required_minor=0]
fi
[required=`echo $required. | sed 's/[0-9]*[^0-9]//'`]
[required_patch=`echo $required | sed 's/[^0-9].*//'`]
if test -z "$required_patch" ; then
[required_patch=0]
fi
# Calculate the available version number components
[available=$swig_version]
[available_major=`echo $available | sed 's/[^0-9].*//'`]
if test -z "$available_major" ; then
[available_major=0]
fi
[available=`echo $available | sed 's/[0-9]*[^0-9]//'`]
[available_minor=`echo $available | sed 's/[^0-9].*//'`]
if test -z "$available_minor" ; then
[available_minor=0]
fi
[available=`echo $available | sed 's/[0-9]*[^0-9]//'`]
[available_patch=`echo $available | sed 's/[^0-9].*//'`]
if test -z "$available_patch" ; then
[available_patch=0]
fi
# Convert the version tuple into a single number for easier comparison.
# Using base 100 should be safe since SWIG internally uses BCD values
# to encode its version number.
required_swig_vernum=`expr $required_major \* 10000 \
\+ $required_minor \* 100 \+ $required_patch`
available_swig_vernum=`expr $available_major \* 10000 \
\+ $available_minor \* 100 \+ $available_patch`
if test $available_swig_vernum -lt $required_swig_vernum; then
AC_MSG_WARN([SWIG version >= $1 is required. You have $swig_version.])
SWIG=''
m4_ifval([$3],[$3],[])
else
AC_MSG_CHECKING([for SWIG library])
SWIG_LIB=`$SWIG -swiglib`
AC_MSG_RESULT([$SWIG_LIB])
m4_ifval([$2],[$2],[])
fi
else
AC_MSG_WARN([cannot determine SWIG version])
SWIG=''
m4_ifval([$3],[$3],[])
fi
fi
AC_SUBST([SWIG_LIB])
])

199
org/examples.org
View File

@ -1,199 +0,0 @@
#+TITLE: Code examples
#+SETUPFILE: ../tools/theme.setup
#+INCLUDE: ../tools/lib.org
In this section, we present examples of usage of QMCkl.
For simplicity, we assume that the wave function parameters are stores
in a [[https://github.com/TREX-CoE/trexio][TREXIO]] file.
* Checking errors
All QMCkl functions return an error code. A convenient way to handle
errors is to write an error-checking function that displays the
error in text format and exits the program.
#+NAME: qmckl_check_error
#+begin_src f90
subroutine qmckl_check_error(rc, message)
use qmckl
implicit none
integer(qmckl_exit_code), intent(in) :: rc
character(len=*) , intent(in) :: message
character(len=128) :: str_buffer
if (rc /= QMCKL_SUCCESS) then
print *, message
call qmckl_string_of_error(rc, str_buffer)
print *, str_buffer
call exit(rc)
end if
end subroutine qmckl_check_error
#+end_src
* Computing an atomic orbital on a grid
:PROPERTIES:
:header-args: :tangle ao_grid.f90
:END:
The following program, in Fortran, computes the values of an atomic
orbital on a regular 3-dimensional grid. The 100^3 grid points are
automatically defined, such that the molecule fits in a box with 5
atomic units in the borders.
This program uses the ~qmckl_check_error~ function defined above.
To use this program, run
#+begin_src bash :tangle no
$ ao_grid <trexio_file> <AO_id> <point_num>
#+end_src
#+begin_src f90 :noweb yes
<<qmckl_check_error>>
program ao_grid
use qmckl
implicit none
integer(qmckl_context) :: qmckl_ctx ! QMCkl context
integer(qmckl_exit_code) :: rc ! Exit code of QMCkl functions
character(len=128) :: trexio_filename
character(len=128) :: str_buffer
integer :: ao_id
integer :: point_num_x
integer(c_int64_t) :: nucl_num
double precision, allocatable :: nucl_coord(:,:)
integer(c_int64_t) :: point_num
integer(c_int64_t) :: ao_num
integer(c_int64_t) :: ipoint, i, j, k
double precision :: x, y, z, dr(3)
double precision :: rmin(3), rmax(3)
double precision, allocatable :: points(:,:)
double precision, allocatable :: ao_vgl(:,:,:)
#+end_src
Start by fetching the command-line arguments:
#+begin_src f90
if (iargc() /= 3) then
print *, 'Syntax: ao_grid <trexio_file> <AO_id> <point_num>'
call exit(-1)
end if
call getarg(1, trexio_filename)
call getarg(2, str_buffer)
read(str_buffer, *) ao_id
call getarg(3, str_buffer)
read(str_buffer, *) point_num_x
if (point_num_x < 0 .or. point_num_x > 300) then
print *, 'Error: 0 < point_num < 300'
call exit(-1)
end if
#+end_src
Create the QMCkl context and initialize it with the wave function
present in the TREXIO file:
#+begin_src f90
qmckl_ctx = qmckl_context_create()
rc = qmckl_trexio_read(qmckl_ctx, trexio_filename, 1_8*len(trim(trexio_filename)))
call qmckl_check_error(rc, 'Read TREXIO')
#+end_src
We need to check that ~ao_id~ is in the range, so we get the total
number of AOs from QMCkl:
#+begin_src f90
rc = qmckl_get_ao_basis_ao_num(qmckl_ctx, ao_num)
call qmckl_check_error(rc, 'Getting ao_num')
if (ao_id < 0 .or. ao_id > ao_num) then
print *, 'Error: 0 < ao_id < ', ao_num
call exit(-1)
end if
#+end_src
Now we will compute the limits of the box in which the molecule fits.
For that, we first need to ask QMCkl the coordinates of nuclei.
#+begin_src f90
rc = qmckl_get_nucleus_num(qmckl_ctx, nucl_num)
call qmckl_check_error(rc, 'Get nucleus num')
allocate( nucl_coord(3, nucl_num) )
rc = qmckl_get_nucleus_coord(qmckl_ctx, 'N', nucl_coord, 3_8*nucl_num)
call qmckl_check_error(rc, 'Get nucleus coord')
#+end_src
We now compute the coordinates of opposite points of the box, and
the distance between points along the 3 directions:
#+begin_src f90
rmin(1) = minval( nucl_coord(1,:) ) - 5.d0
rmin(2) = minval( nucl_coord(2,:) ) - 5.d0
rmin(3) = minval( nucl_coord(3,:) ) - 5.d0
rmax(1) = maxval( nucl_coord(1,:) ) + 5.d0
rmax(2) = maxval( nucl_coord(2,:) ) + 5.d0
rmax(3) = maxval( nucl_coord(3,:) ) + 5.d0
dr(1:3) = (rmax(1:3) - rmin(1:3)) / dble(point_num_x-1)
#+end_src
We now produce the list of point coordinates where the AO will be
evaluated:
#+begin_src f90
point_num = point_num_x**3
allocate( points(point_num, 3) )
ipoint=0
z = rmin(3)
do k=1,point_num_x
y = rmin(2)
do j=1,point_num_x
x = rmin(1)
do i=1,point_num_x
ipoint = ipoint+1
points(ipoint,1) = x
points(ipoint,2) = y
points(ipoint,3) = z
x = x + dr(1)
end do
y = y + dr(2)
end do
z = z + dr(3)
end do
#+end_src
We give the points to QMCkl:
#+begin_src f90
rc = qmckl_set_point(qmckl_ctx, 'T', points, point_num)
call qmckl_check_error(rc, 'Setting points')
#+end_src
We allocate the space required to retrieve the values, gradients and
Laplacian of all AOs, and ask to retrieve the values of the
AOs computed at the point positions.
#+begin_src f90
allocate( ao_vgl(ao_num, 5, point_num) )
rc = qmckl_get_ao_basis_ao_vgl(qmckl_ctx, ao_vgl, ao_num*5_8*point_num)
call qmckl_check_error(rc, 'Setting points')
#+end_src
We finally print the value of the AO:
#+begin_src f90
do ipoint=1, point_num
print '(3(F16.10,X),E20.10)', points(ipoint, 1:3), ao_vgl(ao_id,1,ipoint)
end do
#+end_src
#+begin_src f90
deallocate( nucl_coord, points, ao_vgl )
end program ao_grid
#+end_src

1198
org/qmckl_ao.org
View File

File diff suppressed because it is too large Load Diff

6
org/qmckl_blas.org
View File

@ -84,8 +84,8 @@ are not intended to be passed to external codes.
#+begin_src c :comments org :tangle (eval h_private_type) :exports none
typedef struct qmckl_vector {
int64_t size;
double* restrict data;
int64_t size;
} qmckl_vector;
#+end_src
@ -160,8 +160,8 @@ qmckl_vector_free( qmckl_context context,
#+begin_src c :comments org :tangle (eval h_private_type) :exports none
typedef struct qmckl_matrix {
int64_t size[2];
double* restrict data;
int64_t size[2];
} qmckl_matrix;
#+end_src
@ -245,9 +245,9 @@ qmckl_matrix_free( qmckl_context context,
#define QMCKL_TENSOR_ORDER_MAX 16
typedef struct qmckl_tensor {
double* restrict data;
int64_t order;
int64_t size[QMCKL_TENSOR_ORDER_MAX];
double* restrict data;
} qmckl_tensor;
#+end_src

4
org/qmckl_context.org
View File

@ -169,7 +169,7 @@ qmckl_context qmckl_context_check(const qmckl_context context) {
if (context == QMCKL_NULL_CONTEXT)
return QMCKL_NULL_CONTEXT;
const qmckl_context_struct* const ctx = (const qmckl_context_struct*) context;
const qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
/* Try to access memory */
if (ctx->tag != VALID_TAG) {
@ -267,7 +267,7 @@ qmckl_context qmckl_context_create() {
{
ctx->tag = VALID_TAG;
const qmckl_context context = (const qmckl_context) ctx;
const qmckl_context context = (qmckl_context) ctx;
assert ( qmckl_context_check(context) != QMCKL_NULL_CONTEXT );
qmckl_exit_code rc;

113
org/qmckl_determinant.org
View File

@ -182,7 +182,7 @@ qmckl_exit_code qmckl_init_determinant(qmckl_context context) {
return false;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
ctx->det.uninitialized = (1 << 6) - 1;
@ -216,7 +216,7 @@ bool qmckl_determinant_provided(const qmckl_context context) {
return false;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
return ctx->det.provided;
@ -238,7 +238,7 @@ char qmckl_get_determinant_type (const qmckl_context context) {
return (char) 0;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int32_t mask = 1;
@ -256,7 +256,7 @@ int64_t qmckl_get_determinant_walk_num (const qmckl_context context) {
return (int64_t) 0;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int32_t mask = 1 << 1;
@ -274,7 +274,7 @@ int64_t qmckl_get_determinant_det_num_alpha (const qmckl_context context) {
return (int64_t) 0;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int32_t mask = 1 << 2;
@ -292,7 +292,7 @@ int64_t qmckl_get_determinant_det_num_beta (const qmckl_context context) {
return (int64_t) 0;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int32_t mask = 1 << 3;
@ -310,7 +310,7 @@ int64_t* qmckl_get_determinant_mo_index_alpha (const qmckl_context context) {
return NULL;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int32_t mask = 1 << 4;
@ -328,7 +328,7 @@ int64_t* qmckl_get_determinant_mo_index_beta (const qmckl_context context) {
return NULL;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int32_t mask = 1 << 5;
@ -363,7 +363,7 @@ if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
#+end_src
#+NAME:post2
@ -525,7 +525,7 @@ qmckl_exit_code qmckl_finalize_determinant(qmckl_context context) {
NULL);
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
qmckl_exit_code rc;
@ -596,7 +596,7 @@ qmckl_exit_code qmckl_get_det_vgl_alpha(qmckl_context context, double * const de
rc = qmckl_provide_det_vgl_alpha(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = 5 * ctx->det.det_num_alpha * ctx->det.walk_num *
@ -623,7 +623,7 @@ qmckl_exit_code qmckl_get_det_vgl_beta(qmckl_context context, double * const det
rc = qmckl_provide_det_vgl_beta(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = 5 * ctx->det.det_num_beta * ctx->det.walk_num *
@ -649,7 +649,7 @@ qmckl_exit_code qmckl_provide_det_vgl_alpha(qmckl_context context) {
return QMCKL_NULL_CONTEXT;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if(!(ctx->nucleus.provided)) {
@ -748,7 +748,7 @@ qmckl_exit_code qmckl_provide_det_vgl_beta(qmckl_context context) {
return QMCKL_NULL_CONTEXT;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if(!(ctx->nucleus.provided)) {
@ -1134,36 +1134,28 @@ end function qmckl_compute_det_vgl_beta_f
#+begin_src c :tangle (eval c_test) :exports none
#define walk_num chbrclf_walk_num
#define elec_num chbrclf_elec_num
#define shell_num chbrclf_shell_num
#define ao_num chbrclf_ao_num
int64_t elec_up_num = chbrclf_elec_up_num;
int64_t elec_dn_num = chbrclf_elec_dn_num;
double* elec_coord = &(chbrclf_elec_coord[0][0][0]);
const int64_t nucl_num = chbrclf_nucl_num;
const double* nucl_charge = chbrclf_charge;
const double* nucl_coord = &(chbrclf_nucl_coord[0][0]);
rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num);
rc = qmckl_set_electron_num (context, chbrclf_elec_up_num, chbrclf_elec_dn_num);
assert (rc == QMCKL_SUCCESS);
rc = qmckl_set_electron_walk_num (context, walk_num);
rc = qmckl_set_electron_walk_num (context, chbrclf_walk_num);
assert (rc == QMCKL_SUCCESS);
assert(qmckl_electron_provided(context));
rc = qmckl_set_electron_coord (context, 'N', elec_coord, walk_num*elec_num*3);
rc = qmckl_set_electron_coord (context, 'N', elec_coord, chbrclf_walk_num*chbrclf_elec_num*3);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_num (context, nucl_num);
rc = qmckl_set_nucleus_num (context, chbrclf_nucl_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), nucl_num*3);
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), chbrclf_nucl_num*3);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_charge(context, nucl_charge, nucl_num);
rc = qmckl_set_nucleus_charge(context, nucl_charge, chbrclf_nucl_num);
assert(rc == QMCKL_SUCCESS);
assert(qmckl_nucleus_provided(context));
@ -1195,27 +1187,27 @@ rc = qmckl_set_ao_basis_prim_num (context, chbrclf_prim_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_nucleus_index (context, nucleus_index, nucl_num);
rc = qmckl_set_ao_basis_nucleus_index (context, nucleus_index, chbrclf_nucl_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_nucleus_shell_num (context, nucleus_shell_num, nucl_num);
rc = qmckl_set_ao_basis_nucleus_shell_num (context, nucleus_shell_num, chbrclf_nucl_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_shell_ang_mom (context, shell_ang_mom, shell_num);
rc = qmckl_set_ao_basis_shell_ang_mom (context, shell_ang_mom, chbrclf_shell_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_shell_factor (context, shell_factor, shell_num);
rc = qmckl_set_ao_basis_shell_factor (context, shell_factor, chbrclf_shell_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_shell_prim_num (context, shell_prim_num, shell_num);
rc = qmckl_set_ao_basis_shell_prim_num (context, shell_prim_num, chbrclf_shell_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_shell_prim_index (context, shell_prim_index, shell_num);
rc = qmckl_set_ao_basis_shell_prim_index (context, shell_prim_index, chbrclf_shell_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
@ -1239,14 +1231,13 @@ assert(rc == QMCKL_SUCCESS);
assert(qmckl_ao_basis_provided(context));
double ao_vgl[walk_num*elec_num][5][chbrclf_ao_num];
double ao_vgl[chbrclf_walk_num*chbrclf_elec_num][5][chbrclf_ao_num];
rc = qmckl_get_ao_basis_ao_vgl(context, &(ao_vgl[0][0][0]), (int64_t) 5*walk_num*elec_num*chbrclf_ao_num);
rc = qmckl_get_ao_basis_ao_vgl(context, &(ao_vgl[0][0][0]), (int64_t) 5*chbrclf_walk_num*chbrclf_elec_num*chbrclf_ao_num);
assert (rc == QMCKL_SUCCESS);
/* Set up MO data */
const int64_t mo_num = chbrclf_mo_num;
rc = qmckl_set_mo_basis_mo_num(context, mo_num);
rc = qmckl_set_mo_basis_mo_num(context, chbrclf_mo_num);
assert (rc == QMCKL_SUCCESS);
const double * mo_coefficient = &(chbrclf_mo_coef[0]);
@ -1256,31 +1247,31 @@ assert (rc == QMCKL_SUCCESS);
assert(qmckl_mo_basis_provided(context));
double mo_vgl[walk_num*elec_num][5][chbrclf_mo_num];
rc = qmckl_get_mo_basis_mo_vgl(context, &(mo_vgl[0][0][0]), 5*walk_num*elec_num*chbrclf_mo_num);
double mo_vgl[chbrclf_walk_num*chbrclf_elec_num][5][chbrclf_mo_num];
rc = qmckl_get_mo_basis_mo_vgl(context, &(mo_vgl[0][0][0]), 5*chbrclf_walk_num*chbrclf_elec_num*chbrclf_mo_num);
assert (rc == QMCKL_SUCCESS);
/* Set up determinant data */
const int64_t det_num_alpha = 1;
const int64_t det_num_beta = 1;
int64_t mo_index_alpha[det_num_alpha][walk_num][elec_up_num];
int64_t mo_index_beta[det_num_alpha][walk_num][elec_dn_num];
#define det_num_alpha 1
#define det_num_beta 1
int64_t mo_index_alpha[det_num_alpha][chbrclf_walk_num][chbrclf_elec_up_num];
int64_t mo_index_beta[det_num_alpha][chbrclf_walk_num][chbrclf_elec_dn_num];
int i, j, k;
for(k = 0; k < det_num_alpha; ++k)
for(i = 0; i < walk_num; ++i)
for(j = 0; j < elec_up_num; ++j)
for(i = 0; i < chbrclf_walk_num; ++i)
for(j = 0; j < chbrclf_elec_up_num; ++j)
mo_index_alpha[k][i][j] = j + 1;
for(k = 0; k < det_num_beta; ++k)
for(i = 0; i < walk_num; ++i)
for(j = 0; j < elec_up_num; ++j)
for(i = 0; i < chbrclf_walk_num; ++i)
for(j = 0; j < chbrclf_elec_up_num; ++j)
mo_index_beta[k][i][j] = j + 1;
rc = qmckl_set_determinant_type (context, typ);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_determinant_walk_num (context, walk_num);
rc = qmckl_set_determinant_walk_num (context, chbrclf_walk_num);
assert (rc == QMCKL_SUCCESS);
rc = qmckl_set_determinant_det_num_alpha (context, det_num_alpha);
@ -1297,8 +1288,8 @@ assert (rc == QMCKL_SUCCESS);
// Get slater-determinant
double det_vgl_alpha[det_num_alpha][walk_num][5][elec_up_num][elec_up_num];
double det_vgl_beta[det_num_beta][walk_num][5][elec_dn_num][elec_dn_num];
double det_vgl_alpha[det_num_alpha][chbrclf_walk_num][5][chbrclf_elec_up_num][chbrclf_elec_up_num];
double det_vgl_beta[det_num_beta][chbrclf_walk_num][5][chbrclf_elec_dn_num][chbrclf_elec_dn_num];
rc = qmckl_get_det_vgl_alpha(context, &(det_vgl_alpha[0][0][0][0][0]));
assert (rc == QMCKL_SUCCESS);
@ -1347,7 +1338,7 @@ qmckl_exit_code qmckl_get_det_inv_matrix_alpha(qmckl_context context, double * c
rc = qmckl_provide_det_inv_matrix_alpha(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = ctx->det.det_num_alpha * ctx->det.walk_num * ctx->electron.up_num * ctx->electron.up_num;
@ -1376,7 +1367,7 @@ qmckl_exit_code qmckl_get_det_inv_matrix_beta(qmckl_context context, double * co
rc = qmckl_provide_det_inv_matrix_beta(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = ctx->det.det_num_alpha * ctx->det.walk_num * ctx->electron.down_num * ctx->electron.down_num;
@ -1405,7 +1396,7 @@ qmckl_exit_code qmckl_get_det_adj_matrix_alpha(qmckl_context context, double * c
rc = qmckl_provide_det_inv_matrix_alpha(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = ctx->det.det_num_alpha * ctx->det.walk_num * ctx->electron.up_num * ctx->electron.up_num;
@ -1434,7 +1425,7 @@ qmckl_exit_code qmckl_get_det_adj_matrix_beta(qmckl_context context, double * co
rc = qmckl_provide_det_inv_matrix_beta(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = ctx->det.det_num_alpha * ctx->det.walk_num * ctx->electron.down_num * ctx->electron.down_num;
@ -1463,7 +1454,7 @@ qmckl_exit_code qmckl_get_det_alpha(qmckl_context context, double * const det_va
rc = qmckl_provide_det_inv_matrix_alpha(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = ctx->det.det_num_alpha * ctx->det.walk_num;
@ -1492,7 +1483,7 @@ qmckl_exit_code qmckl_get_det_beta(qmckl_context context, double * const det_val
rc = qmckl_provide_det_inv_matrix_beta(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = ctx->det.det_num_alpha * ctx->det.walk_num;
@ -1517,7 +1508,7 @@ qmckl_exit_code qmckl_provide_det_inv_matrix_alpha(qmckl_context context) {
return QMCKL_NULL_CONTEXT;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if(!(ctx->nucleus.provided)) {
@ -1640,7 +1631,7 @@ qmckl_exit_code qmckl_provide_det_inv_matrix_beta(qmckl_context context) {
return QMCKL_NULL_CONTEXT;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if(!(ctx->nucleus.provided)) {
@ -2047,8 +2038,8 @@ end function qmckl_compute_det_inv_matrix_beta_f
#+begin_src c :tangle (eval c_test) :exports none
// Get adjoint of the slater-determinant
double det_inv_matrix_alpha[det_num_alpha][walk_num][elec_up_num][elec_up_num];
double det_inv_matrix_beta[det_num_beta][walk_num][elec_dn_num][elec_dn_num];
double det_inv_matrix_alpha[det_num_alpha][chbrclf_walk_num][chbrclf_elec_up_num][chbrclf_elec_up_num];
double det_inv_matrix_beta[det_num_beta][chbrclf_walk_num][chbrclf_elec_dn_num][chbrclf_elec_dn_num];
rc = qmckl_get_det_inv_matrix_alpha(context, &(det_inv_matrix_alpha[0][0][0][0]));
assert (rc == QMCKL_SUCCESS);

6
org/qmckl_distance.org
View File

@ -133,7 +133,7 @@ integer function qmckl_distance_sq_f(context, transa, transb, m, n, &
if (transb == 'N' .or. transb == 'n') then
continue
else if (transa == 'T' .or. transa == 't') then
else if (transb == 'T' .or. transb == 't') then
transab = transab + 2
else
transab = -100
@ -533,7 +533,7 @@ integer function qmckl_distance_f(context, transa, transb, m, n, &
if (transb == 'N' .or. transb == 'n') then
continue
else if (transa == 'T' .or. transa == 't') then
else if (transb == 'T' .or. transb == 't') then
transab = transab + 2
else
transab = -100
@ -1314,7 +1314,7 @@ integer function qmckl_distance_rescaled_deriv_e_f(context, transa, transb, m, n
if (transb == 'N' .or. transb == 'n') then
continue
else if (transa == 'T' .or. transa == 't') then
else if (transb == 'T' .or. transb == 't') then
transab = transab + 2
else
transab = -100

82
org/qmckl_electron.org
View File

@ -97,8 +97,8 @@ int main() {
| ~ee_distance_rescaled_deriv_e_date~ | ~uint64_t~ | Last modification date of the electron-electron distance derivatives |
| ~ee_pot~ | ~double[walk_num]~ | Electron-electron rescaled distances derivatives |
| ~ee_pot_date~ | ~uint64_t~ | Last modification date of the electron-electron distance derivatives |
| ~en_pot~ | double[walk_num] | Electron-nucleus potential energy |
| ~en_pot_date~ | int64_t | Date when the electron-nucleus potential energy was computed |
| ~en_pot~ | ~double[walk_num]~ | Electron-nucleus potential energy |
| ~en_pot_date~ | ~int64_t~ | Date when the electron-nucleus potential energy was computed |
| ~en_distance_rescaled~ | ~double[walk_num][nucl_num][num]~ | Electron-nucleus distances |
| ~en_distance_rescaled_date~ | ~uint64_t~ | Last modification date of the electron-electron distances |
| ~en_distance_rescaled_deriv_e~ | ~double[walk_num][4][nucl_num][num]~ | Electron-electron rescaled distances derivatives |
@ -157,7 +157,7 @@ qmckl_exit_code qmckl_init_electron(qmckl_context context) {
return false;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
ctx->electron.uninitialized = (1 << 2) - 1;
@ -182,7 +182,7 @@ bool qmckl_electron_provided(const qmckl_context context) {
return false;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
return ctx->electron.provided;
@ -228,7 +228,7 @@ qmckl_get_electron_num (const qmckl_context context, int64_t* const num) {
"num is a null pointer");
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int32_t mask = 1 << 0;
@ -256,7 +256,7 @@ qmckl_get_electron_up_num (const qmckl_context context, int64_t* const up_num) {
"up_num is a null pointer");
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int32_t mask = 1 << 0;
@ -284,7 +284,7 @@ qmckl_get_electron_down_num (const qmckl_context context, int64_t* const down_nu
"down_num is a null pointer");
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int32_t mask = 1 << 0;
@ -323,7 +323,7 @@ qmckl_get_electron_walk_num (const qmckl_context context, int64_t* const walk_nu
"walk_num is a null pointer");
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int32_t mask = 1 << 1;
@ -360,7 +360,7 @@ qmckl_get_electron_rescale_factor_ee (const qmckl_context context, double* const
"rescale_factor_kappa_ee is a null pointer");
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
assert (ctx->electron.rescale_factor_kappa_ee > 0.0);
@ -383,7 +383,7 @@ qmckl_get_electron_rescale_factor_en (const qmckl_context context, double* const
"rescale_factor_kappa_en is a null pointer");
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
assert (ctx->electron.rescale_factor_kappa_en > 0.0);
@ -448,7 +448,7 @@ qmckl_get_electron_coord (const qmckl_context context,
return QMCKL_INVALID_CONTEXT;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if (!ctx->electron.provided) {
@ -489,7 +489,7 @@ if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
#+end_src
#+NAME:post2
@ -718,7 +718,7 @@ qmckl_set_electron_coord(qmckl_context context,
ctx->electron.coord_old = ctx->electron.coord_new ;
qmckl_exit_code rc;
rc = qmckl_set_point(context, transp, coord, size_max/3);
rc = qmckl_set_point(context, transp, size_max/3, coord, size_max);
assert (rc == QMCKL_SUCCESS);
ctx->electron.coord_new = ctx->point.coord ;
@ -897,7 +897,7 @@ qmckl_exit_code qmckl_get_electron_ee_distance(qmckl_context context, double* co
rc = qmckl_provide_ee_distance(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = ctx->electron.num * ctx->electron.num * ctx->electron.walk_num;
@ -921,7 +921,7 @@ qmckl_exit_code qmckl_provide_ee_distance(qmckl_context context)
return QMCKL_NULL_CONTEXT;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
@ -1138,7 +1138,7 @@ qmckl_exit_code qmckl_get_electron_ee_distance_rescaled(qmckl_context context, d
rc = qmckl_provide_ee_distance_rescaled(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = ctx->electron.num * ctx->electron.num * ctx->electron.walk_num;
@ -1162,7 +1162,7 @@ qmckl_exit_code qmckl_provide_ee_distance_rescaled(qmckl_context context)
return QMCKL_NULL_CONTEXT;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
@ -1218,7 +1218,7 @@ qmckl_exit_code qmckl_provide_ee_distance_rescaled(qmckl_context context)
| ~elec_num~ | ~int64_t~ | in | Number of electrons |
| ~rescale_factor_kappa_ee~ | ~double~ | in | Factor to rescale ee distances |
| ~walk_num~ | ~int64_t~ | in | Number of walkers |
| ~coord~ | ~double[walk_num][3][elec_num]~ | in | Electron coordinates |
| ~coord~ | ~double[3][walk_num][elec_num]~ | in | Electron coordinates |
| ~ee_distance~ | ~double[walk_num][elec_num][elec_num]~ | out | Electron-electron rescaled distances |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
@ -1231,7 +1231,7 @@ integer function qmckl_compute_ee_distance_rescaled_f(context, elec_num, rescale
integer*8 , intent(in) :: elec_num
double precision , intent(in) :: rescale_factor_kappa_ee
integer*8 , intent(in) :: walk_num
double precision , intent(in) :: coord(elec_num,3,walk_num)
double precision , intent(in) :: coord(elec_num,walk_num,3)
double precision , intent(out) :: ee_distance_rescaled(elec_num,elec_num,walk_num)
integer*8 :: k
@ -1357,7 +1357,7 @@ assert(fabs(ee_distance_rescaled[elec_num*elec_num+1]-0.9985724058042633) < 1.e-
#+end_src
** Electron-electron rescaled distance gradients and laplacian with respect to electron coords
** Electron-electron rescaled distance gradients and Laplacian with respect to electron coords
The rescaled distances which is given as $R = (1 - \exp{-\kappa r})/\kappa$
needs to be perturbed with respect to the electorn coordinates.
@ -1384,7 +1384,7 @@ qmckl_exit_code qmckl_get_electron_ee_distance_rescaled_deriv_e(qmckl_context co
rc = qmckl_provide_ee_distance_rescaled_deriv_e(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = 4 * ctx->electron.num * ctx->electron.num * ctx->electron.walk_num;
@ -1408,7 +1408,7 @@ qmckl_exit_code qmckl_provide_ee_distance_rescaled_deriv_e(qmckl_context context
return QMCKL_NULL_CONTEXT;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
@ -1464,7 +1464,7 @@ qmckl_exit_code qmckl_provide_ee_distance_rescaled_deriv_e(qmckl_context context
| ~elec_num~ | ~int64_t~ | in | Number of electrons |
| ~rescale_factor_kappa_ee~ | ~double~ | in | Factor to rescale ee distances |
| ~walk_num~ | ~int64_t~ | in | Number of walkers |
| ~coord~ | ~double[walk_num][3][elec_num]~ | in | Electron coordinates |
| ~coord~ | ~double[3][walk_num][elec_num]~ | in | Electron coordinates |
| ~ee_distance_deriv_e~ | ~double[walk_num][4][elec_num][elec_num]~ | out | Electron-electron rescaled distance derivatives |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
@ -1477,7 +1477,7 @@ integer function qmckl_compute_ee_distance_rescaled_deriv_e_f(context, elec_num,
integer*8 , intent(in) :: elec_num
double precision , intent(in) :: rescale_factor_kappa_ee
integer*8 , intent(in) :: walk_num
double precision , intent(in) :: coord(elec_num,3,walk_num)
double precision , intent(in) :: coord(elec_num,walk_num,3)
double precision , intent(out) :: ee_distance_rescaled_deriv_e(4,elec_num,elec_num,walk_num)
integer*8 :: k
@ -1501,8 +1501,8 @@ integer function qmckl_compute_ee_distance_rescaled_deriv_e_f(context, elec_num,
do k=1,walk_num
info = qmckl_distance_rescaled_deriv_e(context, 'T', 'T', elec_num, elec_num, &
coord(1,1,k), elec_num, &
coord(1,1,k), elec_num, &
coord(1,k,1), elec_num*walk_num, &
coord(1,k,1), elec_num*walk_num, &
ee_distance_rescaled_deriv_e(1,1,1,k), elec_num, rescale_factor_kappa_ee)
if (info /= QMCKL_SUCCESS) then
exit
@ -1613,7 +1613,7 @@ qmckl_exit_code qmckl_get_electron_ee_potential(qmckl_context context, double* c
rc = qmckl_provide_ee_potential(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = ctx->electron.walk_num * sizeof(double);
@ -1637,7 +1637,7 @@ qmckl_exit_code qmckl_provide_ee_potential(qmckl_context context)
return QMCKL_NULL_CONTEXT;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if (!ctx->electron.provided) return QMCKL_NOT_PROVIDED;
@ -1818,7 +1818,7 @@ qmckl_exit_code qmckl_get_electron_en_distance(qmckl_context context, double* di
rc = qmckl_provide_en_distance(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num;
@ -1842,7 +1842,7 @@ qmckl_exit_code qmckl_provide_en_distance(qmckl_context context)
return QMCKL_NULL_CONTEXT;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if (!(ctx->nucleus.provided)) {
@ -1905,7 +1905,7 @@ qmckl_exit_code qmckl_provide_en_distance(qmckl_context context)
| ~elec_num~ | ~int64_t~ | in | Number of electrons |
| ~nucl_num~ | ~int64_t~ | in | Number of nuclei |
| ~walk_num~ | ~int64_t~ | in | Number of walkers |
| ~elec_coord~ | ~double[walk_num][3][elec_num]~ | in | Electron coordinates |
| ~elec_coord~ | ~double[3][walk_num][elec_num]~ | in | Electron coordinates |
| ~nucl_coord~ | ~double[3][elec_num]~ | in | Nuclear coordinates |
| ~en_distance~ | ~double[walk_num][nucl_num][elec_num]~ | out | Electron-nucleus distances |
@ -2097,7 +2097,7 @@ qmckl_exit_code qmckl_get_electron_en_distance_rescaled(qmckl_context context, d
rc = qmckl_provide_en_distance_rescaled(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num;
@ -2121,7 +2121,7 @@ qmckl_exit_code qmckl_provide_en_distance_rescaled(qmckl_context context)
return QMCKL_NULL_CONTEXT;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if (!(ctx->nucleus.provided)) {
@ -2183,7 +2183,7 @@ qmckl_exit_code qmckl_provide_en_distance_rescaled(qmckl_context context)
| ~nucl_num~ | ~int64_t~ | in | Number of nuclei |
| ~rescale_factor_kappa_en~ | ~double~ | in | The factor for rescaled distances |
| ~walk_num~ | ~int64_t~ | in | Number of walkers |
| ~elec_coord~ | ~double[walk_num][3][elec_num]~ | in | Electron coordinates |
| ~elec_coord~ | ~double[3][walk_num][elec_num]~ | in | Electron coordinates |
| ~nucl_coord~ | ~double[3][elec_num]~ | in | Nuclear coordinates |
| ~en_distance_rescaled~ | ~double[walk_num][nucl_num][elec_num]~ | out | Electron-nucleus distances |
@ -2385,7 +2385,7 @@ qmckl_exit_code qmckl_get_electron_en_distance_rescaled_deriv_e(qmckl_context co
rc = qmckl_provide_en_distance_rescaled_deriv_e(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = 4 * ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num;
@ -2409,7 +2409,7 @@ qmckl_exit_code qmckl_provide_en_distance_rescaled_deriv_e(qmckl_context context
return QMCKL_NULL_CONTEXT;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if (!(ctx->nucleus.provided)) {
@ -2471,9 +2471,9 @@ qmckl_exit_code qmckl_provide_en_distance_rescaled_deriv_e(qmckl_context context
| ~nucl_num~ | ~int64_t~ | in | Number of nuclei |
| ~rescale_factor_kappa_en~ | ~double~ | in | The factor for rescaled distances |
| ~walk_num~ | ~int64_t~ | in | Number of walkers |
| ~elec_coord~ | ~double[walk_num][3][elec_num]~ | in | Electron coordinates |
| ~elec_coord~ | ~double[3][walk_num][elec_num]~ | in | Electron coordinates |
| ~nucl_coord~ | ~double[3][elec_num]~ | in | Nuclear coordinates |
| ~en_distance_rescaled_deriv_e~ | ~double[walk_num][4][nucl_num][elec_num]~ | out | Electron-nucleus distance derivatives |
| ~en_distance_rescaled_deriv_e~ | ~double[walk_num][nucl_num][elec_num][4]~ | out | Electron-nucleus distance derivatives |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_en_distance_rescaled_deriv_e_f(context, elec_num, nucl_num, &
@ -2656,7 +2656,7 @@ qmckl_exit_code qmckl_get_electron_en_potential(qmckl_context context, double* c
rc = qmckl_provide_en_potential(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = ctx->electron.walk_num * sizeof(double);
@ -2680,7 +2680,7 @@ qmckl_exit_code qmckl_provide_en_potential(qmckl_context context)
return QMCKL_NULL_CONTEXT;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if (!ctx->electron.provided) return QMCKL_NOT_PROVIDED;
@ -2843,7 +2843,7 @@ assert (rc == QMCKL_SUCCESS);
*** Compute :noexport:
# begin_src f90 :comments org :tangle (eval f) :noweb yes
# begin_src f90 :comments org :tangle (eval f) :noweb yes
subroutine draw_init_points
implicit none
BEGIN_DOC

120
org/qmckl_error.org
View File

@ -239,7 +239,7 @@ for (text, code, message) in table:
message = message.replace("'",'"')
result += [ f"""case {text}:
return {message};
break;""" ]
""" ]
return '\n'.join(result)
#+end_src
@ -247,89 +247,91 @@ return '\n'.join(result)
#+RESULTS: cases
#+begin_example
case QMCKL_SUCCESS:
return "Success";
break;
return "Success";
case QMCKL_INVALID_ARG_1:
return "Invalid argument 1";
break;
return "Invalid argument 1";
case QMCKL_INVALID_ARG_2:
return "Invalid argument 2";
break;
return "Invalid argument 2";
case QMCKL_INVALID_ARG_3:
return "Invalid argument 3";
break;
return "Invalid argument 3";
case QMCKL_INVALID_ARG_4:
return "Invalid argument 4";
break;
return "Invalid argument 4";
case QMCKL_INVALID_ARG_5:
return "Invalid argument 5";
break;
return "Invalid argument 5";
case QMCKL_INVALID_ARG_6:
return "Invalid argument 6";
break;
return "Invalid argument 6";
case QMCKL_INVALID_ARG_7:
return "Invalid argument 7";
break;
return "Invalid argument 7";
case QMCKL_INVALID_ARG_8:
return "Invalid argument 8";
break;
return "Invalid argument 8";
case QMCKL_INVALID_ARG_9:
return "Invalid argument 9";
break;
return "Invalid argument 9";
case QMCKL_INVALID_ARG_10:
return "Invalid argument 10";
break;
return "Invalid argument 10";
case QMCKL_INVALID_ARG_11:
return "Invalid argument 11";
break;
return "Invalid argument 11";
case QMCKL_INVALID_ARG_12:
return "Invalid argument 12";
break;
return "Invalid argument 12";
case QMCKL_INVALID_ARG_13:
return "Invalid argument 13";
break;
return "Invalid argument 13";
case QMCKL_INVALID_ARG_14:
return "Invalid argument 14";
break;
return "Invalid argument 14";
case QMCKL_INVALID_ARG_15:
return "Invalid argument 15";
break;
return "Invalid argument 15";
case QMCKL_INVALID_ARG_16:
return "Invalid argument 16";
break;
return "Invalid argument 16";
case QMCKL_INVALID_ARG_17:
return "Invalid argument 17";
break;
return "Invalid argument 17";
case QMCKL_INVALID_ARG_18:
return "Invalid argument 18";
break;
return "Invalid argument 18";
case QMCKL_INVALID_ARG_19:
return "Invalid argument 19";
break;
return "Invalid argument 19";
case QMCKL_INVALID_ARG_20:
return "Invalid argument 20";
break;
return "Invalid argument 20";
case QMCKL_FAILURE:
return "Failure";
break;
return "Failure";
case QMCKL_ERRNO:
return strerror(errno);
break;
return strerror(errno);
case QMCKL_INVALID_CONTEXT:
return "Invalid context";
break;
return "Invalid context";
case QMCKL_ALLOCATION_FAILED:
return "Allocation failed";
break;
return "Allocation failed";
case QMCKL_DEALLOCATION_FAILED:
return "De-allocation failed";
break;
return "De-allocation failed";
case QMCKL_NOT_PROVIDED:
return "Not provided";
break;
return "Not provided";
case QMCKL_OUT_OF_BOUNDS:
return "Index out of bounds";
case QMCKL_INVALID_EXIT_CODE:
return "Invalid exit code";
break;
return "Invalid exit code";
#+end_example
# Source
@ -414,7 +416,7 @@ qmckl_set_error(qmckl_context context,
qmckl_lock(context);
{
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL); /* Impossible because the context is valid. */
ctx->error.exit_code = exit_code;
@ -460,7 +462,7 @@ qmckl_get_error(qmckl_context context,
qmckl_lock(context);
{
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL); /* Impossible because the context is valid. */
/* Turn off annoying GCC warning */

360
org/qmckl_examples.org Normal file
View File

@ -0,0 +1,360 @@
#+TITLE: Code examples
#+SETUPFILE: ../tools/theme.setup
#+INCLUDE: ../tools/lib.org
In this section, we present examples of usage of QMCkl.
For simplicity, we assume that the wave function parameters are stored
in a [[https://github.com/TREX-CoE/trexio][TREXIO]] file.
* Python
** Check numerically that MOs are orthonormal
In this example, we will compute numerically the overlap
between the molecular orbitals:
\[
S_{ij} = \int \phi_i(\mathbf{r}) \phi_j(\mathbf{r})
\text{d}\mathbf{r} \sim \sum_{k=1}^{N} \phi_i(\mathbf{r}_k)
\phi_j(\mathbf{r}_k) \delta \mathbf{r}
\]
\[
S_{ij} = \langle \phi_i | \phi_j \rangle
\sim \sum_{k=1}^{N} \langle \phi_i | \mathbf{r}_k \rangle
\langle \mathbf{r}_k | \phi_j \rangle
\]
#+begin_src python :exports code
import numpy as np
import qmckl
#+end_src
#+RESULTS:
First, we create a context for the QMCkl calculation, and load the
wave function stored in =h2o_5z.h5= inside it. It is a Hartree-Fock
determinant for the water molecule in the cc-pV5Z basis set.
#+begin_src python :exports code
trexio_filename = "..//share/qmckl/test_data/h2o_5z.h5"
context = qmckl.context_create()
qmckl.trexio_read(context, trexio_filename)
#+end_src
#+RESULTS:
: None
We now define the grid points $\mathbf{r}_k$ as a regular grid around the
molecule.
We fetch the nuclear coordinates from the context,
#+begin_src python :exports code
nucl_num = qmckl.get_nucleus_num(context)
nucl_charge = qmckl.get_nucleus_charge(context, nucl_num)
nucl_coord = qmckl.get_nucleus_coord(context, 'N', nucl_num*3)
nucl_coord = np.reshape(nucl_coord, (3, nucl_num))
for i in range(nucl_num):
print("%d %+f %+f %+f"%(int(nucl_charge[i]),
nucl_coord[i,0],
nucl_coord[i,1],
nucl_coord[i,2]) )
#+end_src
#+begin_example
8 +0.000000 +0.000000 +0.000000
1 -1.430429 +0.000000 -1.107157
1 +1.430429 +0.000000 -1.107157
#+end_example
and compute the coordinates of the grid points:
#+begin_src python :exports code
nx = ( 120, 120, 120 )
shift = np.array([5.,5.,5.])
point_num = nx[0] * nx[1] * nx[2]
rmin = np.array( list([ np.min(nucl_coord[:,a]) for a in range(3) ]) )
rmax = np.array( list([ np.max(nucl_coord[:,a]) for a in range(3) ]) )
linspace = [ None for i in range(3) ]
step = [ None for i in range(3) ]
for a in range(3):
linspace[a], step[a] = np.linspace(rmin[a]-shift[a],
rmax[a]+shift[a],
num=nx[a],
retstep=True)
dr = step[0] * step[1] * step[2]
#+end_src
#+RESULTS:
Now the grid is ready, we can create the list of grid points
$\mathbf{r}_k$ on which the MOs $\phi_i$ will be evaluated, and
transfer them to the QMCkl context:
#+begin_src python :exports code
point = []
for x in linspace[0]:
for y in linspace[1]:
for z in linspace[2]:
point += [ [x, y, z] ]
point = np.array(point)
point_num = len(point)
qmckl.set_point(context, 'N', point_num, np.reshape(point, (point_num*3)))
#+end_src
#+RESULTS:
: None
Then, we evaluate all the MOs at the grid points (and time the execution),
and thus obtain the matrix $M_{ki} = \langle \mathbf{r}_k | \phi_i \rangle =
\phi_i(\mathbf{r}_k)$.
#+begin_src python :exports code
import time
mo_num = qmckl.get_mo_basis_mo_num(context)
before = time.time()
mo_value = qmckl.get_mo_basis_mo_value(context, point_num*mo_num)
after = time.time()
mo_value = np.reshape( mo_value, (point_num, mo_num) )
print("Number of MOs: ", mo_num)
print("Number of grid points: ", point_num)
print("Execution time : ", (after - before), "seconds")
#+end_src
#+begin_example
Number of MOs: 201
Number of grid points: 1728000
Execution time : 3.511528968811035 seconds
#+end_example
and finally we compute the overlap between all the MOs as
$M^\dagger M$.
#+begin_src python :exports code
overlap = mo_value.T @ mo_value * dr
print (overlap)
#+end_src
#+begin_example
[[ 9.88693941e-01 2.34719693e-03 -1.50518232e-08 ... 3.12084178e-09
-5.81064929e-10 3.70130091e-02]
[ 2.34719693e-03 9.99509628e-01 3.18930040e-09 ... -2.46888958e-10
-1.06064273e-09 -7.65567973e-03]
[-1.50518232e-08 3.18930040e-09 9.99995073e-01 ... -5.84882580e-06
-1.21598117e-06 4.59036468e-08]
...
[ 3.12084178e-09 -2.46888958e-10 -5.84882580e-06 ... 1.00019107e+00
-2.03342837e-04 -1.36954855e-08]
[-5.81064929e-10 -1.06064273e-09 -1.21598117e-06 ... -2.03342837e-04
9.99262427e-01 1.18264754e-09]
[ 3.70130091e-02 -7.65567973e-03 4.59036468e-08 ... -1.36954855e-08
1.18264754e-09 8.97215950e-01]]
#+end_example
* Fortran
** Checking errors
All QMCkl functions return an error code. A convenient way to handle
errors is to write an error-checking function that displays the
error in text format and exits the program.
#+NAME: qmckl_check_error
#+begin_src f90
subroutine qmckl_check_error(rc, message)
use qmckl
implicit none
integer(qmckl_exit_code), intent(in) :: rc
character(len=*) , intent(in) :: message
character(len=128) :: str_buffer
if (rc /= QMCKL_SUCCESS) then
print *, message
call qmckl_string_of_error(rc, str_buffer)
print *, str_buffer
call exit(rc)
end if
end subroutine qmckl_check_error
#+end_src
** Computing an atomic orbital on a grid
:PROPERTIES:
:header-args: :tangle ao_grid.f90
:END:
The following program, in Fortran, computes the values of an atomic
orbital on a regular 3-dimensional grid. The 100^3 grid points are
automatically defined, such that the molecule fits in a box with 5
atomic units in the borders.
This program uses the ~qmckl_check_error~ function defined above.
To use this program, run
#+begin_src bash :tangle no :exports code
$ ao_grid <trexio_file> <AO_id> <point_num>
#+end_src
#+begin_src f90 :noweb yes
<<qmckl_check_error>>
program ao_grid
use qmckl
implicit none
integer(qmckl_context) :: qmckl_ctx ! QMCkl context
integer(qmckl_exit_code) :: rc ! Exit code of QMCkl functions
character(len=128) :: trexio_filename
character(len=128) :: str_buffer
integer :: ao_id
integer :: point_num_x
integer(c_int64_t) :: nucl_num
double precision, allocatable :: nucl_coord(:,:)
integer(c_int64_t) :: point_num
integer(c_int64_t) :: ao_num
integer(c_int64_t) :: ipoint, i, j, k
double precision :: x, y, z, dr(3)
double precision :: rmin(3), rmax(3)
double precision, allocatable :: points(:,:)
double precision, allocatable :: ao_vgl(:,:,:)
#+end_src
Start by fetching the command-line arguments:
#+begin_src f90
if (iargc() /= 3) then
print *, 'Syntax: ao_grid <trexio_file> <AO_id> <point_num>'
call exit(-1)
end if
call getarg(1, trexio_filename)
call getarg(2, str_buffer)
read(str_buffer, *) ao_id
call getarg(3, str_buffer)
read(str_buffer, *) point_num_x
if (point_num_x < 0 .or. point_num_x > 300) then
print *, 'Error: 0 < point_num < 300'
call exit(-1)
end if
#+end_src
Create the QMCkl context and initialize it with the wave function
present in the TREXIO file:
#+begin_src f90
qmckl_ctx = qmckl_context_create()
rc = qmckl_trexio_read(qmckl_ctx, trexio_filename, 1_8*len(trim(trexio_filename)))
call qmckl_check_error(rc, 'Read TREXIO')
#+end_src
We need to check that ~ao_id~ is in the range, so we get the total
number of AOs from QMCkl:
#+begin_src f90
rc = qmckl_get_ao_basis_ao_num(qmckl_ctx, ao_num)
call qmckl_check_error(rc, 'Getting ao_num')
if (ao_id < 0 .or. ao_id > ao_num) then
print *, 'Error: 0 < ao_id < ', ao_num
call exit(-1)
end if
#+end_src
Now we will compute the limits of the box in which the molecule fits.
For that, we first need to ask QMCkl the coordinates of nuclei.
#+begin_src f90
rc = qmckl_get_nucleus_num(qmckl_ctx, nucl_num)
call qmckl_check_error(rc, 'Get nucleus num')
allocate( nucl_coord(3, nucl_num) )
rc = qmckl_get_nucleus_coord(qmckl_ctx, 'N', nucl_coord, 3_8*nucl_num)
call qmckl_check_error(rc, 'Get nucleus coord')
#+end_src
We now compute the coordinates of opposite points of the box, and
the distance between points along the 3 directions:
#+begin_src f90
rmin(1) = minval( nucl_coord(1,:) ) - 5.d0
rmin(2) = minval( nucl_coord(2,:) ) - 5.d0
rmin(3) = minval( nucl_coord(3,:) ) - 5.d0
rmax(1) = maxval( nucl_coord(1,:) ) + 5.d0
rmax(2) = maxval( nucl_coord(2,:) ) + 5.d0
rmax(3) = maxval( nucl_coord(3,:) ) + 5.d0
dr(1:3) = (rmax(1:3) - rmin(1:3)) / dble(point_num_x-1)
#+end_src
We now produce the list of point coordinates where the AO will be
evaluated:
#+begin_src f90
point_num = point_num_x**3
allocate( points(point_num, 3) )
ipoint=0
z = rmin(3)
do k=1,point_num_x
y = rmin(2)
do j=1,point_num_x
x = rmin(1)
do i=1,point_num_x
ipoint = ipoint+1
points(ipoint,1) = x
points(ipoint,2) = y
points(ipoint,3) = z
x = x + dr(1)
end do
y = y + dr(2)
end do
z = z + dr(3)
end do
#+end_src
We give the points to QMCkl:
#+begin_src f90
rc = qmckl_set_point(qmckl_ctx, 'T', point_num, points, size(points)*1_8 )
call qmckl_check_error(rc, 'Setting points')
#+end_src
We allocate the space required to retrieve the values, gradients and
Laplacian of all AOs, and ask to retrieve the values of the
AOs computed at the point positions.
#+begin_src f90
allocate( ao_vgl(ao_num, 5, point_num) )
rc = qmckl_get_ao_basis_ao_vgl(qmckl_ctx, ao_vgl, ao_num*5_8*point_num)
call qmckl_check_error(rc, 'Setting points')
#+end_src
We finally print the value and Laplacian of the AO:
#+begin_src f90
do ipoint=1, point_num
print '(3(F10.6,X),2(E20.10,X))', points(ipoint, 1:3), ao_vgl(ao_id,1,ipoint), ao_vgl(ao_id,5,ipoint)
end do
#+end_src
#+begin_src f90
deallocate( nucl_coord, points, ao_vgl )
end program ao_grid
#+end_src

3324
org/qmckl_jastrow.org
View File

File diff suppressed because it is too large Load Diff

99
org/qmckl_local_energy.org
View File

@ -226,7 +226,7 @@ qmckl_exit_code qmckl_get_kinetic_energy(qmckl_context context, double * const k
rc = qmckl_provide_kinetic_energy(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = ctx->electron.walk_num * sizeof(double);
@ -250,7 +250,7 @@ qmckl_exit_code qmckl_provide_kinetic_energy(qmckl_context context) {
return QMCKL_NULL_CONTEXT;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if(!(ctx->nucleus.provided)) {
@ -549,37 +549,28 @@ end function qmckl_compute_kinetic_energy_f
*** Test
#+begin_src c :tangle (eval c_test) :exports none
#define walk_num chbrclf_walk_num
#define elec_num chbrclf_elec_num
#define shell_num chbrclf_shell_num
#define ao_num chbrclf_ao_num
int64_t elec_up_num = chbrclf_elec_up_num;
int64_t elec_dn_num = chbrclf_elec_dn_num;
double* elec_coord = &(chbrclf_elec_coord[0][0][0]);
const int64_t nucl_num = chbrclf_nucl_num;
const double* nucl_charge = chbrclf_charge;
const double* nucl_coord = &(chbrclf_nucl_coord[0][0]);
rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num);
rc = qmckl_set_electron_num (context, chbrclf_elec_up_num, chbrclf_elec_dn_num);
assert (rc == QMCKL_SUCCESS);
rc = qmckl_set_electron_walk_num (context, walk_num);
rc = qmckl_set_electron_walk_num (context, chbrclf_walk_num);
assert (rc == QMCKL_SUCCESS);
assert(qmckl_electron_provided(context));
rc = qmckl_set_electron_coord (context, 'N', elec_coord, walk_num*elec_num*3);
rc = qmckl_set_electron_coord (context, 'N', elec_coord, chbrclf_walk_num*chbrclf_elec_num*3);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_num (context, nucl_num);
rc = qmckl_set_nucleus_num (context, chbrclf_nucl_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), nucl_num*3);
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), chbrclf_nucl_num*3);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_charge(context, nucl_charge, nucl_num);
rc = qmckl_set_nucleus_charge(context, nucl_charge, chbrclf_nucl_num);
assert(rc == QMCKL_SUCCESS);
assert(qmckl_nucleus_provided(context));
@ -611,11 +602,11 @@ rc = qmckl_set_ao_basis_prim_num (context, chbrclf_prim_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_nucleus_index (context, nucleus_index, nucl_num);
rc = qmckl_set_ao_basis_nucleus_index (context, nucleus_index, chbrclf_nucl_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_nucleus_shell_num (context, nucleus_shell_num, nucl_num);
rc = qmckl_set_ao_basis_nucleus_shell_num (context, nucleus_shell_num, chbrclf_nucl_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
@ -655,10 +646,10 @@ assert(rc == QMCKL_SUCCESS);
assert(qmckl_ao_basis_provided(context));
double ao_vgl[walk_num*elec_num][5][chbrclf_ao_num];
double ao_vgl[chbrclf_walk_num*chbrclf_elec_num][5][chbrclf_ao_num];
rc = qmckl_get_ao_basis_ao_vgl(context, &(ao_vgl[0][0][0]),
(int64_t) 5*walk_num*elec_num*chbrclf_ao_num);
(int64_t) 5*chbrclf_walk_num*chbrclf_elec_num*chbrclf_ao_num);
assert (rc == QMCKL_SUCCESS);
/* Set up MO data */
@ -673,31 +664,31 @@ assert (rc == QMCKL_SUCCESS);
assert(qmckl_mo_basis_provided(context));
double mo_vgl[walk_num*elec_num][5][chbrclf_mo_num];
rc = qmckl_get_mo_basis_mo_vgl(context, &(mo_vgl[0][0][0]), 5*walk_num*elec_num*chbrclf_mo_num);
double mo_vgl[chbrclf_walk_num*chbrclf_elec_num][5][chbrclf_mo_num];
rc = qmckl_get_mo_basis_mo_vgl(context, &(mo_vgl[0][0][0]), 5*chbrclf_walk_num*chbrclf_elec_num*chbrclf_mo_num);
assert (rc == QMCKL_SUCCESS);
/* Set up determinant data */
const int64_t det_num_alpha = 1;
const int64_t det_num_beta = 1;
int64_t mo_index_alpha[det_num_alpha][walk_num][elec_up_num];
int64_t mo_index_beta[det_num_alpha][walk_num][elec_dn_num];
#define det_num_alpha 1
#define det_num_beta 1
int64_t mo_index_alpha[det_num_alpha][chbrclf_walk_num][chbrclf_elec_up_num];
int64_t mo_index_beta[det_num_alpha][chbrclf_walk_num][chbrclf_elec_dn_num];
int i, j, k;
for(k = 0; k < det_num_alpha; ++k)
for(i = 0; i < walk_num; ++i)
for(j = 0; j < elec_up_num; ++j)
for(i = 0; i < chbrclf_walk_num; ++i)
for(j = 0; j < chbrclf_elec_up_num; ++j)
mo_index_alpha[k][i][j] = j + 1;
for(k = 0; k < det_num_beta; ++k)
for(i = 0; i < walk_num; ++i)
for(j = 0; j < elec_up_num; ++j)
for(i = 0; i < chbrclf_walk_num; ++i)
for(j = 0; j < chbrclf_elec_up_num; ++j)
mo_index_beta[k][i][j] = j + 1;
rc = qmckl_set_determinant_type (context, typ);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_determinant_walk_num (context, walk_num);
rc = qmckl_set_determinant_walk_num (context, chbrclf_walk_num);
assert (rc == QMCKL_SUCCESS);
rc = qmckl_set_determinant_det_num_alpha (context, det_num_alpha);
@ -714,8 +705,8 @@ assert (rc == QMCKL_SUCCESS);
// Get alpha determinant
double det_vgl_alpha[det_num_alpha][walk_num][5][elec_up_num][elec_up_num];
double det_vgl_beta[det_num_beta][walk_num][5][elec_dn_num][elec_dn_num];
double det_vgl_alpha[det_num_alpha][chbrclf_walk_num][5][chbrclf_elec_up_num][chbrclf_elec_up_num];
double det_vgl_beta[det_num_beta][chbrclf_walk_num][5][chbrclf_elec_dn_num][chbrclf_elec_dn_num];
rc = qmckl_get_det_vgl_alpha(context, &(det_vgl_alpha[0][0][0][0][0]));
assert (rc == QMCKL_SUCCESS);
@ -725,8 +716,8 @@ assert (rc == QMCKL_SUCCESS);
// Get adjoint of the slater-determinant
double det_inv_matrix_alpha[det_num_alpha][walk_num][elec_up_num][elec_up_num];
double det_inv_matrix_beta[det_num_beta][walk_num][elec_dn_num][elec_dn_num];
double det_inv_matrix_alpha[det_num_alpha][chbrclf_walk_num][chbrclf_elec_up_num][chbrclf_elec_up_num];
double det_inv_matrix_beta[det_num_beta][chbrclf_walk_num][chbrclf_elec_dn_num][chbrclf_elec_dn_num];
rc = qmckl_get_det_inv_matrix_alpha(context, &(det_inv_matrix_alpha[0][0][0][0]));
assert (rc == QMCKL_SUCCESS);
@ -736,7 +727,7 @@ assert (rc == QMCKL_SUCCESS);
// Calculate the Kinetic energy
double kinetic_energy[walk_num];
double kinetic_energy[chbrclf_walk_num];
rc = qmckl_get_kinetic_energy(context, &(kinetic_energy[0]));
assert (rc == QMCKL_SUCCESS);
@ -799,7 +790,7 @@ qmckl_exit_code qmckl_get_potential_energy(qmckl_context context, double * const
rc = qmckl_provide_potential_energy(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = ctx->electron.walk_num * sizeof(double);
@ -822,7 +813,7 @@ qmckl_exit_code qmckl_provide_potential_energy(qmckl_context context) {
return QMCKL_NULL_CONTEXT;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
qmckl_exit_code rc;
@ -1034,7 +1025,7 @@ end function qmckl_compute_potential_energy_f
#+begin_src c :tangle (eval c_test) :exports none
// Calculate the Potential energy
double potential_energy[walk_num];
double potential_energy[chbrclf_walk_num];
rc = qmckl_get_potential_energy(context, &(potential_energy[0]));
assert (rc == QMCKL_SUCCESS);
@ -1058,11 +1049,11 @@ E_L = KE + PE
*** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_local_energy(qmckl_context context, double* const local_energy);
qmckl_exit_code qmckl_get_local_energy(qmckl_context context, double* const local_energy, const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_local_energy(qmckl_context context, double * const local_energy) {
qmckl_exit_code qmckl_get_local_energy(qmckl_context context, double * const local_energy, const int64_t size_max) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
@ -1083,11 +1074,17 @@ qmckl_exit_code qmckl_get_local_energy(qmckl_context context, double * const loc
rc = qmckl_provide_local_energy(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = ctx->electron.walk_num * sizeof(double);
memcpy(local_energy, ctx->local_energy.e_local, sze);
const int64_t sze = ctx->electron.walk_num;
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_local_energy",
"input array too small");
}
memcpy(local_energy, ctx->local_energy.e_local, sze * sizeof(double));
return QMCKL_SUCCESS;
}
@ -1106,7 +1103,7 @@ qmckl_exit_code qmckl_provide_local_energy(qmckl_context context) {
return QMCKL_NULL_CONTEXT;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if(!(ctx->nucleus.provided)) {
@ -1290,9 +1287,9 @@ end function qmckl_compute_local_energy_f
#+begin_src c :tangle (eval c_test) :exports none
// Calculate the Local energy
double local_energy[walk_num];
double local_energy[chbrclf_walk_num];
rc = qmckl_get_local_energy(context, &(local_energy[0]));
rc = qmckl_get_local_energy(context, &(local_energy[0]), chbrclf_walk_num);
assert (rc == QMCKL_SUCCESS);
#+end_src
@ -1339,7 +1336,7 @@ qmckl_exit_code qmckl_get_drift_vector(qmckl_context context, double * const dri
rc = qmckl_provide_drift_vector(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = ctx->electron.walk_num * ctx->electron.num * 3 * sizeof(double);
@ -1362,7 +1359,7 @@ qmckl_exit_code qmckl_provide_drift_vector(qmckl_context context) {
return QMCKL_NULL_CONTEXT;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if(!(ctx->nucleus.provided)) {
@ -1645,7 +1642,7 @@ end function qmckl_compute_drift_vector_f
#+begin_src c :tangle (eval c_test) :exports none
// Calculate the Drift vector
double drift_vector[walk_num][elec_num][3];
double drift_vector[chbrclf_walk_num][chbrclf_elec_num][3];
rc = qmckl_get_drift_vector(context, &(drift_vector[0][0][0]));
assert (rc == QMCKL_SUCCESS);

4
org/qmckl_memory.org
View File

@ -116,7 +116,7 @@ void* qmckl_malloc(qmckl_context context, const qmckl_memory_info_struct info) {
assert (qmckl_context_check(context) != QMCKL_NULL_CONTEXT);
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
/* Allocate memory and zero it */
void * pointer = malloc(info.size);
@ -217,7 +217,7 @@ qmckl_exit_code qmckl_free(qmckl_context context, void * const ptr) {
"NULL pointer");
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
qmckl_lock(context);
{

550
org/qmckl_mo.org
View File

@ -92,10 +92,12 @@ int main() {
Computed data:
|---------------+--------------------------+-------------------------------------------------------------------------------------|
| ~mo_vgl~ | ~[point_num][5][mo_num]~ | Value, gradients, Laplacian of the MOs at point positions |
| ~mo_vgl_date~ | ~uint64_t~ | Late modification date of Value, gradients, Laplacian of the MOs at point positions |
|---------------+--------------------------+-------------------------------------------------------------------------------------|
|-----------------+--------------------------+-------------------------------------------------------------------------------------|
| ~mo_value~ | ~[point_num][mo_num]~ | Value of the MOs at point positions |
| ~mo_value_date~ | ~uint64_t~ | Late modification date of the value of the MOs at point positions |
| ~mo_vgl~ | ~[point_num][5][mo_num]~ | Value, gradients, Laplacian of the MOs at point positions |
| ~mo_vgl_date~ | ~uint64_t~ | Late modification date of Value, gradients, Laplacian of the MOs at point positions |
|-----------------+--------------------------+-------------------------------------------------------------------------------------|
** Data structure
@ -106,7 +108,9 @@ typedef struct qmckl_mo_basis_struct {
double * restrict coefficient_t;
double * restrict mo_vgl;
double * restrict mo_value;
uint64_t mo_vgl_date;
uint64_t mo_value_date;
int32_t uninitialized;
bool provided;
@ -131,7 +135,7 @@ qmckl_exit_code qmckl_init_mo_basis(qmckl_context context) {
return false;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
ctx->mo_basis.uninitialized = (1 << 2) - 1;
@ -158,10 +162,9 @@ qmckl_get_mo_basis_mo_num (const qmckl_context context,
QMCKL_INVALID_CONTEXT,
"qmckl_get_mo_basis_mo_num",
NULL);
return (int64_t) 0;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int32_t mask = 1;
@ -200,7 +203,7 @@ qmckl_get_mo_basis_coefficient (const qmckl_context context,
NULL);
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int32_t mask = 1 << 1;
@ -248,7 +251,7 @@ bool qmckl_mo_basis_provided(const qmckl_context context) {
return false;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
return ctx->mo_basis.provided;
@ -257,10 +260,9 @@ bool qmckl_mo_basis_provided(const qmckl_context context) {
#+end_src
*** Fortran interfaces
#+begin_src f90 :tangle (eval fh_func) :comments org
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_get_mo_basis_mo_num (context, &
mo_num) bind(C)
@ -280,7 +282,7 @@ interface
implicit none
integer (c_int64_t) , intent(in) , value :: context
double precision, intent(out) :: coefficient(*)
integer (c_int64_t) , intent(in) , value :: size_max
integer (c_int64_t) , intent(in), value :: size_max
end function qmckl_get_mo_basis_coefficient
end interface
@ -302,7 +304,7 @@ if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
#+end_src
#+NAME:post
@ -383,7 +385,7 @@ qmckl_exit_code qmckl_finalize_mo_basis(qmckl_context context) {
NULL);
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
@ -421,7 +423,464 @@ qmckl_exit_code qmckl_finalize_mo_basis(qmckl_context context) {
* Computation
** Computation of MOs
** Computation of MOs: values only
*** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code
qmckl_get_mo_basis_mo_value(qmckl_context context,
double* const mo_value,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_mo_basis_mo_value(qmckl_context context,
double* const mo_value,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_exit_code rc;
rc = qmckl_provide_ao_value(context);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_mo_value(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
const int64_t sze = ctx->point.num * ctx->mo_basis.mo_num;
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_mo_basis_mo_value",
"input array too small");
}
memcpy(mo_value, ctx->mo_basis.mo_value, sze * sizeof(double));
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_get_mo_basis_mo_value (context, &
mo_value, size_max) bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
double precision, intent(out) :: mo_value(*)
integer (c_int64_t) , intent(in) , value :: size_max
end function qmckl_get_mo_basis_mo_value
end interface
#+end_src
Uses the given array to compute the values.
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code
qmckl_get_mo_basis_mo_value_inplace (qmckl_context context,
double* const mo_value,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_mo_basis_mo_value_inplace (qmckl_context context,
double* const mo_value,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT,
"qmckl_get_mo_basis_mo_value",
NULL);
}
qmckl_exit_code rc;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
const int64_t sze = ctx->mo_basis.mo_num * ctx->point.num;
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_mo_basis_mo_value",
"input array too small");
}
rc = qmckl_context_touch(context);
if (rc != QMCKL_SUCCESS) return rc;
double* old_array = ctx->mo_basis.mo_value;
ctx->mo_basis.mo_value = mo_value;
rc = qmckl_provide_mo_value(context);
if (rc != QMCKL_SUCCESS) return rc;
ctx->mo_basis.mo_value = old_array;
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_get_mo_basis_mo_value_inplace (context, &
mo_value, size_max) bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
double precision, intent(out) :: mo_value(*)
integer (c_int64_t) , intent(in) , value :: size_max
end function qmckl_get_mo_basis_mo_value_inplace
end interface
#+end_src
*** Provide
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_provide_mo_value(qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_provide_mo_value(qmckl_context context)
{
qmckl_exit_code rc;
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if (!ctx->ao_basis.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_ao_basis",
NULL);
}
rc = qmckl_provide_ao_value(context);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_ao_value",
NULL);
}
if (!ctx->mo_basis.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_mo_basis",
NULL);
}
/* Compute if necessary */
if (ctx->point.date > ctx->mo_basis.mo_value_date) {
/* Allocate array */
if (ctx->mo_basis.mo_value == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->point.num * ctx->mo_basis.mo_num * sizeof(double);
double* mo_value = (double*) qmckl_malloc(context, mem_info);
if (mo_value == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_mo_basis_mo_value",
NULL);
}
ctx->mo_basis.mo_value = mo_value;
}
if (ctx->mo_basis.mo_vgl_date == ctx->point.date) {
// mo_vgl has been computed at this step: Just copy the data.
double * v = &(ctx->mo_basis.mo_value[0]);
double * vgl = &(ctx->mo_basis.mo_vgl[0]);
for (int i=0 ; i<ctx->point.num ; ++i) {
for (int k=0 ; k<ctx->mo_basis.mo_num ; ++k) {
v[k] = vgl[k];
}
v += ctx->mo_basis.mo_num;
vgl += ctx->mo_basis.mo_num * 5;
}
} else {
rc = qmckl_compute_mo_basis_mo_value(context,
ctx->ao_basis.ao_num,
ctx->mo_basis.mo_num,
ctx->point.num,
ctx->mo_basis.coefficient_t,
ctx->ao_basis.ao_value,
ctx->mo_basis.mo_value);
if (rc != QMCKL_SUCCESS) {
return rc;
}
}
ctx->mo_basis.mo_value_date = ctx->date;
}
return QMCKL_SUCCESS;
}
#+end_src
*** Compute
:PROPERTIES:
:Name: qmckl_compute_mo_basis_mo_value
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
#+NAME: qmckl_mo_basis_mo_value_args
| Variable | Type | In/Out | Description |
|-----------------+-----------------------------+--------+-------------------------------------------------|
| ~context~ | ~qmckl_context~ | in | Global state |
| ~ao_num~ | ~int64_t~ | in | Number of AOs |
| ~mo_num~ | ~int64_t~ | in | Number of MOs |
| ~point_num~ | ~int64_t~ | in | Number of points |
| ~coefficient_t~ | ~double[mo_num][ao_num]~ | in | Transpose of the AO to MO transformation matrix |
| ~ao_value~ | ~double[point_num][ao_num]~ | in | Value of the AOs |
| ~mo_value~ | ~double[point_num][mo_num]~ | out | Value of the MOs |
The matrix of AO values is very sparse, so we use a sparse-dense
matrix multiplication instead of a dgemm, as exposed in
https://dx.doi.org/10.1007/978-3-642-38718-0_14.
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_mo_basis_mo_value_doc_f(context, &
ao_num, mo_num, point_num, &
coefficient_t, ao_value, mo_value) &
result(info)
use qmckl
implicit none
integer(qmckl_context), intent(in) :: context
integer*8 , intent(in) :: ao_num, mo_num
integer*8 , intent(in) :: point_num
double precision , intent(in) :: ao_value(ao_num,point_num)
double precision , intent(in) :: coefficient_t(mo_num,ao_num)
double precision , intent(out) :: mo_value(mo_num,point_num)
integer*8 :: i,j,k
double precision :: c1, c2, c3, c4, c5
integer*8 :: LDA, LDB, LDC
info = QMCKL_SUCCESS
if (.True.) then ! fast algorithm
do j=1,point_num
mo_value(:,j) = 0.d0
do k=1,ao_num
if (ao_value(k,j) /= 0.d0) then
c1 = ao_value(k,j)
do i=1,mo_num
mo_value(i,j) = mo_value(i,j) + coefficient_t(i,k) * c1
end do
end if
end do
end do
else ! dgemm
LDA = size(coefficient_t,1)
LDB = size(ao_value,1)
LDC = size(mo_value,1)
info = qmckl_dgemm(context,'N', 'N', mo_num, point_num, ao_num, 1.d0, &
coefficient_t, LDA, ao_value, LDB, &
0.d0, mo_value, LDC)
end if
end function qmckl_compute_mo_basis_mo_value_doc_f
#+end_src
#+CALL: generate_c_header(table=qmckl_mo_basis_mo_value_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_value"))
#+RESULTS:
#+begin_src c :tangle (eval h_func) :comments org
qmckl_exit_code qmckl_compute_mo_basis_mo_value (
const qmckl_context context,
const int64_t ao_num,
const int64_t mo_num,
const int64_t point_num,
const double* coefficient_t,
const double* ao_value,
double* const mo_value );
#+end_src
#+CALL: generate_c_header(table=qmckl_mo_basis_mo_value_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_value_doc"))
#+RESULTS:
#+begin_src c :tangle (eval h_func) :comments org
qmckl_exit_code qmckl_compute_mo_basis_mo_value_doc (
const qmckl_context context,
const int64_t ao_num,
const int64_t mo_num,
const int64_t point_num,
const double* coefficient_t,
const double* ao_value,
double* const mo_value );
#+end_src
#+CALL: generate_c_interface(table=qmckl_mo_basis_mo_value_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_value_doc"))
#+RESULTS:
#+begin_src f90 :tangle (eval f) :comments org :exports none
integer(c_int32_t) function qmckl_compute_mo_basis_mo_value_doc &
(context, ao_num, mo_num, point_num, coefficient_t, ao_value, mo_value) &
bind(C) result(info)
use, intrinsic :: iso_c_binding
implicit none
integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(in) , value :: ao_num
integer (c_int64_t) , intent(in) , value :: mo_num
integer (c_int64_t) , intent(in) , value :: point_num
real (c_double ) , intent(in) :: coefficient_t(ao_num,mo_num)
real (c_double ) , intent(in) :: ao_value(ao_num,point_num)
real (c_double ) , intent(out) :: mo_value(mo_num,point_num)
integer(c_int32_t), external :: qmckl_compute_mo_basis_mo_value_doc_f
info = qmckl_compute_mo_basis_mo_value_doc_f &
(context, ao_num, mo_num, point_num, coefficient_t, ao_value, mo_value)
end function qmckl_compute_mo_basis_mo_value_doc
#+end_src
#+begin_src c :tangle (eval c) :comments org
qmckl_exit_code
qmckl_compute_mo_basis_mo_value (const qmckl_context context,
const int64_t ao_num,
const int64_t mo_num,
const int64_t point_num,
const double* coefficient_t,
const double* ao_value,
double* const mo_value )
{
#ifdef HAVE_HPC
return qmckl_compute_mo_basis_mo_value_hpc (context, ao_num, mo_num, point_num, coefficient_t, ao_value, mo_value);
#else
return qmckl_compute_mo_basis_mo_value_doc (context, ao_num, mo_num, point_num, coefficient_t, ao_value, mo_value);
#endif
}
#+end_src
*** HPC version
#+begin_src c :tangle (eval h_func) :comments org
#ifdef HAVE_HPC
qmckl_exit_code
qmckl_compute_mo_basis_mo_value_hpc (const qmckl_context context,
const int64_t ao_num,
const int64_t mo_num,
const int64_t point_num,
const double* coefficient_t,
const double* ao_value,
double* const mo_value );
#endif
#+end_src
#+begin_src c :tangle (eval c) :comments org
#ifdef HAVE_HPC
qmckl_exit_code
qmckl_compute_mo_basis_mo_value_hpc (const qmckl_context context,
const int64_t ao_num,
const int64_t mo_num,
const int64_t point_num,
const double* restrict coefficient_t,
const double* restrict ao_value,
double* restrict const mo_value )
{
assert (context != QMCKL_NULL_CONTEXT);
#ifdef HAVE_OPENMP
#pragma omp parallel for
#endif
for (int64_t ipoint=0 ; ipoint < point_num ; ++ipoint) {
double* restrict const vgl1 = &(mo_value[ipoint*mo_num]);
const double* restrict avgl1 = &(ao_value[ipoint*ao_num]);
for (int64_t i=0 ; i<mo_num ; ++i) {
vgl1[i] = 0.;
}
int64_t nidx=0;
int64_t idx[ao_num];
double av1[ao_num];
for (int64_t k=0 ; k<ao_num ; ++k) {
if (avgl1[k] != 0.) {
idx[nidx] = k;
av1[nidx] = avgl1[k];
++nidx;
}
}
int64_t n;
for (n=0 ; n < nidx-4 ; n+=4) {
const double* restrict ck1 = coefficient_t + idx[n ]*mo_num;
const double* restrict ck2 = coefficient_t + idx[n+1]*mo_num;
const double* restrict ck3 = coefficient_t + idx[n+2]*mo_num;
const double* restrict ck4 = coefficient_t + idx[n+3]*mo_num;
const double a11 = av1[n ];
const double a21 = av1[n+1];
const double a31 = av1[n+2];
const double a41 = av1[n+3];
#ifdef HAVE_OPENMP
#pragma omp simd
#endif
for (int64_t i=0 ; i<mo_num ; ++i) {
vgl1[i] = vgl1[i] + ck1[i] * a11 + ck2[i] * a21 + ck3[i] * a31 + ck4[i] * a41;
}
}
const int64_t n0 = n < 0 ? 0 : n;
for (int64_t m=n0 ; m < nidx ; m+=1) {
const double* restrict ck = coefficient_t + idx[m]*mo_num;
const double a1 = av1[m];
#ifdef HAVE_OPENMP
#pragma omp simd
#endif
for (int64_t i=0 ; i<mo_num ; ++i) {
vgl1[i] += ck[i] * a1;
}
}
}
return QMCKL_SUCCESS;
}
#endif
#+end_src
** Computation of MOs: values, gradient, Laplacian
*** Get
@ -451,7 +910,7 @@ qmckl_get_mo_basis_mo_vgl(qmckl_context context,
rc = qmckl_provide_mo_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
const int64_t sze = ctx->point.num * 5 * ctx->mo_basis.mo_num;
@ -507,7 +966,7 @@ qmckl_get_mo_basis_mo_vgl_inplace (qmckl_context context,
qmckl_exit_code rc;
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
const int64_t sze = ctx->mo_basis.mo_num * 5 * ctx->point.num;
@ -563,7 +1022,7 @@ qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context)
return QMCKL_NULL_CONTEXT;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if (!ctx->ao_basis.provided) {
@ -767,7 +1226,6 @@ qmckl_compute_mo_basis_mo_vgl (const qmckl_context context,
}
#+end_src
*** HPC version
@ -795,6 +1253,8 @@ qmckl_compute_mo_basis_mo_vgl_hpc (const qmckl_context context,
const double* restrict ao_vgl,
double* restrict const mo_vgl )
{
assert (context != QMCKL_NULL_CONTEXT);
#ifdef HAVE_OPENMP
#pragma omp parallel for
#endif
@ -827,7 +1287,6 @@ qmckl_compute_mo_basis_mo_vgl_hpc (const qmckl_context context,
double av4[ao_num];
double av5[ao_num];
for (int64_t k=0 ; k<ao_num ; ++k) {
const double* restrict ck1 = coefficient_t + k*mo_num;
if (avgl1[k] != 0.) {
idx[nidx] = k;
av1[nidx] = avgl1[k];
@ -842,7 +1301,6 @@ qmckl_compute_mo_basis_mo_vgl_hpc (const qmckl_context context,
int64_t n;
for (n=0 ; n < nidx-4 ; n+=4) {
int64_t k = idx[n];
const double* restrict ck1 = coefficient_t + idx[n ]*mo_num;
const double* restrict ck2 = coefficient_t + idx[n+1]*mo_num;
const double* restrict ck3 = coefficient_t + idx[n+2]*mo_num;
@ -886,13 +1344,13 @@ qmckl_compute_mo_basis_mo_vgl_hpc (const qmckl_context context,
}
const int64_t n0 = n < 0 ? 0 : n;
for (int64_t n=n0 ; n < nidx ; n+=1) {
const double* restrict ck = coefficient_t + idx[n]*mo_num;
const double a1 = av1[n];
const double a2 = av2[n];
const double a3 = av3[n];
const double a4 = av4[n];
const double a5 = av5[n];
for (int64_t m=n0 ; m < nidx ; m+=1) {
const double* restrict ck = coefficient_t + idx[m]*mo_num;
const double a1 = av1[m];
const double a2 = av2[m];
const double a3 = av3[m];
const double a4 = av4[m];
const double a5 = av5[m];
#ifdef HAVE_OPENMP
#pragma omp simd
@ -911,9 +1369,9 @@ qmckl_compute_mo_basis_mo_vgl_hpc (const qmckl_context context,
#endif
#+end_src
*** Test
** Test
#+begin_src python :results output :exports none
#+begin_src python :results output :exports none
import numpy as np
def f(a,x,y):
@ -973,9 +1431,9 @@ print ( "[2][1][15][14] : %25.15e"% df(a,x,y,2))
print ( "[3][1][15][14] : %25.15e"% df(a,x,y,3))
print ( "[4][1][15][14] : %25.15e"% lf(a,x,y))
#+end_src
#+end_src
#+begin_src c :tangle (eval c_test) :exports none
#+begin_src c :tangle (eval c_test) :exports none
{
#define walk_num chbrclf_walk_num
#define elec_num chbrclf_elec_num
@ -1100,10 +1558,36 @@ assert (rc == QMCKL_SUCCESS);
assert(qmckl_mo_basis_provided(context));
rc = qmckl_context_touch(context);
assert (rc == QMCKL_SUCCESS);
double mo_value[walk_num*elec_num][chbrclf_mo_num];
rc = qmckl_get_mo_basis_mo_value(context, &(mo_value[0][0]), walk_num*elec_num*chbrclf_mo_num);
assert (rc == QMCKL_SUCCESS);
double mo_vgl[walk_num*elec_num][5][chbrclf_mo_num];
rc = qmckl_get_mo_basis_mo_vgl(context, &(mo_vgl[0][0][0]), walk_num*elec_num*5*chbrclf_mo_num);
assert (rc == QMCKL_SUCCESS);
for (int i=0 ; i< walk_num*elec_num; ++i) {
for (int k=0 ; k< chbrclf_mo_num ; ++k) {
assert(fabs(mo_vgl[i][0][k] - mo_value[i][k]) < 1.e-12) ;
}
}
rc = qmckl_context_touch(context);
assert (rc == QMCKL_SUCCESS);
rc = qmckl_get_mo_basis_mo_value(context, &(mo_value[0][0]), walk_num*elec_num*chbrclf_mo_num);
assert (rc == QMCKL_SUCCESS);
for (int i=0 ; i< walk_num*elec_num; ++i) {
for (int k=0 ; k< chbrclf_mo_num ; ++k) {
assert(fabs(mo_vgl[i][0][k] - mo_value[i][k]) < 1.e-12) ;
}
}
// Test overlap of MO
//double point_x[10];
//double point_y[10];
@ -1161,7 +1645,7 @@ printf(" mo_vgl mo_vgl[1][26][224] %25.15e\n", mo_vgl[2][1][3]);
printf("\n");
}
#+end_src
#+end_src
* End of files :noexport:

63
org/qmckl_nucleus.org
View File

@ -125,7 +125,7 @@ qmckl_exit_code qmckl_init_nucleus(qmckl_context context) {
return false;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
ctx->nucleus.uninitialized = (1 << 3) - 1;
@ -167,7 +167,7 @@ qmckl_get_nucleus_num (const qmckl_context context, int64_t* const num) {
"num is a null pointer");
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int32_t mask = 1 << 0;
@ -226,7 +226,7 @@ qmckl_get_nucleus_charge (const qmckl_context context,
"charge is a null pointer");
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int32_t mask = 1 << 1;
@ -293,7 +293,7 @@ qmckl_get_nucleus_rescale_factor (const qmckl_context context,
"rescale_factor_kappa is a null pointer");
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
assert (ctx->nucleus.rescale_factor_kappa > 0.0);
@ -351,7 +351,7 @@ qmckl_get_nucleus_coord (const qmckl_context context,
"coord is a null pointer");
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int32_t mask = 1 << 2;
@ -410,7 +410,7 @@ bool qmckl_nucleus_provided(const qmckl_context context) {
return false;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
return ctx->nucleus.provided;
@ -425,7 +425,7 @@ if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
#+end_src
#+NAME:post2
@ -672,7 +672,6 @@ end interface
** Test
#+begin_src c :tangle (eval c_test)
const int64_t nucl_num = chbrclf_nucl_num;
const double* nucl_charge = chbrclf_charge;
const double* nucl_coord = &(chbrclf_nucl_coord[0][0]);
const double nucl_rescale_factor_kappa = 2.0;
@ -688,13 +687,13 @@ rc = qmckl_get_nucleus_num (context, &n);
assert(rc == QMCKL_NOT_PROVIDED);
rc = qmckl_set_nucleus_num (context, nucl_num);
rc = qmckl_set_nucleus_num (context, chbrclf_nucl_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_nucleus_provided(context));
rc = qmckl_get_nucleus_num (context, &n);
assert(rc == QMCKL_SUCCESS);
assert(n == nucl_num);
assert(n == chbrclf_nucl_num);
double k;
rc = qmckl_get_nucleus_rescale_factor (context, &k);
@ -709,41 +708,41 @@ rc = qmckl_get_nucleus_rescale_factor (context, &k);
assert(rc == QMCKL_SUCCESS);
assert(k == nucl_rescale_factor_kappa);
double nucl_coord2[3*nucl_num];
double nucl_coord2[3*chbrclf_nucl_num];
rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2, 3*nucl_num);
rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2, 3*chbrclf_nucl_num);
assert(rc == QMCKL_NOT_PROVIDED);
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), 3*nucl_num);
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), 3*chbrclf_nucl_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_nucleus_provided(context));
rc = qmckl_get_nucleus_coord (context, 'N', nucl_coord2, 3*nucl_num);
rc = qmckl_get_nucleus_coord (context, 'N', nucl_coord2, 3*chbrclf_nucl_num);
assert(rc == QMCKL_SUCCESS);
for (size_t k=0 ; k<3 ; ++k) {
for (int64_t i=0 ; i<nucl_num ; ++i) {
assert( nucl_coord[nucl_num*k+i] == nucl_coord2[3*i+k] );
for (int64_t i=0 ; i<chbrclf_nucl_num ; ++i) {
assert( nucl_coord[chbrclf_nucl_num*k+i] == nucl_coord2[3*i+k] );
}
}
rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2, 3*nucl_num);
rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2, 3*chbrclf_nucl_num);
assert(rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i<3*nucl_num ; ++i) {
for (int64_t i=0 ; i<3*chbrclf_nucl_num ; ++i) {
assert( nucl_coord[i] == nucl_coord2[i] );
}
double nucl_charge2[nucl_num];
double nucl_charge2[chbrclf_nucl_num];
rc = qmckl_get_nucleus_charge(context, nucl_charge2, nucl_num);
rc = qmckl_get_nucleus_charge(context, nucl_charge2, chbrclf_nucl_num);
assert(rc == QMCKL_NOT_PROVIDED);
rc = qmckl_set_nucleus_charge(context, nucl_charge, nucl_num);
rc = qmckl_set_nucleus_charge(context, nucl_charge, chbrclf_nucl_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_get_nucleus_charge(context, nucl_charge2, nucl_num);
rc = qmckl_get_nucleus_charge(context, nucl_charge2, chbrclf_nucl_num);
assert(rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i<nucl_num ; ++i) {
for (int64_t i=0 ; i<chbrclf_nucl_num ; ++i) {
assert( nucl_charge[i] == nucl_charge2[i] );
}
assert(qmckl_nucleus_provided(context));
@ -784,7 +783,7 @@ qmckl_get_nucleus_nn_distance(qmckl_context context,
qmckl_exit_code rc = qmckl_provide_nn_distance(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
const int64_t sze = ctx->nucleus.num * ctx->nucleus.num;
@ -828,7 +827,7 @@ qmckl_exit_code qmckl_provide_nn_distance(qmckl_context context)
return (char) 0;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if (!ctx->nucleus.provided) return QMCKL_NOT_PROVIDED;
@ -940,10 +939,10 @@ qmckl_exit_code qmckl_compute_nn_distance (
assert(qmckl_nucleus_provided(context));
double distance[nucl_num*nucl_num];
rc = qmckl_get_nucleus_nn_distance(context, distance, nucl_num*nucl_num);
double distance[chbrclf_nucl_num*chbrclf_nucl_num];
rc = qmckl_get_nucleus_nn_distance(context, distance, chbrclf_nucl_num*chbrclf_nucl_num);
assert(distance[0] == 0.);
assert(distance[1] == distance[nucl_num]);
assert(distance[1] == distance[chbrclf_nucl_num]);
assert(fabs(distance[1]-2.070304721365169) < 1.e-12);
#+end_src
@ -973,7 +972,7 @@ qmckl_get_nucleus_nn_distance_rescaled(qmckl_context context,
qmckl_exit_code rc = qmckl_provide_nn_distance_rescaled(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
const int64_t sze = ctx->nucleus.num * ctx->nucleus.num;
@ -1019,7 +1018,7 @@ qmckl_exit_code qmckl_provide_nn_distance_rescaled(qmckl_context context)
return (char) 0;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if (!ctx->nucleus.provided) return QMCKL_NOT_PROVIDED;
@ -1167,7 +1166,7 @@ qmckl_exit_code qmckl_get_nucleus_repulsion(qmckl_context context, double* const
qmckl_exit_code rc = qmckl_provide_nucleus_repulsion(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
*energy = ctx->nucleus.repulsion;
@ -1203,7 +1202,7 @@ qmckl_exit_code qmckl_provide_nucleus_repulsion(qmckl_context context)
return (char) 0;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
qmckl_exit_code rc;

8
org/qmckl_numprec.org
View File

@ -141,7 +141,7 @@ qmckl_exit_code qmckl_set_numprec_precision(const qmckl_context context, const i
"precision > 53");
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
/* This should be always true because the context is valid */
assert (ctx != NULL);
@ -185,7 +185,7 @@ int qmckl_get_numprec_precision(const qmckl_context context) {
"");
}
const qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
const qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
return ctx->numprec.precision;
}
#+end_src
@ -232,7 +232,7 @@ qmckl_exit_code qmckl_set_numprec_range(const qmckl_context context, const int r
"range > 11");
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
/* This should be always true because the context is valid */
assert (ctx != NULL);
@ -275,7 +275,7 @@ int qmckl_get_numprec_range(const qmckl_context context) {
"");
}
const qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
const qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
return ctx->numprec.range;
}
#+end_src

32
org/qmckl_point.org
View File

@ -81,7 +81,7 @@ int main() {
|----------+----------------+-------------------------------------------|
| ~num~ | ~int64_t~ | Total number of points |
| ~date~ | ~uint64_t~ | Last modification date of the coordinates |
| ~coord~ | ~qmckl_matrix~ | ~num~ \times 3 matrix3 |
| ~coord~ | ~qmckl_matrix~ | ~num~ \times 3 matrix |
We consider that the matrix is stored 'transposed' and 'normal'
corresponds to the 3 \times ~num~ matrix.
@ -108,7 +108,7 @@ qmckl_exit_code qmckl_init_point(qmckl_context context) {
return false;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
memset(&(ctx->point), 0, sizeof(qmckl_point_struct));
@ -148,7 +148,7 @@ qmckl_get_point_num (const qmckl_context context, int64_t* const num) {
"num is a null pointer");
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
assert (ctx->point.num >= (int64_t) 0);
@ -202,7 +202,7 @@ qmckl_get_point(const qmckl_context context,
"coord is a null pointer");
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int64_t point_num = ctx->point.num;
@ -263,8 +263,9 @@ end interface
#+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code qmckl_set_point (qmckl_context context,
const char transp,
const int64_t num,
const double* coord,
const int64_t num);
const int64_t size_max);
#+end_src
Copy a sequence of ~num~ points $(x,y,z)$ into the context.
@ -273,14 +274,22 @@ qmckl_exit_code qmckl_set_point (qmckl_context context,
qmckl_exit_code
qmckl_set_point (qmckl_context context,
const char transp,
const int64_t num,
const double* coord,
const int64_t num)
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
if (size_max < 3*num) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_4,
"qmckl_set_point",
"Array too small");
}
if (transp != 'N' && transp != 'T') {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
@ -295,7 +304,7 @@ qmckl_set_point (qmckl_context context,
"coord is a NULL pointer");
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
qmckl_exit_code rc;
@ -349,15 +358,16 @@ qmckl_set_point (qmckl_context context,
#+begin_src f90 :comments org :tangle (eval fh_func) :noweb yes
interface
integer(c_int32_t) function qmckl_set_point(context, &
transp, coord, num) bind(C)
transp, num, 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 :: num
real (c_double ) , intent(in) :: coord(*)
integer (c_int64_t) , intent(in) , value :: size_max
end function
end interface
#+end_src
@ -380,7 +390,7 @@ 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);
rc = qmckl_set_point (context, 'N', point_num, coord, (point_num*3));
assert(rc == QMCKL_SUCCESS);
int64_t n;
@ -404,7 +414,7 @@ for (int64_t i=0 ; i<point_num ; ++i) {
assert( coord[3*i+2] == coord2[i+point_num*2] );
}
rc = qmckl_set_point (context, 'T', coord2, point_num);
rc = qmckl_set_point (context, 'T', point_num, coord2, (point_num*3));
assert(rc == QMCKL_SUCCESS);
rc = qmckl_get_point (context, 'N', coord3, (point_num*3));

2
org/qmckl_sherman_morrison_woodbury.org
View File

@ -965,7 +965,7 @@ qmckl_exit_code qmckl_sherman_morrison_smw32s(const qmckl_context context,
rc = qmckl_woodbury_3(context, LDS, Dim, Updates_3block, Updates_index_3block, breakdown, Slater_inv, determinant);
if (rc != 0) { // Send the entire block to slagel_splitting
uint64_t l = 0;
rc = qmckl_slagel_splitting(LDS, Dim, 3, Updates_3block, Updates_index_3block,
(void) qmckl_slagel_splitting(LDS, Dim, 3, Updates_3block, Updates_index_3block,
breakdown, Slater_inv, later_updates + (Dim * later), later_index + later, &l, determinant);
later = later + l;
}

1
org/table_of_contents
View File

@ -18,3 +18,4 @@ qmckl_utils.org
qmckl_trexio.org
qmckl_tests.org
qmckl_verificarlo.org
qmckl_examples.org

42
python/README.md Normal file
View File

@ -0,0 +1,42 @@
# Python API of the QMCkl library
## Requirements
- `setuptools`
- `numpy`
- `swig` (>= 4.0)
## Manual installation
1. Install the QMCkl library (see upstream instructions)
2. `./manual_install_qmckl.sh` which should do the following
3. Copy the produced `_qmckl.so` and `qmckl.py` files into your working directory and do not forget to `import qmckl` in your Python scripts
The second step executes the following under the hood:
1. `./build_qmckl.sh`
2. `<c-compiler> -I/usr/include/python3.8 -c -fPIC qmckl_wrap.c` to compile the wrapper code into an object file using the `<c-compiler>` (replace with your C compiler, e.g. `gcc`) on your machine
3. `<c-compiler> -shared qmckl_wrap.o -lqmckl -o _qmckl.so` to produce the final C extension (this requires the `qmckl` library to be installed and present in the linking paths together with all its dependencies like `trexio`)
## Python-ic installation (recommended)
1. Install the QMCkl library (see upstream instructions)
2. `./pip_install_qmckl.sh`
The last step runs `./build_qmckl.sh`, copies the result into the `qmckl/` directory and
then runs `pip install .` to install the `qmckl` Python package in your environment.
## SWIG pre-processing
Both aforementioned steps call `build_qmckl.sh` script which does the following pre-processing for SWIG
1. Copy the latest `qmckl.h` file fron `include/` into the `src/` directory
2. `python process_header.py` to generate `qmckl_include.i` list of SWIG patterns
3. `swig -python -py3 -builtin -threads -o qmckl_wrap.c qmckl.i` to generate the SWIG wrapper code in C and `qmckl.py` module in Python.
**Note:** for this to work three files have to be present in the working directory: `qmckl.i`, `qmckl_include.i` and `numpy.i`.

25
python/manual_install_qmckl.sh Executable file
View File

@ -0,0 +1,25 @@
#!/bin/bash
set -x
set -e
# swig pre-processing
./build_qmckl.sh
cd src/
# compile the wrapper code
cc -c -fPIC `pkg-config --cflags qmckl` -I/usr/include/python3.8 qmckl_wrap.c -o qmckl_wrap.o
# link against the previously installed QMCkl library (as detected by pkg-config)
cc -shared qmckl_wrap.o `pkg-config --libs qmckl` -o _qmckl.so
cd ..
# copy the produced files into the test dir
cp src/_qmckl.so src/qmckl.py test/
# run tests
cd test/
python test_api.py
cd ..

7
python/pyproject.toml Normal file
View File

@ -0,0 +1,7 @@
[build-system]
requires = [
"setuptools>=42",
"wheel",
"numpy>=1.17.3"
]
build-backend = "setuptools.build_meta"

2
python/requirements.txt Normal file
View File

@ -0,0 +1,2 @@
setuptools>=42
numpy>=1.17.3

75
python/setup.py Normal file
View File

@ -0,0 +1,75 @@
#!/usr/bin/env python3
"""
setup.py file for qmckl package
"""
import os, sys
from setuptools import setup, Extension
from os.path import join
# Read the long description
with open("README.md", "r") as fh:
long_description = fh.read()
# this was recommended to solve the problem of the missing numpy header files
try:
import numpy
except ImportError:
raise Exception("numpy Python package cannot be imported.")
numpy_includedir = numpy.get_include()
# Define the name of the Python package
MODULE_NAME = "qmckl"
# derive the QMCkl libdir and includedir
QMCKL_LIBDIR = os.environ.get("QMCKL_LIBDIR", None)
QMCKL_INCLUDEDIR = os.environ.get("QMCKL_INCLUDEDIR", None)
libdir_undefined = QMCKL_LIBDIR is None or QMCKL_LIBDIR==""
includedir_undefined = QMCKL_INCLUDEDIR is None or QMCKL_INCLUDEDIR==""
# Define qmckl extension module based on SWIG interface file (requires qmckl.h)
qmckl_module = Extension(name = "._" + MODULE_NAME,
sources = [ join("src", MODULE_NAME + "_wrap.c") ],
include_dirs = [numpy_includedir, QMCKL_INCLUDEDIR],
#library_dirs = [QMCKL_LIBDIR],
runtime_library_dirs = [QMCKL_LIBDIR],
libraries = ["qmckl"],
extra_compile_args = ["-Wall"],
extra_link_args = ["-L" + QMCKL_LIBDIR],
depends = [ join("src", "qmckl.h") ],
language = "c"
)
setup(name = MODULE_NAME,
version = "0.2.0",
author = "TREX-CoE",
author_email = "posenitskiy@irsamc.ups-tlse.fr",
description = """Python API of the QMCkl library""",
long_description = long_description,
long_description_content_type = "text/markdown",
ext_modules = [qmckl_module],
py_modules = [MODULE_NAME],
url = "https://github.com/TREX-CoE/qmckl",
license = "BSD",
classifiers=[
"Intended Audience :: Science/Research",
"Intended Audience :: Developers",
"Topic :: Scientific/Engineering",
"Programming Language :: C",
"Programming Language :: Python",
"Programming Language :: Python :: 3",
"Programming Language :: Python :: 3 :: Only",
"Programming Language :: Python :: Implementation :: CPython",
"License :: OSI Approved :: BSD License",
"Operating System :: POSIX",
"Operating System :: Unix",
"Operating System :: MacOS"
],
python_requires = ">=3.0",
install_requires = ["numpy>=1.17.3"]
)

3183
python/src/numpy.i Normal file
View File

File diff suppressed because it is too large Load Diff

188
python/src/process_header.py Normal file
View File

@ -0,0 +1,188 @@
#!/usr/bin/env python3
import os
import sys
qmckl_h = sys.argv[1]
collect = False
process = False
get_name = False
block = []
res_str = ''
func_name = ''
arrays = {}
numbers = {}
qmckl_public_api = []
qmckl_errors = []
with open(qmckl_h, 'r') as f_in:
for line in f_in:
# get the errors but without the type cast because SWIG does not recognize it
if '#define' in line and 'qmckl_exit_code' in line:
qmckl_errors.append(line.strip().replace('(qmckl_exit_code)',''))
continue
if get_name:
words = line.strip().split()
if '(' in words[0]:
func_name = words[0].split('(')[0]
else:
func_name = words[0]
if 'get' in func_name or 'set' in func_name:
qmckl_public_api.append(func_name)
get_name = False
if 'qmckl_exit_code' in line:
words = line.strip().split()
if len(words) > 1 and 'qmckl_exit_code' in words[0]:
# this means that the function name is on the same line as `qmckl_exit_code`
func_name = words[1].split('(')[0]
if 'get' in func_name or 'set' in func_name:
qmckl_public_api.append(func_name)
elif len(words) == 1:
# this means that the function name is the first element on the next line
get_name = True
#continue # do not `continue` here otherwise collect is not True for some functions
# process functions - oneliners (for arrays)
if 'size_max' in line and ';' in line:
tmp_list = line.split(',')
for i,s in enumerate(tmp_list):
if 'size_max' in s:
end_str = tmp_list[i].replace(';','').replace('\n','')
pattern = f"({tmp_list[i-1]} ,{end_str}"
datatype = tmp_list[i-1].replace('const','').replace('*','').split()[0]
arrays[func_name] = {
'datatype' : datatype,
'pattern' : pattern
}
#if 'qmckl_get_jastrow_type_nucl_vector' in func_name:
# print(line)
# print(pattern)
continue
# if size_max is not provided then the function should deal with numbers or string
#elif 'num' in line and 'get' in func_name:
elif ';' in line and 'get' in func_name:
# special case
if 'size_max' in line:
continue
#print(line)
tmp_str = line.split(',')[-1].strip()
pattern = tmp_str.replace(')','').replace(';','')
datatype = pattern.replace('const','').replace('*','').split()[0]
numbers[func_name] = {
'datatype' : datatype,
'pattern' : pattern
}
continue
# for multilne functions - append line by line to the list
else:
block.append(line)
collect = True
continue
# if size_max is encountered within the multiline function
if 'size_max' in line and collect:
#if 'qmckl_get_electron_rescale_factor_en' in func_name:
# print("LOL")
# this will not work for 2-line functions where array argument is on the same line as
# func name and size_max argument is on the next line
if not 'qmckl_exit_code' in block[-1] and not '*/' in line:
pattern = '(' + block[-1].strip() + line.strip().replace(';','')
datatype = pattern.replace('const','').replace('*','').replace('(','').split()[0]
collect = False
block = []
arrays[func_name] = {
'datatype' : datatype,
'pattern' : pattern
}
continue
#if 'num' in line and 'get' in func_name and not 'qmckl_get' in line and collect:
if 'get' in func_name and not 'qmckl_get' in line and collect and ';' in line:
#print(func_name)
#print(line)
pattern = line.replace(';','').replace(')','').strip()
datatype = pattern.replace('const','').replace('*','').split()[0]
collect = False
block = []
numbers[func_name] = {
'datatype' : datatype,
'pattern' : pattern
}
continue
# stop/continue multiline function analyzer
if collect and ')' in line:
collect = False
block = []
continue
else:
block.append(line)
continue
# remove buggy qmckl_get_electron_rescale_factor_en key
#arrays.pop('qmckl_get_electron_rescale_factor_en')
processed = list(arrays.keys()) + list(numbers.keys())
#for pub_func in qmckl_public_api:
#if pub_func not in processed and 'set' not in pub_func:
#print("TODO", pub_func)
#print(v['datatype'])
#for k,v in numbers.items():
# print(v)
with open("qmckl_include.i", 'w') as f_out:
# write the list of errors as constants without the type cast
for e in qmckl_errors:
line = e.replace('#define', '%constant qmckl_exit_code').replace('(','=').replace(')',';')
f_out.write(line + '\n')
swig_type = ''
for v in numbers.values():
if 'int' in v['datatype']:
swig_type = 'int'
elif 'float' in v['datatype'] or 'double' in v['datatype']:
swig_type = 'float'
elif 'char' in v['datatype'] or 'bool' in v['datatype']:
#print('SWIG, skipping', v['datatype'], v['pattern'])
continue
else:
raise TypeError(f"Unknown datatype for swig conversion: {v['datatype']}")
f_out.write(f"%apply {swig_type} *OUTPUT {{ {v['pattern']} }};\n")
for k,v in arrays.items():
if 'char' in v['datatype']:
#print("String type", k, v)
pass
if len(v['pattern'].split(',')) != 2:
print('Problemo', k, v)
continue
if 'get' in k:
f_out.write(f"%apply ( {v['datatype']}* ARGOUT_ARRAY1 , int64_t DIM1 ) {{ {v['pattern']} }};\n")
elif 'set' in k:
f_out.write(f"%apply ( {v['datatype']}* IN_ARRAY1 , int64_t DIM1 ) {{ {v['pattern']} }};\n")
#else:
#print("HOW-TO ?", k)

88
python/src/qmckl.i Normal file
View File

@ -0,0 +1,88 @@
%module qmckl
/* Define SWIGWORDSIZE in order to properly align long integers on 64-bit system */
#define SWIGWORDSIZE64
%{
#define SWIG_FILE_WITH_INIT
/* Include the headers in the wrapper code */
#include "qmckl.h"
%}
/*
* Get rid of the function prefixes, as the scripting language will use
* the module's namespace.
*/
%rename("%(strip:[qmckl_])s") "";
/* Include stdint to recognize types from stdint.h */
%include <stdint.i>
/* Include typemaps to play with input/output re-casting (e.g. C pointers) */
%include typemaps.i
%apply int *OUTPUT { qmckl_exit_code *exit_code};
/* Avoid passing file_name length as an additiona argument */
%apply (char *STRING, int LENGTH) { (const char* file_name, const int64_t size_max) };
/* For functions that return strings */
%include <cstring.i>
%cstring_bounded_output(char* function_name, 1024);
%cstring_bounded_output(char* message, 1024);
/* This block is needed make SWIG convert NumPy arrays to/from from the C pointer and size_max argument.
NOTE: `numpy.i` interface file is not part of SWIG but it is included in the numpy distribution (under numpy/tools/swig/numpy.i)
*/
%include "numpy.i"
%init %{
import_array();
%}
/* Typemaps below change the type of numpy array dimensions from int to int64_t */
%numpy_typemaps(double, NPY_DOUBLE, int64_t)
%numpy_typemaps(float, NPY_FLOAT, int64_t)
%numpy_typemaps(int32_t, NPY_INT32, int64_t)
%numpy_typemaps(int64_t, NPY_INT64, int64_t)
/* Include typemaps generated by the process_header.py script */
%include "qmckl_include.i"
/* Handle properly get_point */
/* exception.i is a generic (language-independent) module */
%include "exception.i"
/* Error handling */
%typemap(out) qmckl_exit_code %{
if ($1 != QMCKL_SUCCESS) {
SWIG_exception(SWIG_RuntimeError, qmckl_string_of_error($1));
}
$result = Py_None;
Py_INCREF(Py_None); /* Py_None is a singleton so increment its reference if used. */
%}
/* More swig-y solution (e.g. compatible beyond Python) BUT it does not consume the qmckl_exit_code output as the solution above
TODO: the sizeof() check below if a dummy workaround
It is good to skip exception raise for functions like context_create and others, but might fail
if sizeof(result) == sizeof(qmckl_exit_code), e.g. for functions that return non-zero integers or floats
*/
/*
%exception {
$action
if (result != QMCKL_SUCCESS && sizeof(result) == sizeof(qmckl_exit_code)) {
SWIG_exception_fail(SWIG_RuntimeError, qmckl_string_of_error(result));
}
}
*/
/* The exception handling above does not work for void functions like lock/unlock so exclude them for now */
/*
%ignore qmckl_lock;
%ignore qmckl_unlock;
*/
/* Parse the header files to generate wrappers */
%include "qmckl.h"

BIN
python/test/data/Alz_small.h5 Normal file
View File

Binary file not shown.

47404
python/test/data/data.py Normal file
View File

File diff suppressed because it is too large Load Diff

51
python/test/test_api.py Normal file
View File

@ -0,0 +1,51 @@
"""
This is the test of the Python API of the QMCkl library.
It is the `bench_mos.c` C code adapted from the `bench`
repo and translated into Python with some modifications.
"""
from os.path import join
import time
import qmckl as pq
from data.data import coord
walk_num = 100
elec_num = 158
ITERMAX = 10
ctx = pq.context_create()
try:
pq.trexio_read(ctx, 'fake.h5')
except RuntimeError:
print('Error handling check: passed')
fname = join('data', 'Alz_small.h5')
pq.trexio_read(ctx, fname)
print('trexio_read: passed')
pq.set_electron_walk_num(ctx, walk_num)
mo_num = pq.get_mo_basis_mo_num(ctx)
assert mo_num == 404
pq.set_electron_coord(ctx, 'T', coord)
size_max = 5*walk_num*elec_num*mo_num
mo_vgl = pq.get_mo_basis_mo_vgl(ctx, size_max)
assert mo_vgl.size == size_max
start = time.clock_gettime_ns(time.CLOCK_REALTIME)
for _ in range(ITERMAX):
mo_vgl_in = pq.get_mo_basis_mo_vgl_inplace(ctx, size_max)
end = time.clock_gettime_ns(time.CLOCK_REALTIME)
print(f'Time for the calculation of 1 step : {(end-start)*.000001/ITERMAX} ms')