diff --git a/src/README.org b/src/README.org index 7ea4694..abe0663 100644 --- a/src/README.org +++ b/src/README.org @@ -1,58 +1,53 @@ #+TITLE: QMCkl source code documentation #+EXPORT_FILE_NAME: index.html -#+HTML_HEAD: -#+HTML_HEAD: -#+HTML_HEAD: -#+HTML_HEAD: -#+HTML_HEAD: -#+HTML_HEAD: +#+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 C 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 @@ -206,4 +204,4 @@ rm ${nb}.md * Documentation - + diff --git a/src/qmckl_ao.org b/src/qmckl_ao.org index 506beb8..13d28ee 100644 --- a/src/qmckl_ao.org +++ b/src/qmckl_ao.org @@ -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 diff --git a/src/qmckl_context.org b/src/qmckl_context.org index e407466..75a4bc8 100644 --- a/src/qmckl_context.org +++ b/src/qmckl_context.org @@ -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= @@ -18,46 +18,82 @@ MunitResult test_qmckl_context() { #+END_SRC -*** 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 - - The tag is used internally to check if the memory domain pointed by - a pointer is a valid context. +**** Basis set data structure - #+BEGIN_SRC C :comments org :tangle qmckl_context.c + Data structure for the info related to the atomic orbitals + basis set. + + #+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 ; ishell_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 ; ishell_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 ; iexponent [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: