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

Added AO basis in context. Tests to do

Browse Source
This commit is contained in:
Anthony Scemama 2020-11-14 18:27:38 +01:00
parent 8a53306a63
commit 6b797bd5d4
3 changed files with 494 additions and 180 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

184
src/README.org
View File

@ -1,58 +1,53 @@
#+TITLE: QMCkl source code documentation
#+EXPORT_FILE_NAME: index.html
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="http://www.pirilampo.org/styles/readtheorg/css/htmlize.css"/>
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="http://www.pirilampo.org/styles/readtheorg/css/readtheorg.css"/>
#+HTML_HEAD: <script src="https://ajax.googleapis.com/ajax/libs/jquery/2.1.3/jquery.min.js"></script>
#+HTML_HEAD: <script src="https://maxcdn.bootstrapcdn.com/bootstrap/3.3.4/js/bootstrap.min.js"></script>
#+HTML_HEAD: <script type="text/javascript" src="http://www.pirilampo.org/styles/lib/js/jquery.stickytableheaders.js"></script>
#+HTML_HEAD: <script type="text/javascript" src="http://www.pirilampo.org/styles/readtheorg/js/readtheorg.js"></script>
#+SETUPFILE: https://fniessen.github.io/org-html-themes/setup/theme-readtheorg.setup
* Introduction
The ultimate goal of QMCkl is to provide a high-performance
implementation of the main kernels of QMC. In this particular
repository, we focus on the definition of the API and the tests,
and on a /pedagogical/ presentation of the algorithms. We expect the
HPC experts to use this repository as a reference for re-writing
The ultimate goal of QMCkl is to provide a high-performance
implementation of the main kernels of QMC. In this particular
repository, we focus on the definition of the API and the tests, and
on a /pedagogical/ presentation of the algorithms. We expect the
HPC experts to use this repository as a reference for re-writing
optimized libraries.
Literate programming is particularly adapted in this context.
Source files are written in [[https://karl-voit.at/2017/09/23/orgmode-as-markup-only/][org-mode]] format, to provide useful
Literate programming is particularly adapted in this context.
Source files are written in [[https://karl-voit.at/2017/09/23/orgmode-as-markup-only/][org-mode]] format, to provide useful
comments and LaTex formulas close to the code. There exists multiple
possibilities to convert org-mode files into different formats such as
HTML or pdf.
For a tutorial on literate programming with org-mode, follow
[[http://www.howardism.org/Technical/Emacs/literate-programming-tutorial.html][this link]].
possibilities to convert org-mode files into different formats such
as HTML or pdf. For a tutorial on literate programming with
org-mode, follow [[http://www.howardism.org/Technical/Emacs/literate-programming-tutorial.html][this link]].
The code is extracted from the org files using Emacs as a command-line
tool in the =Makefile=, and then the produced files are compiled.
The code is extracted from the org files using Emacs as a
command-line tool in the =Makefile=, and then the produced files are
compiled.
** Language used
Fortran is one of the most common languages used by the community,
and is simple enough to make the algorithms readable. Hence we
propose in this pedagogical implementation of QMCkl to use Fortran
to express the algorithms. For specific internal functions where
Fortran is one of the most common languages used by the community,
and is simple enough to make the algorithms readable. Hence we
propose in this pedagogical implementation of QMCkl to use Fortran
to express the algorithms. For specific internal functions where
the C language is more natural, C is used.
As Fortran modules generate compiler-dependent files, the use of
As Fortran modules generate compiler-dependent files, the use of
modules is restricted to the internal use of the library, otherwise
the compliance with C is violated.
The external dependencies should be kept as small as possible, so
external libraries should be used /only/ if their used is strongly
The external dependencies should be kept as small as possible, so
external libraries should be used /only/ if their used is strongly
justified.
** Source code editing
Any text editor can be used to edit org-mode files. For a better
user experience Emacs is recommended.
For users hating Emacs, it is good to know that Emacs can behave
like Vim when switched into ``Evil'' mode. There also exists
[[https://www.spacemacs.org][Spacemacs]] which helps the transition for Vim users.
Any text editor can be used to edit org-mode files. For a better
user experience Emacs is recommended. For users hating Emacs, it
is good to know that Emacs can behave like Vim when switched into
``Evil'' mode. There also exists [[https://www.spacemacs.org][Spacemacs]] which helps the
transition for Vim users.
For users with a preference for Jupyter notebooks, the following
For users with a preference for Jupyter notebooks, the following
script can convert jupyter notebooks to org-mode files:
#+BEGIN_SRC sh tangle: nb_to_org.sh
@ -72,24 +67,25 @@ rm ${nb}.md
** Writing in Fortran
The Fortran source files should provide a C interface using
=iso_c_binding=. The name of the Fortran source files should end
with =_f.f90= to be properly handled by the Makefile.
The names of the functions defined in fortran should be the same as
those exposed in the API suffixed by =_f=.
Fortran interface files should also be written in a file with a
=.fh= extension.
The Fortran source files should provide a C interface using
=iso_c_binding=. The name of the Fortran source files should end
with =_f.f90= to be properly handled by the Makefile. The names of
the functions defined in fortran should be the same as those
exposed in the API suffixed by =_f=. Fortran interface files
should also be written in the =qmckl_f.f90= file.
For more guidelines on using Fortran to generate a C interface, see
[[http://fortranwiki.org/fortran/show/Generating+C+Interfaces][this link]]
[[http://fortranwiki.org/fortran/show/Generating+C+Interfaces][this link]].
** Coding style
# TODO: decide on a coding style
To improve readability, we maintain a consistent coding style in the library.
To improve readability, we maintain a consistent coding style in
the library.
- For C source files, we will use __(decide on a coding style)__
- For Fortran source files, we will use __(decide on a coding style)__
- For Fortran source files, we will use __(decide on a coding
style)__
Coding style can be automatically checked with [[https://clang.llvm.org/docs/ClangFormat.html][clang-format]].
@ -99,103 +95,105 @@ rm ${nb}.md
- deal with memory transfers between CPU and accelerators
- use different levels of floating-point precision
We chose a multi-layered design with low-level and high-level
We chose a multi-layered design with low-level and high-level
functions (see below).
*** Naming conventions
Use =qmckl_= as a prefix for all exported functions and variables.
All exported header files should have a filename with the prefix
All exported header files should have a filename with the prefix
=qmckl_=.
If the name of the org-mode file is =xxx.org=, the name of the
If the name of the org-mode file is =xxx.org=, the name of the
produced C files should be =xxx.c= and =xxx.h= and the name of the
produced Fortran files should be =xxx.f90=
Arrays are in uppercase and scalars are in lowercase.
In the names of the variables and functions, only the singular
form is allowed.
*** Application programming interface
The application programming interface (API) is designed to be
compatible with the C programming language (not C++), to ensure
that the library will be easily usable in any language.
This implies that only the following data types are allowed in the API:
The application programming interface (API) is designed to be
compatible with the C programming language (not C++), to ensure
that the library will be easily usable in /any/ language. This
implies that only the following data types are allowed in the API:
- 32-bit and 64-bit floats and arrays (=real= and =double=)
- 32-bit and 64-bit integers and arrays (=int32_t= and =int64_t=)
- Pointers should be represented as 64-bit integers (even on
- Pointers should be represented as 64-bit integers (even on
32-bit architectures)
- ASCII strings are represented as a pointers to a character arrays
and terminated by a zero character (C convention).
- ASCII strings are represented as a pointers to a character
arrays and terminated by a zero character (C convention).
Complex numbers can be represented by an array of 2 floats.
# TODO : Link to repositories for bindings
To facilitate the use in other languages than C, we provide some
To facilitate the use in other languages than C, we provide some
bindings in other languages in other repositories.
*** Global state
Global variables should be avoided in the library, because it is
possible that one single program needs to use multiple instances of
the library. To solve this problem we propose to use a pointer to a
=context= variable, built by the library with the
Global variables should be avoided in the library, because it is
possible that one single program needs to use multiple instances
of the library. To solve this problem we propose to use a pointer
to a =context= variable, built by the library with the
=qmckl_context_create= function. The =context= contains the global
state of the library, and is used as the first argument of many
state of the library, and is used as the first argument of many
QMCkl functions.
Modifying the state is done by setters and getters, prefixed
by =qmckl_context_set_= an =qmckl_context_get_=.
When a context variable is modified by a setter, a copy of the old
data structure is made and updated, and the pointer to the new data
structure is returned, such that the old contexts can still be
accessed.
It is also possible to modify the state in an impure fashion, using
the =qmckl_context_update_= functions.
The context and its old versions can be destroyed with
=qmckl_context_destroy=.
The internal structure of the context is not specified, to give a
maximum of freedom to the different implementations. Modifying
the state is done by setters and getters, prefixed by
=qmckl_context_set_= an =qmckl_context_get_=. When a context
variable is modified by a setter, a copy of the old data structure
is made and updated, and the pointer to the new data structure is
returned, such that the old contexts can still be accessed. It is
also possible to modify the state in an impure fashion, using the
=qmckl_context_update_= functions. The context and its old
versions can be destroyed with =qmckl_context_destroy=.
*** Low-level functions
Low-level functions are very simple functions which are leaves of the
function call tree (they don't call any other QMCkl function).
Low-level functions are very simple functions which are leaves of
the function call tree (they don't call any other QMCkl function).
This functions are /pure/, and unaware of the QMCkl =context=. They are
not allowed to allocate/deallocate memory, and if they need
temporary memory it should be provided in input.
These functions are /pure/, and unaware of the QMCkl
=context=. They are not allowed to allocate/deallocate memory, and
if they need temporary memory it should be provided in input.
*** High-level functions
High-level functions are at the top of the function call tree.
They are able to choose which lower-level function to call
High-level functions are at the top of the function call tree.
They are able to choose which lower-level function to call
depending on the required precision, and do the corresponding type
conversions.
These functions are also responsible for allocating temporary
storage, to simplify the use of accelerators.
conversions. These functions are also responsible for allocating
temporary storage, to simplify the use of accelerators.
The high-level functions should be pure, unless the introduction of
non-purity is justified. All the side effects should be made in the
=context= variable.
The high-level functions should be pure, unless the introduction
of non-purity is justified. All the side effects should be made in
the =context= variable.
# TODO : We need an identifier for impure functions
*** Numerical precision
The number of bits of precision required for a function should be
given as an input of low-level computational functions. This input will
be used to define the values of the different thresholds that might
be used to avoid computing unnecessary noise.
High-level functions will use the precision specified in the
=context= variable.
The number of bits of precision required for a function should be
given as an input of low-level computational functions. This input
will be used to define the values of the different thresholds that
might be used to avoid computing unnecessary noise. High-level
functions will use the precision specified in the =context=
variable.
** Algorithms
Reducing the scaling of an algorithm usually implies also reducing
its arithmetic complexity (number of flops per byte). Therefore,
for small sizes \(\mathcal{O}(N^3)\) and \(\mathcal{O}(N^2)\) algorithms
are better adapted than linear scaling algorithms.
As QMCkl is a general purpose library, multiple algorithms should
be implemented adapted to different problem sizes.
Reducing the scaling of an algorithm usually implies also reducing
its arithmetic complexity (number of flops per byte). Therefore,
for small sizes \(\mathcal{O}(N^3)\) and \(\mathcal{O}(N^2)\)
algorithms are better adapted than linear scaling algorithms. As
QMCkl is a general purpose library, multiple algorithms should be
implemented adapted to different problem sizes.
** Rules for the API

106
src/qmckl_ao.org
View File

@ -1,7 +1,7 @@
** Atomic Orbitals
This files contains all the routines for the computation of the
This files contains all the routines for the computation of the
values, gradients and Laplacian of the atomic basis functions.
3 files are produced:
@ -25,9 +25,12 @@ MunitResult test_qmckl_ao() {
P_l(\mathbf{r},\mathbf{R}_i) = (x-X_i)^a (y-Y_i)^b (z-Z_i)^c
\]
\begin{eqnarray*}
\frac{\partial }{\partial x} P_l\left(\mathbf{r},\mathbf{R}_i \right) & = & a (x-X_i)^{a-1} (y-Y_i)^b (z-Z_i)^c \\
\frac{\partial }{\partial y} P_l\left(\mathbf{r},\mathbf{R}_i \right) & = & b (x-X_i)^a (y-Y_i)^{b-1} (z-Z_i)^c \\
\frac{\partial }{\partial z} P_l\left(\mathbf{r},\mathbf{R}_i \right) & = & c (x-X_i)^a (y-Y_i)^b (z-Z_i)^{c-1} \\
\frac{\partial }{\partial x} P_l\left(\mathbf{r},\mathbf{R}_i \right) &
= & a (x-X_i)^{a-1} (y-Y_i)^b (z-Z_i)^c \\
\frac{\partial }{\partial y} P_l\left(\mathbf{r},\mathbf{R}_i \right) &
= & b (x-X_i)^a (y-Y_i)^{b-1} (z-Z_i)^c \\
\frac{\partial }{\partial z} P_l\left(\mathbf{r},\mathbf{R}_i \right) &
= & c (x-X_i)^a (y-Y_i)^b (z-Z_i)^{c-1} \\
\end{eqnarray*}
\begin{eqnarray*}
\left( \frac{\partial }{\partial x^2} +
@ -39,9 +42,9 @@ MunitResult test_qmckl_ao() {
&& c(c-1) (x-X_i)^a (y-Y_i)^b (z-Z_i)^{c-1}
\end{eqnarray*}
**** =qmckl_ao_powers=
**** =qmckl_ao_power=
Computes all the powers of the =n= input data up to the given
Computes all the powers of the =n= input data up to the given
maximum value given in input for each of the $n$ points:
\[ P_{ij} = X_j^i \]
@ -66,7 +69,7 @@ MunitResult test_qmckl_ao() {
***** Header
#+BEGIN_SRC C :tangle qmckl.h
qmckl_exit_code qmckl_ao_powers(const qmckl_context context,
qmckl_exit_code qmckl_ao_power(const qmckl_context context,
const int64_t n,
const double *X, const int32_t *LMAX,
const double *P, const int64_t LDP);
@ -74,7 +77,7 @@ qmckl_exit_code qmckl_ao_powers(const qmckl_context context,
***** Source
#+BEGIN_SRC f90 :tangle qmckl_ao.f90
integer function qmckl_ao_powers_f(context, n, X, LMAX, P, ldp) result(info)
integer function qmckl_ao_power_f(context, n, X, LMAX, P, ldp) result(info)
implicit none
integer*8 , intent(in) :: context
integer*8 , intent(in) :: n
@ -104,12 +107,12 @@ integer function qmckl_ao_powers_f(context, n, X, LMAX, P, ldp) result(info)
end do
end do
end function qmckl_ao_powers_f
end function qmckl_ao_power_f
#+END_SRC
***** C interface :noexport:
#+BEGIN_SRC f90 :tangle qmckl_ao.f90
integer(c_int32_t) function qmckl_ao_powers(context, n, X, LMAX, P, ldp) &
integer(c_int32_t) function qmckl_ao_power(context, n, X, LMAX, P, ldp) &
bind(C) result(info)
use, intrinsic :: iso_c_binding
implicit none
@ -120,14 +123,14 @@ integer(c_int32_t) function qmckl_ao_powers(context, n, X, LMAX, P, ldp) &
real (c_double) , intent(out) :: P(ldp,n)
integer (c_int64_t) , intent(in) , value :: ldp
integer, external :: qmckl_ao_powers_f
info = qmckl_ao_powers_f(context, n, X, LMAX, P, ldp)
end function qmckl_ao_powers
integer, external :: qmckl_ao_power_f
info = qmckl_ao_power_f(context, n, X, LMAX, P, ldp)
end function qmckl_ao_power
#+END_SRC
#+BEGIN_SRC f90 :tangle qmckl_f.f90
interface
integer(c_int32_t) function qmckl_ao_powers(context, n, X, LMAX, P, ldp) bind(C)
integer(c_int32_t) function qmckl_ao_power(context, n, X, LMAX, P, ldp) bind(C)
use, intrinsic :: iso_c_binding
integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(in) , value :: n
@ -135,13 +138,13 @@ end function qmckl_ao_powers
real (c_double) , intent(in) :: X(n)
integer (c_int32_t) , intent(in) :: LMAX(n)
real (c_double) , intent(out) :: P(ldp,n)
end function qmckl_ao_powers
end function qmckl_ao_power
end interface
#+END_SRC
***** Test :noexport:
#+BEGIN_SRC f90 :tangle test_qmckl_ao_f.f90
integer(c_int32_t) function test_qmckl_ao_powers(context) bind(C)
integer(c_int32_t) function test_qmckl_ao_power(context) bind(C)
use qmckl
implicit none
@ -165,10 +168,10 @@ integer(c_int32_t) function test_qmckl_ao_powers(context) bind(C)
LMAX(j) = 1 + int(mod(j, 5),4)
end do
test_qmckl_ao_powers = qmckl_ao_powers(context, n, X, LMAX, P, LDP)
if (test_qmckl_ao_powers /= 0) return
test_qmckl_ao_power = qmckl_ao_power(context, n, X, LMAX, P, LDP)
if (test_qmckl_ao_power /= 0) return
test_qmckl_ao_powers = -1
test_qmckl_ao_power = -1
do j=1,n
do i=1,LMAX(j)
@ -180,14 +183,14 @@ integer(c_int32_t) function test_qmckl_ao_powers(context) bind(C)
end do
end do
test_qmckl_ao_powers = 0
test_qmckl_ao_power = 0
deallocate(X,P,LMAX)
end function test_qmckl_ao_powers
end function test_qmckl_ao_power
#+END_SRC
#+BEGIN_SRC C :tangle test_qmckl_ao.c
int test_qmckl_ao_powers(qmckl_context context);
munit_assert_int(0, ==, test_qmckl_ao_powers(context));
int test_qmckl_ao_power(qmckl_context context);
munit_assert_int(0, ==, test_qmckl_ao_power(context));
#+END_SRC
@ -248,7 +251,7 @@ integer function qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL,
real*8 :: Y(3)
integer :: lmax_array(3)
real*8 :: pows(-2:lmax,3)
integer, external :: qmckl_ao_powers_f
integer, external :: qmckl_ao_power_f
double precision :: xy, yz, xz
double precision :: da, db, dc, dd
@ -281,11 +284,11 @@ integer function qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL,
pows(-2:-1,1:3) = 0.d0
pows(0,1:3) = 1.d0
lmax_array(1:3) = lmax
info = qmckl_ao_powers_f(context, 1_8, Y(1), (/lmax/), pows(1,1), size(pows,1,kind=8))
info = qmckl_ao_power_f(context, 1_8, Y(1), (/lmax/), pows(1,1), size(pows,1,kind=8))
if (info /= 0) return
info = qmckl_ao_powers_f(context, 1_8, Y(2), (/lmax/), pows(1,2), size(pows,1,kind=8))
info = qmckl_ao_power_f(context, 1_8, Y(2), (/lmax/), pows(1,2), size(pows,1,kind=8))
if (info /= 0) return
info = qmckl_ao_powers_f(context, 1_8, Y(3), (/lmax/), pows(1,3), size(pows,1,kind=8))
info = qmckl_ao_power_f(context, 1_8, Y(3), (/lmax/), pows(1,3), size(pows,1,kind=8))
if (info /= 0) return
@ -482,7 +485,7 @@ munit_assert_int(0, ==, test_qmckl_ao_polynomial_vgl(context));
*** Gaussian basis functions
**** =qmckl_ao_gaussians_vgl=
**** =qmckl_ao_gaussian_vgl=
Computes the values, gradients and Laplacians at a given point of
=n= Gaussian functions centered at the same point:
@ -516,7 +519,7 @@ munit_assert_int(0, ==, test_qmckl_ao_polynomial_vgl(context));
***** Header
#+BEGIN_SRC C :tangle qmckl.h
qmckl_exit_code qmckl_ao_gaussians_vgl(const qmckl_context context,
qmckl_exit_code qmckl_ao_gaussian_vgl(const qmckl_context context,
const double *X, const double *R,
const int64_t *n, const int64_t *A,
const double *VGL, const int64_t ldv);
@ -524,7 +527,7 @@ qmckl_exit_code qmckl_ao_gaussians_vgl(const qmckl_context context,
***** Source
#+BEGIN_SRC f90 :tangle qmckl_ao.f90
integer function qmckl_ao_gaussians_vgl_f(context, X, R, n, A, VGL, ldv) result(info)
integer function qmckl_ao_gaussian_vgl_f(context, X, R, n, A, VGL, ldv) result(info)
implicit none
integer*8 , intent(in) :: context
real*8 , intent(in) :: X(3), R(3)
@ -582,12 +585,12 @@ integer function qmckl_ao_gaussians_vgl_f(context, X, R, n, A, VGL, ldv) result(
VGL(i,5) = (t * A(i) - 6.d0) * VGL(i,5)
end do
end function qmckl_ao_gaussians_vgl_f
end function qmckl_ao_gaussian_vgl_f
#+END_SRC
***** C interface :noexport:
#+BEGIN_SRC f90 :tangle qmckl_ao.f90
integer(c_int32_t) function qmckl_ao_gaussians_vgl(context, X, R, n, A, VGL, ldv) &
integer(c_int32_t) function qmckl_ao_gaussian_vgl(context, X, R, n, A, VGL, ldv) &
bind(C) result(info)
use, intrinsic :: iso_c_binding
implicit none
@ -598,14 +601,14 @@ integer(c_int32_t) function qmckl_ao_gaussians_vgl(context, X, R, n, A, VGL, ldv
real (c_double) , intent(out) :: VGL(ldv,5)
integer (c_int64_t) , intent(in) , value :: ldv
integer, external :: qmckl_ao_gaussians_vgl_f
info = qmckl_ao_gaussians_vgl_f(context, X, R, n, A, VGL, ldv)
end function qmckl_ao_gaussians_vgl
integer, external :: qmckl_ao_gaussian_vgl_f
info = qmckl_ao_gaussian_vgl_f(context, X, R, n, A, VGL, ldv)
end function qmckl_ao_gaussian_vgl
#+END_SRC
#+BEGIN_SRC f90 :tangle qmckl_f.f90
interface
integer(c_int32_t) function qmckl_ao_gaussians_vgl(context, X, R, n, A, VGL, ldv) &
integer(c_int32_t) function qmckl_ao_gaussian_vgl(context, X, R, n, A, VGL, ldv) &
bind(C)
use, intrinsic :: iso_c_binding
integer (c_int64_t) , intent(in) , value :: context
@ -613,12 +616,12 @@ end function qmckl_ao_gaussians_vgl
integer (c_int64_t) , intent(in) , value :: n
real (c_double) , intent(in) :: X(3), R(3), A(n)
real (c_double) , intent(out) :: VGL(ldv,5)
end function qmckl_ao_gaussians_vgl
end function qmckl_ao_gaussian_vgl
end interface
#+END_SRC
***** Test :noexport:
#+BEGIN_SRC f90 :tangle test_qmckl_ao_f.f90
integer(c_int32_t) function test_qmckl_ao_gaussians_vgl(context) bind(C)
integer(c_int32_t) function test_qmckl_ao_gaussian_vgl(context) bind(C)
use qmckl
implicit none
@ -645,49 +648,48 @@ integer(c_int32_t) function test_qmckl_ao_gaussians_vgl(context) bind(C)
end do
test_qmckl_ao_gaussians_vgl = &
qmckl_ao_gaussians_vgl(context, X, R, n, A, VGL, ldv)
if (test_qmckl_ao_gaussians_vgl /= 0) return
test_qmckl_ao_gaussian_vgl = &
qmckl_ao_gaussian_vgl(context, X, R, n, A, VGL, ldv)
if (test_qmckl_ao_gaussian_vgl /= 0) return
test_qmckl_ao_gaussians_vgl = -1
test_qmckl_ao_gaussian_vgl = -1
do i=1,n
test_qmckl_ao_gaussians_vgl = -11
test_qmckl_ao_gaussian_vgl = -11
if (dabs(1.d0 - VGL(i,1) / (&
dexp(-A(i) * r2) &
)) > epsilon ) return
test_qmckl_ao_gaussians_vgl = -12
test_qmckl_ao_gaussian_vgl = -12
if (dabs(1.d0 - VGL(i,2) / (&
-2.d0 * A(i) * Y(1) * dexp(-A(i) * r2) &
)) > epsilon ) return
test_qmckl_ao_gaussians_vgl = -13
test_qmckl_ao_gaussian_vgl = -13
if (dabs(1.d0 - VGL(i,3) / (&
-2.d0 * A(i) * Y(2) * dexp(-A(i) * r2) &
)) > epsilon ) return
test_qmckl_ao_gaussians_vgl = -14
test_qmckl_ao_gaussian_vgl = -14
if (dabs(1.d0 - VGL(i,4) / (&
-2.d0 * A(i) * Y(3) * dexp(-A(i) * r2) &
)) > epsilon ) return
test_qmckl_ao_gaussians_vgl = -15
test_qmckl_ao_gaussian_vgl = -15
if (dabs(1.d0 - VGL(i,5) / (&
A(i) * (4.d0*r2*A(i) - 6.d0) * dexp(-A(i) * r2) &
)) > epsilon ) return
end do
test_qmckl_ao_gaussians_vgl = 0
test_qmckl_ao_gaussian_vgl = 0
deallocate(VGL)
end function test_qmckl_ao_gaussians_vgl
end function test_qmckl_ao_gaussian_vgl
#+END_SRC
#+BEGIN_SRC C :tangle test_qmckl_ao.c
int test_qmckl_ao_gaussians_vgl(qmckl_context context);
munit_assert_int(0, ==, test_qmckl_ao_gaussians_vgl(context));
#+END_SRC
int test_qmckl_ao_gaussian_vgl(qmckl_context context);
munit_assert_int(0, ==, test_qmckl_ao_gaussian_vgl(context));
#+END_SRC

374
src/qmckl_context.org
View File

@ -1,7 +1,7 @@
** Context
This file is written in C because it is more natural to express the context in
C than in Fortran.
This file is written in C because it is more natural to express the
context in C than in Fortran.
2 files are produced:
- a source file : =qmckl_context.c=
@ -21,43 +21,79 @@ MunitResult test_qmckl_context() {
*** Context
The context variable is a handle for the state of the library, and
is stored in the following data structure, which can't be seen
outside of the library. To simplify compatibility with other
is stored in the following data structure, which can't be seen
outside of the library. To simplify compatibility with other
languages, the pointer to the internal data structure is converted
into a 64-bit signed integer, defined in the =qmckl_context= type.
A value of 0 for the context is equivalent to a =NULL= pointer.
#+BEGIN_SRC C :comments org :tangle qmckl.h
# The following code block should be kept to insert comments into
# the qmckl.h file
#+BEGIN_SRC C :comments org :tangle qmckl.h :export none
#+END_SRC
***** Source
**** Basis set data structure
The tag is used internally to check if the memory domain pointed by
a pointer is a valid context.
Data structure for the info related to the atomic orbitals
basis set.
#+BEGIN_SRC C :comments org :tangle qmckl_context.c
#+BEGIN_SRC C :comments org :tangle qmckl_context.c
typedef struct qmckl_ao_basis_struct {
int64_t shell_num;
int64_t prim_num;
int64_t * shell_center;
int32_t * shell_ang_mom;
double * shell_factor;
double * exponent ;
double * coefficient ;
int64_t * shell_prim_num;
char type;
} qmckl_ao_basis_struct;
#+END_SRC
**** Source
The tag is used internally to check if the memory domain pointed
by a pointer is a valid context.
#+BEGIN_SRC C :comments org :tangle qmckl_context.c
typedef struct qmckl_context_struct {
struct qmckl_context_struct * prev;
/* Molecular system */
// struct qmckl_nucleus_struct * nucleus;
// struct qmckl_electron_struct * electron;
struct qmckl_ao_basis_struct * ao_basis;
// struct qmckl_mo_struct * mo;
// struct qmckl_determinant_struct * det;
/* Numerical precision */
uint32_t tag;
int32_t precision;
int32_t range;
} qmckl_context_struct;
#define VALID_TAG 0xBEEFFACE
#define INVALID_TAG 0xDEADBEEF
#+END_SRC
#+END_SRC
***** Test :noexport:
#+BEGIN_SRC C :tangle test_qmckl_context.c
**** Test :noexport:
#+BEGIN_SRC C :tangle test_qmckl_context.c
qmckl_context context;
qmckl_context new_context;
#+END_SRC
#+END_SRC
**** =qmckl_context_check=
Checks if the domain pointed by the pointer is a valid context.
Returns the input =qmckl_context= if the context is valid, 0 otherwise.
Checks if the domain pointed by the pointer is a valid context.
Returns the input =qmckl_context= if the context is valid, 0
otherwise.
#+BEGIN_SRC C :comments org :tangle qmckl.h
qmckl_context qmckl_context_check(const qmckl_context context) ;
@ -98,6 +134,7 @@ qmckl_context qmckl_context_create() {
}
context->prev = NULL;
context->ao_basis = NULL;
context->precision = QMCKL_DEFAULT_PRECISION;
context->range = QMCKL_DEFAULT_RANGE;
context->tag = VALID_TAG;
@ -153,6 +190,7 @@ qmckl_context qmckl_context_copy(const qmckl_context context) {
}
new_context->prev = old_context;
new_context->ao_basis = old_context->ao_basis;
new_context->precision = old_context->precision;
new_context->range = old_context->range;
new_context->tag = VALID_TAG;
@ -231,9 +269,9 @@ munit_assert_int64(qmckl_context_previous((qmckl_context) 0), ==, (qmckl_context
- Fails if the 0-valued context is given in argument
- Fails if the the pointer is not a valid context
#+BEGIN_SRC C :comments org :tangle qmckl.h
#+BEGIN_SRC C :comments org :tangle qmckl.h
qmckl_exit_code qmckl_context_destroy(qmckl_context context);
#+END_SRC
#+END_SRC
***** Source
#+BEGIN_SRC C :tangle qmckl_context.c
@ -271,16 +309,296 @@ munit_assert_int64(qmckl_context_check(new_context), ==, (qmckl_context) 0);
munit_assert_int64(qmckl_context_destroy((qmckl_context) 0), ==, QMCKL_FAILURE);
#+END_SRC
*** Basis set
For H_2 with the following basis set,
#+BEGIN_EXAMPLE
HYDROGEN
S 5
1 3.387000E+01 6.068000E-03
2 5.095000E+00 4.530800E-02
3 1.159000E+00 2.028220E-01
4 3.258000E-01 5.039030E-01
5 1.027000E-01 3.834210E-01
S 1
1 3.258000E-01 1.000000E+00
S 1
1 1.027000E-01 1.000000E+00
P 1
1 1.407000E+00 1.000000E+00
P 1
1 3.880000E-01 1.000000E+00
D 1
1 1.057000E+00 1.0000000
#+END_EXAMPLE
we have:
#+BEGIN_EXAMPLE
type = 'G'
shell_num = 12
prim_num = 20
SHELL_CENTER = [1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2]
SHELL_ANG_MOM = ['S', 'S', 'S', 'P', 'P', 'D', 'S', 'S', 'S', 'P', 'P', 'D']
SHELL_PRIM_NUM = [5, 1, 1, 1, 1, 1, 5, 1, 1, 1, 1, 1]
prim_index = [1, 6, 7, 8, 9, 10, 11, 16, 17, 18, 19, 20]
EXPONENT = [ 33.87, 5.095, 1.159, 0.3258, 0.1027, 0.3258, 0.1027,
1.407, 0.388, 1.057, 33.87, 5.095, 1.159, 0.3258, 0.1027,
0.3258, 0.1027, 1.407, 0.388, 1.057]
COEFFICIENT = [ 0.006068, 0.045308, 0.202822, 0.503903, 0.383421,
1.0, 1.0, 1.0, 1.0, 1.0, 0.006068, 0.045308, 0.202822,
0.503903, 0.383421, 1.0, 1.0, 1.0, 1.0, 1.0]
#+END_EXAMPLE
**** =qmckl_context_update_ao_basis=
Updates the data describing the AO basis set into the context.
| =type= | Gaussian or Slater |
| =shell_num= | Number of shells |
| =prim_num= | Total number of primitives |
| =SHELL_CENTER(shell_num)= | Id of the nucleus on which the shell is centered |
| =SHELL_ANG_MOM(shell_num)= | Id of the nucleus on which the shell is centered |
| =SHELL_FACTOR(shell_num)= | Normalization factor for the shell |
| =SHELL_PRIM_NUM(shell_num)= | Number of primitives in the shell |
| =SHELL_PRIM_INDEX(shell_num)= | Address of the first primitive of the shelll in the =EXPONENT= array |
| =EXPONENT(prim_num)= | Array of exponents |
| =COEFFICIENT(prim_num)= | Array of coefficients |
#+BEGIN_SRC C :comments org :tangle qmckl.h
qmckl_exit_code
qmckl_context_update_ao_basis(qmckl_context context , const char type,
const int64_t shell_num , const int64_t prim_num,
const int64_t * SHELL_CENTER, const int32_t * SHELL_ANG_MOM,
const double * SHELL_FACTOR, const int64_t * SHELL_PRIM_NUM,
const int64_t * SHELL_PRIM_INDEX,
const double * EXPONENT , const double * COEFFICIENT);
#+END_SRC
***** Source
#+BEGIN_SRC C :tangle qmckl_context.c
qmckl_exit_code
qmckl_context_update_ao_basis(qmckl_context context , const char type,
const int64_t shell_num , const int64_t prim_num,
const int64_t * SHELL_CENTER, const int32_t * SHELL_ANG_MOM,
const double * SHELL_FACTOR, const int64_t * SHELL_PRIM_NUM,
const int64_t * SHELL_PRIM_INDEX,
const double * EXPONENT , const double * COEFFICIENT)
{
int64_t i;
/* Check input */
if (type != 'G' && type != 'S') return QMCKL_FAILURE;
if (shell_num <= 0) return QMCKL_FAILURE;
if (prim_num <= 0) return QMCKL_FAILURE;
if (prim_num < shell_num) return QMCKL_FAILURE;
for (i=0 ; i<shell_num ; i++) {
if (SHELL_CENTER[i] <= 0) return QMCKL_FAILURE;
if (SHELL_PRIM_NUM[i] <= 0) return QMCKL_FAILURE;
if (SHELL_ANG_MOM[i] < 0) return QMCKL_FAILURE;
if (SHELL_PRIM_INDEX[i] < 0) return QMCKL_FAILURE;
}
for (i=0 ; i<prim_num ; i++) {
if (EXPONENT[i] <= 0) return QMCKL_FAILURE;
}
qmckl_context_struct* ctx = (qmckl_context_struct*) context;
if (ctx == NULL) return QMCKL_FAILURE;
qmckl_ao_basis_struct* basis = (qmckl_ao_basis_struct*) malloc (sizeof(qmckl_ao_basis_struct));
if (basis == NULL) return QMCKL_FAILURE;
/* Memory allocations */
basis->shell_center = (int64_t*) malloc (shell_num * sizeof(int64_t));
if (basis->shell_center == NULL) {
free(basis);
return QMCKL_FAILURE;
}
basis->shell_ang_mom = (int32_t*) malloc (shell_num * sizeof(int32_t));
if (basis->shell_ang_mom == NULL) {
free(basis->shell_center);
free(basis);
return QMCKL_FAILURE;
}
basis->shell_prim_num= (int64_t*) malloc (shell_num * sizeof(int64_t));
if (basis->shell_prim_num == NULL) {
free(basis->shell_ang_mom);
free(basis->shell_center);
free(basis);
return QMCKL_FAILURE;
}
basis->shell_factor = (double *) malloc (shell_num * sizeof(double ));
if (basis->shell_factor == NULL) {
free(basis->shell_prim_num);
free(basis->shell_ang_mom);
free(basis->shell_center);
free(basis);
return QMCKL_FAILURE;
}
basis->exponent = (double *) malloc (prim_num * sizeof(double ));
if (basis->exponent == NULL) {
free(basis->shell_factor);
free(basis->shell_prim_num);
free(basis->shell_ang_mom);
free(basis->shell_center);
free(basis);
return QMCKL_FAILURE;
}
basis->coefficient = (double *) malloc (prim_num * sizeof(double ));
if (basis->coefficient == NULL) {
free(basis->exponent);
free(basis->shell_factor);
free(basis->shell_prim_num);
free(basis->shell_ang_mom);
free(basis->shell_center);
free(basis);
return QMCKL_FAILURE;
}
/* Assign data */
basis->type = type;
basis->shell_num = shell_num;
basis->prim_num = prim_num;
for (i=0 ; i<shell_num ; i++) {
basis->shell_center [i] = SHELL_CENTER [i];
basis->shell_ang_mom [i] = SHELL_ANG_MOM [i];
basis->shell_prim_num[i] = SHELL_PRIM_NUM[i];
basis->shell_factor [i] = SHELL_FACTOR [i];
}
for (i=0 ; i<prim_num ; i++) {
basis->exponent [i] = EXPONENT[i];
basis->coefficient[i] = COEFFICIENT[i];
}
ctx->ao_basis = basis;
return QMCKL_SUCCESS;
}
#+END_SRC
***** Fortran interface
#+BEGIN_SRC f90 :tangle qmckl_f.f90
interface
integer (c_int32_t) function qmckl_context_update_ao_basis(context, &
typ, shell_num, prim_num, SHELL_CENTER, SHELL_ANG_MOM, SHELL_FACTOR, &
SHELL_PRIM_NUM, SHELL_PRIM_INDEX, EXPONENT, COEFFICIENT) bind(C)
use, intrinsic :: iso_c_binding
integer (c_int64_t), intent(in), value :: context
character(c_char) , intent(in), value :: typ
integer (c_int64_t), intent(in), value :: shell_num
integer (c_int64_t), intent(in), value :: prim_num
integer (c_int64_t), intent(in) :: SHELL_CENTER(shell_num)
integer (c_int32_t), intent(in) :: SHELL_ANG_MOM(shell_num)
double precision , intent(in) :: SHELL_FACTOR(shell_num)
integer (c_int64_t), intent(in) :: SHELL_PRIM_NUM(shell_num)
integer (c_int64_t), intent(in) :: SHELL_PRIM_INDEX(shell_num)
double precision , intent(in) :: EXPONENT(prim_num)
double precision , intent(in) :: COEFFICIENT(prim_num)
end function qmckl_context_update_ao_basis
end interface
#+END_SRC
***** TODO Test
**** =qmckl_context_set_ao_basis=
Sets the data describing the AO basis set into the context.
| =type= | Gaussian or Slater |
| =shell_num= | Number of shells |
| =prim_num= | Total number of primitives |
| =SHELL_CENTER(shell_num)= | Id of the nucleus on which the shell is centered |
| =SHELL_ANG_MOM(shell_num)= | Id of the nucleus on which the shell is centered |
| =SHELL_FACTOR(shell_num)= | Normalization factor for the shell |
| =SHELL_PRIM_NUM(shell_num)= | Number of primitives in the shell |
| =SHELL_PRIM_INDEX(shell_num)= | Address of the first primitive of the shelll in the =EXPONENT= array |
| =EXPONENT(prim_num)= | Array of exponents |
| =COEFFICIENT(prim_num)= | Array of coefficients |
#+BEGIN_SRC C :comments org :tangle qmckl.h
qmckl_context
qmckl_context_set_ao_basis(const qmckl_context context , const char type,
const int64_t shell_num , const int64_t prim_num,
const int64_t * SHELL_CENTER, const int32_t * SHELL_ANG_MOM,
const double * SHELL_FACTOR, const int64_t * SHELL_PRIM_NUM,
const int64_t * SHELL_PRIM_INDEX,
const double * EXPONENT , const double * COEFFICIENT);
#+END_SRC
***** Source
#+BEGIN_SRC C :tangle qmckl_context.c
qmckl_context
qmckl_context_set_ao_basis(const qmckl_context context , const char type,
const int64_t shell_num , const int64_t prim_num,
const int64_t * SHELL_CENTER, const int32_t * SHELL_ANG_MOM,
const double * SHELL_FACTOR, const int64_t * SHELL_PRIM_NUM,
const int64_t * SHELL_PRIM_INDEX,
const double * EXPONENT , const double * COEFFICIENT)
{
qmckl_context new_context = qmckl_context_copy(context);
if (new_context == 0) return 0;
if (qmckl_context_update_ao_basis(context, type, shell_num, prim_num,
SHELL_CENTER, SHELL_ANG_MOM, SHELL_FACTOR,
SHELL_PRIM_NUM, SHELL_PRIM_INDEX, EXPONENT,
COEFFICIENT
) == QMCKL_FAILURE)
return 0;
return new_context;
}
#+END_SRC
***** Fortran interface
#+BEGIN_SRC f90 :tangle qmckl_f.f90
interface
integer (c_int64_t) function qmckl_context_set_ao_basis(context, &
typ, shell_num, prim_num, SHELL_CENTER, SHELL_ANG_MOM, SHELL_FACTOR, &
SHELL_PRIM_NUM, SHELL_PRIM_INDEX, EXPONENT, COEFFICIENT) bind(C)
use, intrinsic :: iso_c_binding
integer (c_int64_t), intent(in), value :: context
character(c_char) , intent(in), value :: typ
integer (c_int64_t), intent(in), value :: shell_num
integer (c_int64_t), intent(in), value :: prim_num
integer (c_int64_t), intent(in) :: SHELL_CENTER(shell_num)
integer (c_int32_t), intent(in) :: SHELL_ANG_MOM(shell_num)
double precision , intent(in) :: SHELL_FACTOR(shell_num)
integer (c_int64_t), intent(in) :: SHELL_PRIM_NUM(shell_num)
integer (c_int64_t), intent(in) :: SHELL_PRIM_INDEX(shell_num)
double precision , intent(in) :: EXPONENT(prim_num)
double precision , intent(in) :: COEFFICIENT(prim_num)
end function qmckl_context_set_ao_basis
end interface
#+END_SRC
***** TODO Test
*** Precision
The following functions set and get the expected required precision
and range. =precision= should be an integer between 2 and 53, and
=range= should be an integer between 2 and 11.
The following functions set and get the expected required
precision and range. =precision= should be an integer between 2
and 53, and =range= should be an integer between 2 and 11.
The setter functions functions return a new context as a 64-bit integer.
The getter functions return the value, as a 32-bit integer.
The update functions return =QMCKL_SUCCESS= or =QMCKL_FAILURE=.
The setter functions functions return a new context as a 64-bit
integer. The getter functions return the value, as a 32-bit
integer. The update functions return =QMCKL_SUCCESS= or
=QMCKL_FAILURE=.
**** =qmckl_context_update_precision=
Modifies the parameter for the numerical precision in a given context.
@ -350,7 +668,7 @@ qmckl_exit_code qmckl_context_update_range(const qmckl_context context, const in
***** TODO Tests :noexport:
**** =qmckl_context_set_precision=
Returns a copy of the context with a different precision parameter.
#+BEGIN_SRC C :comments or :tangle qmckl.h
#+BEGIN_SRC C :comments org :tangle qmckl.h
qmckl_context qmckl_context_set_precision(const qmckl_context context, const int precision);
#+END_SRC
@ -369,7 +687,7 @@ qmckl_context qmckl_context_set_precision(const qmckl_context context, const int
***** Fortran interface
#+BEGIN_SRC f90 :tangle qmckl_f.f90
interface
integer (c_int32_t) function qmckl_context_set_precision(context, precision) bind(C)
integer (c_int64_t) function qmckl_context_set_precision(context, precision) bind(C)
use, intrinsic :: iso_c_binding
integer (c_int64_t), intent(in), value :: context
integer (c_int32_t), intent(in), value :: precision
@ -399,7 +717,7 @@ qmckl_context qmckl_context_set_range(const qmckl_context context, const int ran
***** Fortran interface
#+BEGIN_SRC f90 :tangle qmckl_f.f90
interface
integer (c_int32_t) function qmckl_context_set_range(context, range) bind(C)
integer (c_int64_t) function qmckl_context_set_range(context, range) bind(C)
use, intrinsic :: iso_c_binding
integer (c_int64_t), intent(in), value :: context
integer (c_int32_t), intent(in), value :: range
@ -486,11 +804,7 @@ double qmckl_context_get_epsilon(const qmckl_context context) {
***** TODO Tests :noexport:
*** Info about the molecular system
**** TODO =qmckl_context_set_nucl_coord=
**** TODO =qmckl_context_set_nucl_charge=
**** TODO =qmckl_context_set_elec_num=
*** End of files :noexport: