TREX/qmckl
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2025-01-08 20:33:40 +01:00
Code Issues Projects Releases Wiki Activity

Merge branch 'TREX-CoE:master' into master

Browse Source
This commit is contained in:
vijay 2022-01-20 14:28:22 +01:00 committed by GitHub
commit 7ce2a98d8b
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
13 changed files with 1563 additions and 268 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
Makefile.am
View File

@ -53,7 +53,7 @@ src_qmckl_f = src/qmckl_f.f90
src_qmckl_fo = src/qmckl_f.o
header_tests = tests/chbrclf.h tests/n2.h
fortrandir = $(datadir)/fortran
fortrandir = $(datadir)/qmckl/fortran
fortran_DATA = $(src_qmckl_f)
QMCKL_TEST_DIR = $(abs_srcdir)/share/qmckl/test_data/

23
configure.ac
View File

@ -176,6 +176,22 @@ AX_BLAS([], [AC_MSG_ERROR([BLAS was not found.])])
## LAPACK
AX_LAPACK([], [AC_MSG_ERROR([LAPACK was not found.])])
# Specific options required with some compilers
case $FC in
ifort*)
FCFLAGS="$FCFLAGS -nofor-main"
;;
gfortran*)
# Order is important here
FCFLAGS="-cpp $FCFLAGS"
;;
esac
# Options.
AC_ARG_ENABLE(debug, [AS_HELP_STRING([--enable-debug],[compile for debugging])], ok=$enableval, ok=no)
@ -282,13 +298,6 @@ AC_ARG_ENABLE([vfc_ci],
esac],[vfc_ci=false])
AM_CONDITIONAL([VFC_CI], [test x$vfc_ci = xtrue])
# Enable Fortran preprocessor
if test "$FC" = "gfortran"; then
AC_MSG_NOTICE(gfortran detected)
# Arguments order is important here
FCFLAGS="-cpp $FCFLAGS"
fi
if test "$FC" = "verificarlo-f"; then
AC_MSG_NOTICE(verificarlo-f detected)
# Arguments order is important here

132
org/qmckl.org
View File

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

7
org/qmckl_ao.org
View File

@ -3312,7 +3312,7 @@ print ( "[7][4][26] : %e"% lf(a,x,y))
assert(qmckl_electron_provided(context));
rc = qmckl_set_electron_coord (context, 'N', elec_coord);
rc = qmckl_set_electron_coord (context, 'N', elec_coord, walk_num*3*elec_num);
assert(rc == QMCKL_SUCCESS);
@ -3726,7 +3726,7 @@ print ( "[1][4][26] : %25.15e"% lf(a,x,y))
assert(qmckl_electron_provided(context));
rc = qmckl_set_electron_coord (context, 'N', elec_coord);
rc = qmckl_set_electron_coord (context, 'N', elec_coord, walk_num*3*elec_num);
assert(rc == QMCKL_SUCCESS);
@ -4752,7 +4752,7 @@ assert (rc == QMCKL_SUCCESS);
assert(qmckl_electron_provided(context));
rc = qmckl_set_electron_coord (context, 'N', elec_coord);
rc = qmckl_set_electron_coord (context, 'N', elec_coord, walk_num*3*elec_num);
assert(rc == QMCKL_SUCCESS);
@ -4843,6 +4843,5 @@ assert( fabs(ao_vgl[1][26][224] - (-3.843864637762753e-09)) < 1.e-14 );
# vim: syntax=c
* TODO [0/1] Missing features :noexport:
- [ ] Error messages to tell what is missing when initializing

537
org/qmckl_blas.org
View File

@ -7,18 +7,543 @@
(org-babel-lob-ingest "../tools/lib.org")
#+end_src
#+begin_src c :tangle (eval h_private_type)
#ifndef QMCKL_BLAS_HPT
#define QMCKL_BLAS_HPT
#+end_src
#+begin_src c :tangle (eval h_private_func)
#ifndef QMCKL_BLAS_HPF
#define QMCKL_BLAS_HPF
#include "qmckl_blas_private_type.h"
#+end_src
#+begin_src c :tangle (eval c)
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#ifdef HAVE_STDINT_H
#include <stdint.h>
#elif HAVE_INTTYPES_H
#include <inttypes.h>
#endif
#include <stdlib.h>
#include <string.h>
#include <stdbool.h>
#include <assert.h>
#include <math.h>
#include "qmckl.h"
#include "qmckl_context_private_type.h"
#include "qmckl_memory_private_type.h"
#include "qmckl_memory_private_func.h"
#include "qmckl_blas_private_type.h"
#include "qmckl_blas_private_func.h"
#+end_src
#+begin_src c :comments link :tangle (eval c_test) :noweb yes
#include "qmckl.h"
#include "assert.h"
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include "qmckl_blas_private_type.h"
#include "qmckl_blas_private_func.h"
int main() {
qmckl_context context;
context = qmckl_context_create();
#+end_src
* Data types
** Vector
| Variable | Type | Description |
|----------+-----------+-------------------------|
| ~size~ | ~int64_t~ | Dimension of the vector |
| ~data~ | ~double*~ | Elements |
#+begin_src c :comments org :tangle (eval h_private_type)
typedef struct qmckl_vector {
int64_t size;
double* data;
} qmckl_vector;
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_vector
qmckl_vector_alloc( qmckl_context context,
const int64_t size);
#+end_src
Allocates a new vector. If the allocation failed the size is zero.
#+begin_src c :comments org :tangle (eval c)
qmckl_vector
qmckl_vector_alloc( qmckl_context context,
const int64_t size)
{
/* Should always be true by contruction */
assert (size > (int64_t) 0);
qmckl_vector result;
result.size = size;
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = size * sizeof(double);
result.data = (double*) qmckl_malloc (context, mem_info);
if (result.data == NULL) {
result.size = (int64_t) 0;
}
return result;
}
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_exit_code
qmckl_vector_free( qmckl_context context,
qmckl_vector vector);
#+end_src
#+begin_src c :comments org :tangle (eval c)
qmckl_exit_code
qmckl_vector_free( qmckl_context context,
qmckl_vector vector)
{
/* Always true */
assert (vector.data != NULL);
qmckl_exit_code rc;
rc = qmckl_free(context, vector.data);
if (rc != QMCKL_SUCCESS) {
return rc;
}
vector.size = (int64_t) 0;
vector.data = NULL;
return QMCKL_SUCCESS;
}
#+end_src
** Matrix
| Variable | Type | Description |
|----------+--------------+-----------------------------|
| ~size~ | ~int64_t[2]~ | Dimension of each component |
| ~data~ | ~double*~ | Elements |
The dimensions use Fortran ordering: two elements differing by one
in the first dimension are consecutive in memory.
#+begin_src c :comments org :tangle (eval h_private_type)
typedef struct qmckl_matrix {
int64_t size[2];
double* data;
} qmckl_matrix;
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_matrix
qmckl_matrix_alloc( qmckl_context context,
const int64_t size1,
const int64_t size2);
#+end_src
Allocates a new matrix. If the allocation failed the sizes are zero.
#+begin_src c :comments org :tangle (eval c)
qmckl_matrix
qmckl_matrix_alloc( qmckl_context context,
const int64_t size1,
const int64_t size2)
{
/* Should always be true by contruction */
assert (size1 * size2 > (int64_t) 0);
qmckl_matrix result;
result.size[0] = size1;
result.size[1] = size2;
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = size1 * size2 * sizeof(double);
result.data = (double*) qmckl_malloc (context, mem_info);
if (result.data == NULL) {
result.size[0] = (int64_t) 0;
result.size[1] = (int64_t) 0;
}
return result;
}
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_exit_code
qmckl_matrix_free( qmckl_context context,
qmckl_matrix matrix);
#+end_src
#+begin_src c :comments org :tangle (eval c)
qmckl_exit_code
qmckl_matrix_free( qmckl_context context,
qmckl_matrix matrix)
{
/* Always true */
assert (matrix.data != NULL);
qmckl_exit_code rc;
rc = qmckl_free(context, matrix.data);
if (rc != QMCKL_SUCCESS) {
return rc;
}
matrix.data = NULL;
matrix.size[0] = (int64_t) 0;
matrix.size[1] = (int64_t) 0;
return QMCKL_SUCCESS;
}
#+end_src
** Tensor
| Variable | Type | Description |
|----------+-----------------------------------+-----------------------------|
| ~order~ | ~int64_t~ | Order of the tensor |
| ~size~ | ~int64_t[QMCKL_TENSOR_ORDER_MAX]~ | Dimension of each component |
| ~data~ | ~double*~ | Elements |
The dimensions use Fortran ordering: two elements differing by one
in the first dimension are consecutive in memory.
#+begin_src c :comments org :tangle (eval h_private_type)
#define QMCKL_TENSOR_ORDER_MAX 16
typedef struct qmckl_tensor {
int64_t order;
int64_t size[QMCKL_TENSOR_ORDER_MAX];
double* data;
} qmckl_tensor;
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_tensor
qmckl_tensor_alloc( qmckl_context context,
const int64_t order,
const int64_t* size);
#+end_src
Allocates memory for a tensor. If the allocation failed, the size
is zero.
#+begin_src c :comments org :tangle (eval c)
qmckl_tensor
qmckl_tensor_alloc( qmckl_context context,
const int64_t order,
const int64_t* size)
{
/* Should always be true by contruction */
assert (order > 0);
assert (order <= QMCKL_TENSOR_ORDER_MAX);
assert (size != NULL);
qmckl_tensor result;
result.order = order;
int64_t prod_size = (int64_t) 1;
for (int64_t i=0 ; i<order ; ++i) {
assert (size[i] > (int64_t) 0);
result.size[i] = size[i];
prod_size *= size[i];
}
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = prod_size * sizeof(double);
result.data = (double*) qmckl_malloc (context, mem_info);
if (result.data == NULL) {
memset(&result, 0, sizeof(qmckl_tensor));
}
return result;
}
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_exit_code
qmckl_tensor_free( qmckl_context context,
qmckl_tensor tensor);
#+end_src
#+begin_src c :comments org :tangle (eval c)
qmckl_exit_code
qmckl_tensor_free( qmckl_context context,
qmckl_tensor tensor)
{
/* Always true */
assert (tensor.data != NULL);
qmckl_exit_code rc;
rc = qmckl_free(context, tensor.data);
if (rc != QMCKL_SUCCESS) {
return rc;
}
memset(&tensor, 0, sizeof(qmckl_tensor));
return QMCKL_SUCCESS;
}
#+end_src
** Reshaping
Reshaping occurs in-place and the pointer to the data is copied.
*** Vector -> Matrix
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_matrix
qmckl_matrix_of_vector(const qmckl_vector vector,
const int64_t size1,
const int64_t size2);
#+end_src
Reshapes a vector into a matrix.
#+begin_src c :comments org :tangle (eval c)
qmckl_matrix
qmckl_matrix_of_vector(const qmckl_vector vector,
const int64_t size1,
const int64_t size2)
{
/* Always true */
assert (size1 * size2 == vector.size);
qmckl_matrix result;
result.size[0] = size1;
result.size[1] = size2;
result.data = vector.data;
return result;
}
#+end_src
*** Vector -> Tensor
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_tensor
qmckl_tensor_of_vector(const qmckl_vector vector,
const int64_t order,
const int64_t* size);
#+end_src
Reshapes a vector into a tensor.
#+begin_src c :comments org :tangle (eval c)
qmckl_tensor
qmckl_tensor_of_vector(const qmckl_vector vector,
const int64_t order,
const int64_t* size)
{
qmckl_tensor result;
int64_t prod_size = 1;
for (int64_t i=0 ; i<order ; ++i) {
result.size[i] = size[i];
prod_size *= size[i];
}
assert (prod_size == vector.size);
result.data = vector.data;
return result;
}
#+end_src
*** Matrix -> Vector
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_vector
qmckl_vector_of_matrix(const qmckl_matrix matrix,
const int64_t size);
#+end_src
Reshapes a matrix into a vector.
#+begin_src c :comments org :tangle (eval c)
qmckl_vector
qmckl_vector_of_matrix(const qmckl_matrix matrix,
const int64_t size)
{
/* Always true */
assert (matrix.size[0] * matrix.size[1] == size);
qmckl_vector result;
result.size = size;
result.data = matrix.data;
return result;
}
#+end_src
*** Matrix -> Tensor
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_tensor
qmckl_tensor_of_matrix(const qmckl_matrix matrix,
const int64_t order,
const int64_t* size);
#+end_src
Reshapes a matrix into a tensor.
#+begin_src c :comments org :tangle (eval c)
qmckl_tensor
qmckl_tensor_of_matrix(const qmckl_matrix matrix,
const int64_t order,
const int64_t* size)
{
qmckl_tensor result;
int64_t prod_size = 1;
for (int64_t i=0 ; i<order ; ++i) {
result.size[i] = size[i];
prod_size *= size[i];
}
assert (prod_size == matrix.size[0] * matrix.size[1]);
result.data = matrix.data;
return result;
}
#+end_src
*** Tensor -> Vector
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_vector
qmckl_vector_of_tensor(const qmckl_tensor tensor,
const int64_t size);
#+end_src
Reshapes a tensor into a vector.
#+begin_src c :comments org :tangle (eval c)
qmckl_vector
qmckl_vector_of_tensor(const qmckl_tensor tensor,
const int64_t size)
{
/* Always true */
int64_t prod_size = (int64_t) 1;
for (int64_t i=0 ; i<tensor.order ; i++) {
prod_size *= tensor.size[i];
}
assert (prod_size == size);
qmckl_vector result;
result.size = size;
result.data = tensor.data;
return result;
}
#+end_src
*** Tensor -> Matrix
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_matrix
qmckl_matrix_of_tensor(const qmckl_tensor tensor,
const int64_t size1,
const int64_t size2);
#+end_src
Reshapes a tensor into a vector.
#+begin_src c :comments org :tangle (eval c)
qmckl_matrix
qmckl_matrix_of_tensor(const qmckl_tensor tensor,
const int64_t size1,
const int64_t size2)
{
/* Always true */
int64_t prod_size = (int64_t) 1;
for (int64_t i=0 ; i<tensor.order ; i++) {
prod_size *= tensor.size[i];
}
assert (prod_size == size1 * size2);
qmckl_matrix result;
result.size[0] = size1;
result.size[1] = size2;
result.data = tensor.data;
return result;
}
#+end_src
** Access macros
#+begin_src c :comments org :tangle (eval h_private_func)
#define qmckl_vec(v, i) v.data[i]
#define qmckl_mat(m, i, j) m.data[(i) + (j)*m.size[0]]
#define qmckl_ten3(t, i, j, k) t.data[(i) + m.size[0]*((j) + size[1]*(k))]
#define qmckl_ten4(t, i, j, k, l) t.data[(i) + m.size[0]*((j) + size[1]*((k) + size[2]*(l)))]
#define qmckl_ten5(t, i, j, k, l, m) t.data[(i) + m.size[0]*((j) + size[1]*((k) + size[2]*((l) + size[3]*(m))))]
#+end_src
** Tests
#+begin_src c :comments link :tangle (eval c_test)
{
int64_t m = 3;
int64_t n = 4;
int64_t p = m*n;
qmckl_vector vec = qmckl_vector_alloc(context, p);
for (int64_t i=0 ; i<p ; ++i)
qmckl_vec(vec, i) = (double) i;
for (int64_t i=0 ; i<p ; ++i)
assert( vec.data[i] == (double) i );
qmckl_matrix mat = qmckl_matrix_of_vector(vec, m, n);
assert (mat.size[0] == m);
assert (mat.size[1] == n);
assert (mat.data == vec.data);
for (int64_t j=0 ; j<n ; ++j)
for (int64_t i=0 ; i<m ; ++i)
assert ( qmckl_mat(mat, i, j) == qmckl_vec(vec, i+j*m)) ;
qmckl_vector vec2 = qmckl_vector_of_matrix(mat, p);
assert (vec2.size == p);
assert (vec2.data == vec.data);
for (int64_t i=0 ; i<p ; ++i)
assert ( qmckl_vec(vec2, i) == qmckl_vec(vec, i) ) ;
qmckl_vector_free(context, vec);
}
#+end_src
* Matrix operations
** ~qmckl_dgemm~
@ -443,8 +968,7 @@ subroutine adjugate4(a,LDA,B,LDB,na,det_l)
double precision :: C(4,4)
call cofactor4(A,LDA,B,4_8,na,det_l)
call cofactor4(A,LDA,C,4_8,na,det_l)
B(1,1) = C(1,1)
B(1,2) = C(2,1)
B(1,3) = C(3,1)
@ -1171,6 +1695,15 @@ assert(QMCKL_SUCCESS == test_qmckl_adjugate(context));
* End of files :noexport:
#+begin_src c :tangle (eval h_private_type)
#endif
#+end_src
#+begin_src c :tangle (eval h_private_func)
#endif
#+end_src
#+begin_src c :comments link :tangle (eval c_test)
assert (qmckl_context_destroy(context) == QMCKL_SUCCESS);
return 0;

8
org/qmckl_context.org
View File

@ -28,6 +28,7 @@ int main() {
#include "qmckl_error_private_type.h"
#include "qmckl_memory_private_type.h"
#include "qmckl_numprec_private_type.h"
#include "qmckl_point_private_type.h"
#include "qmckl_nucleus_private_type.h"
#include "qmckl_electron_private_type.h"
#include "qmckl_ao_private_type.h"
@ -35,6 +36,7 @@ int main() {
#include "qmckl_jastrow_private_type.h"
#include "qmckl_determinant_private_type.h"
#include "qmckl_local_energy_private_type.h"
#include "qmckl_point_private_func.h"
#include "qmckl_nucleus_private_func.h"
#include "qmckl_electron_private_func.h"
#include "qmckl_ao_private_func.h"
@ -121,6 +123,9 @@ typedef struct qmckl_context_struct {
/* Current date */
uint64_t date;
/* Points */
qmckl_point_struct *point;
/* -- Molecular system -- */
qmckl_nucleus_struct nucleus;
qmckl_electron_struct electron;
@ -236,6 +241,9 @@ qmckl_context qmckl_context_create() {
ctx->numprec.precision = QMCKL_DEFAULT_PRECISION;
ctx->numprec.range = QMCKL_DEFAULT_RANGE;
rc = qmckl_init_point(context);
assert (rc == QMCKL_SUCCESS);
rc = qmckl_init_electron(context);
assert (rc == QMCKL_SUCCESS);

142
org/qmckl_determinant.org
View File

@ -106,7 +106,7 @@ int main() {
| ~mo_index_beta~ | ~mo_index[det_num_beta][walk_num][beta_num]~ | Index of MOs for each walker |
Computed data:
|-----------------------------+------------------------------------------------------+-------------------------------------------------------------------------------------------|
| ~up_num~ | ~int64_t~ | Number of number of α electrons |
| ~donwn_num~ | ~int64_t~ | Number of number of β electrons |
@ -128,17 +128,17 @@ int main() {
| ~det_inv_matrix_beta~ | ~[det_num_beta][walk_num][beta_num][beta_num]~ | Inverse of the β electron slater matrix for each determinant of each walker. |
| ~det_inv_matrix_beta_date~ | ~int64_t~ | Date for the Inverse of the β electron slater matrix for each determinant of each walker. |
|-----------------------------+------------------------------------------------------+-------------------------------------------------------------------------------------------|
** Data structure
#+begin_src c :comments org :tangle (eval h_private_type)
typedef struct qmckl_determinant_struct {
char type;
int64_t walk_num;
int64_t det_num_alpha;
int64_t det_num_beta ;
int64_t up_num;
int64_t down_num;
int64_t walk_num;
int64_t det_num_alpha;
int64_t det_num_beta ;
int64_t up_num;
int64_t down_num;
int64_t* mo_index_alpha;
int64_t* mo_index_beta;
@ -171,10 +171,10 @@ typedef struct qmckl_determinant_struct {
Some values are initialized by default, and are not concerned by
this mechanism.
#+begin_src c :comments org :tangle (eval h_private_func)
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_exit_code qmckl_init_determinant(qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c)
qmckl_exit_code qmckl_init_determinant(qmckl_context context) {
@ -190,9 +190,9 @@ qmckl_exit_code qmckl_init_determinant(qmckl_context context) {
return QMCKL_SUCCESS;
}
#+end_src
** Access functions
#+begin_src c :comments org :tangle (eval h_private_func) :exports none
char qmckl_get_determinant_type (const qmckl_context context);
int64_t qmckl_get_determinant_walk_num (const qmckl_context context);
@ -201,7 +201,7 @@ int64_t qmckl_get_determinant_det_num_beta (const qmckl_context context);
int64_t* qmckl_get_determinant_mo_index_alpha (const qmckl_context context);
int64_t* qmckl_get_determinant_mo_index_beta (const qmckl_context context);
#+end_src
When all the data for the slater determinants have been provided, the following
function returns ~true~.
@ -307,7 +307,7 @@ int64_t qmckl_get_determinant_det_num_beta (const qmckl_context context) {
int64_t* qmckl_get_determinant_mo_index_alpha (const qmckl_context context) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return (int64_t) 0;
return NULL;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
@ -325,7 +325,7 @@ int64_t* qmckl_get_determinant_mo_index_alpha (const qmckl_context context) {
int64_t* qmckl_get_determinant_mo_index_beta (const qmckl_context context) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return (int64_t) 0;
return NULL;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
@ -334,7 +334,7 @@ int64_t* qmckl_get_determinant_mo_index_beta (const qmckl_context context) {
int32_t mask = 1 << 5;
if ( (ctx->det.uninitialized & mask) != 0) {
return (int64_t) 0;
return NULL;
}
assert (ctx->det.mo_index_beta != NULL);
@ -342,7 +342,7 @@ int64_t* qmckl_get_determinant_mo_index_beta (const qmckl_context context) {
}
#+end_src
** Initialization functions
To set the basis set, all the following functions need to be
@ -458,7 +458,7 @@ qmckl_exit_code qmckl_set_determinant_mo_index_alpha(qmckl_context context, con
}
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->det.walk_num * ctx->det.det_num_alpha *
mem_info.size = ctx->det.walk_num * ctx->det.det_num_alpha *
ctx->electron.up_num * sizeof(int64_t);
int64_t* new_array = (int64_t*) qmckl_malloc(context, mem_info);
if (new_array == NULL) {
@ -490,7 +490,7 @@ qmckl_exit_code qmckl_set_determinant_mo_index_beta(qmckl_context context, cons
}
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->det.walk_num * ctx->det.det_num_beta *
mem_info.size = ctx->det.walk_num * ctx->det.det_num_beta *
ctx->electron.down_num * sizeof(int64_t);
int64_t* new_array = (int64_t*) qmckl_malloc(context, mem_info);
if (new_array == NULL) {
@ -572,7 +572,7 @@ qmckl_exit_code qmckl_finalize_determinant(qmckl_context context) {
:END:
*** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_det_vgl_alpha(qmckl_context context, double* const det_vgl_alpha);
qmckl_exit_code qmckl_get_det_vgl_beta(qmckl_context context, double* const det_vgl_beta);
@ -580,7 +580,7 @@ qmckl_exit_code qmckl_get_det_vgl_beta(qmckl_context context, double* const det_
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_det_vgl_alpha(qmckl_context context, double * const det_vgl_alpha) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
@ -599,7 +599,7 @@ qmckl_exit_code qmckl_get_det_vgl_alpha(qmckl_context context, double * const de
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
size_t sze = 5 * ctx->det.det_num_alpha * ctx->det.walk_num *
size_t sze = 5 * ctx->det.det_num_alpha * ctx->det.walk_num *
ctx->electron.up_num * ctx->electron.up_num * sizeof(double);
memcpy(det_vgl_alpha, ctx->det.det_vgl_alpha, sze);
@ -607,7 +607,7 @@ qmckl_exit_code qmckl_get_det_vgl_alpha(qmckl_context context, double * const de
}
qmckl_exit_code qmckl_get_det_vgl_beta(qmckl_context context, double * const det_vgl_beta) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
@ -626,7 +626,7 @@ qmckl_exit_code qmckl_get_det_vgl_beta(qmckl_context context, double * const det
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
size_t sze = 5 * ctx->det.det_num_beta * ctx->det.walk_num *
size_t sze = 5 * ctx->det.det_num_beta * ctx->det.walk_num *
ctx->electron.down_num * ctx->electron.down_num * sizeof(double);
memcpy(det_vgl_beta, ctx->det.det_vgl_beta, sze);
@ -644,7 +644,7 @@ qmckl_exit_code qmckl_provide_det_vgl_beta(qmckl_context context);
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_provide_det_vgl_alpha(qmckl_context context) {
qmckl_exit_code rc;
qmckl_exit_code rc;
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
@ -658,14 +658,14 @@ qmckl_exit_code qmckl_provide_det_vgl_alpha(qmckl_context context) {
"qmckl_electron",
NULL);
}
if(!(ctx->electron.provided)) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_electron",
NULL);
}
if (!ctx->ao_basis.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
@ -731,7 +731,7 @@ qmckl_exit_code qmckl_provide_det_vgl_alpha(qmckl_context context) {
QMCKL_FAILURE,
"compute_det_vgl_alpha",
"Not yet implemented");
}
}
if (rc != QMCKL_SUCCESS) {
return rc;
}
@ -757,14 +757,14 @@ qmckl_exit_code qmckl_provide_det_vgl_beta(qmckl_context context) {
"qmckl_electron",
NULL);
}
if(!(ctx->electron.provided)) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_electron",
NULL);
}
if (!ctx->ao_basis.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
@ -806,7 +806,7 @@ qmckl_exit_code qmckl_provide_det_vgl_beta(qmckl_context context) {
ctx->det.det_vgl_beta = det_vgl_beta;
}
qmckl_exit_code rc;
qmckl_exit_code rc;
if (ctx->det.type == 'G') {
rc = qmckl_compute_det_vgl_beta(context,
ctx->det.det_num_beta,
@ -823,7 +823,7 @@ qmckl_exit_code qmckl_provide_det_vgl_beta(qmckl_context context) {
QMCKL_FAILURE,
"compute_det_vgl_beta",
"Not yet implemented");
}
}
if (rc != QMCKL_SUCCESS) {
return rc;
}
@ -854,7 +854,7 @@ qmckl_exit_code qmckl_provide_det_vgl_beta(qmckl_context context) {
| ~mo_num~ | ~int64_t~ | in | Number of MOs |
| ~mo_vgl~ | ~double[5][elec_num][mo_num]~ | in | Value, gradients and Laplacian of the MOs |
| ~det_vgl_alpha~ | ~double[det_num_alpha][walk_num][5][alpha_num][alpha_num]~ | out | Value, gradients and Laplacian of the Det |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_det_vgl_alpha_f(context, &
det_num_alpha, walk_num, alpha_num, beta_num, elec_num, &
@ -919,7 +919,7 @@ end function qmckl_compute_det_vgl_alpha_f
#+end_src
#+CALL: generate_c_header(table=qmckl_compute_det_vgl_alpha_args,rettyp=get_value("CRetType"),fname="qmckl_compute_det_vgl_alpha"))
#+RESULTS:
#+begin_src c :tangle (eval h_func) :comments org
qmckl_exit_code qmckl_compute_det_vgl_alpha (
@ -932,7 +932,7 @@ end function qmckl_compute_det_vgl_alpha_f
const int64_t* mo_index_alpha,
const int64_t mo_num,
const double* mo_vgl,
double* const det_vgl_alpha );
double* const det_vgl_alpha );
#+end_src
#+CALL: generate_c_interface(table=qmckl_compute_det_vgl_alpha_args,rettyp=get_value("CRetType"),fname="qmckl_compute_det_vgl_alpha"))
@ -981,7 +981,7 @@ end function qmckl_compute_det_vgl_alpha_f
end function qmckl_compute_det_vgl_alpha
#+end_src
*** Compute beta
:PROPERTIES:
:Name: qmckl_compute_det_vgl_beta
@ -1002,7 +1002,7 @@ end function qmckl_compute_det_vgl_alpha_f
| ~mo_num~ | ~int64_t~ | in | Number of MOs |
| ~mo_vgl~ | ~double[5][elec_num][mo_num]~ | in | Value, gradients and Laplacian of the MOs |
| ~det_vgl_beta~ | ~double[det_num_beta][walk_num][5][beta_num][beta_num]~ | out | Value, gradients and Laplacian of the Det |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_det_vgl_beta_f(context, &
det_num_beta, walk_num, alpha_num, beta_num, elec_num, &
@ -1080,9 +1080,9 @@ end function qmckl_compute_det_vgl_beta_f
const int64_t* mo_index_beta,
const int64_t mo_num,
const double* mo_vgl,
double* const det_vgl_beta );
double* const det_vgl_beta );
#+end_src
#+CALL: generate_c_interface(table=qmckl_compute_det_vgl_beta_args,rettyp=get_value("CRetType"),fname="qmckl_compute_det_vgl_beta"))
#+RESULTS:
@ -1129,7 +1129,7 @@ end function qmckl_compute_det_vgl_beta_f
end function qmckl_compute_det_vgl_beta
#+end_src
*** Test
#+begin_src c :tangle (eval c_test) :exports none
@ -1146,7 +1146,7 @@ const int64_t nucl_num = chbrclf_nucl_num;
const double* nucl_charge = chbrclf_charge;
const double* nucl_coord = &(chbrclf_nucl_coord[0][0]);
rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num);
rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num);
assert (rc == QMCKL_SUCCESS);
rc = qmckl_set_electron_walk_num (context, walk_num);
@ -1154,7 +1154,7 @@ assert (rc == QMCKL_SUCCESS);
assert(qmckl_electron_provided(context));
rc = qmckl_set_electron_coord (context, 'N', elec_coord);
rc = qmckl_set_electron_coord (context, 'N', elec_coord, walk_num*elec_num*3);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_num (context, nucl_num);
@ -1307,7 +1307,7 @@ rc = qmckl_get_det_vgl_beta(context, &(det_vgl_beta[0][0][0][0][0]));
assert (rc == QMCKL_SUCCESS);
#+end_src
** Inverse of Determinant matrix
:PROPERTIES:
:Name: qmckl_compute_det_inv_matrix
@ -1316,7 +1316,7 @@ assert (rc == QMCKL_SUCCESS);
:END:
*** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_det_inv_matrix_alpha(qmckl_context context, double* const det_inv_matrix_alpha);
qmckl_exit_code qmckl_get_det_inv_matrix_beta(qmckl_context context, double* const det_inv_matrix_beta);
@ -1328,7 +1328,7 @@ qmckl_exit_code qmckl_get_det_beta(qmckl_context context, double* const det_adj_
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_det_inv_matrix_alpha(qmckl_context context, double * const det_inv_matrix_alpha) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
@ -1357,7 +1357,7 @@ qmckl_exit_code qmckl_get_det_inv_matrix_alpha(qmckl_context context, double * c
}
qmckl_exit_code qmckl_get_det_inv_matrix_beta(qmckl_context context, double * const det_inv_matrix_beta) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
@ -1386,7 +1386,7 @@ qmckl_exit_code qmckl_get_det_inv_matrix_beta(qmckl_context context, double * co
}
qmckl_exit_code qmckl_get_det_adj_matrix_alpha(qmckl_context context, double * const det_adj_matrix_alpha) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
@ -1415,7 +1415,7 @@ qmckl_exit_code qmckl_get_det_adj_matrix_alpha(qmckl_context context, double * c
}
qmckl_exit_code qmckl_get_det_adj_matrix_beta(qmckl_context context, double * const det_adj_matrix_beta) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
@ -1444,7 +1444,7 @@ qmckl_exit_code qmckl_get_det_adj_matrix_beta(qmckl_context context, double * co
}
qmckl_exit_code qmckl_get_det_alpha(qmckl_context context, double * const det_value_alpha) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
@ -1473,7 +1473,7 @@ qmckl_exit_code qmckl_get_det_alpha(qmckl_context context, double * const det_va
}
qmckl_exit_code qmckl_get_det_beta(qmckl_context context, double * const det_value_beta) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
@ -1502,7 +1502,7 @@ qmckl_exit_code qmckl_get_det_beta(qmckl_context context, double * const det_val
}
#+end_src
*** Provide
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
@ -1526,14 +1526,14 @@ qmckl_exit_code qmckl_provide_det_inv_matrix_alpha(qmckl_context context) {
"qmckl_electron",
NULL);
}
if(!(ctx->electron.provided)) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_electron",
NULL);
}
if (!ctx->ao_basis.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
@ -1562,7 +1562,7 @@ qmckl_exit_code qmckl_provide_det_inv_matrix_alpha(qmckl_context context) {
if (ctx->det.det_inv_matrix_alpha == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->det.walk_num * ctx->det.det_num_alpha *
mem_info.size = ctx->det.walk_num * ctx->det.det_num_alpha *
ctx->electron.up_num * ctx->electron.up_num * sizeof(double);
double* det_inv_matrix_alpha = (double*) qmckl_malloc(context, mem_info);
@ -1578,7 +1578,7 @@ qmckl_exit_code qmckl_provide_det_inv_matrix_alpha(qmckl_context context) {
if (ctx->det.det_adj_matrix_alpha == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->det.walk_num * ctx->det.det_num_alpha *
mem_info.size = ctx->det.walk_num * ctx->det.det_num_alpha *
ctx->electron.up_num * ctx->electron.up_num * sizeof(double);
double* det_adj_matrix_alpha = (double*) qmckl_malloc(context, mem_info);
@ -1606,7 +1606,7 @@ qmckl_exit_code qmckl_provide_det_inv_matrix_alpha(qmckl_context context) {
ctx->det.det_value_alpha = det_value_alpha;
}
qmckl_exit_code rc;
qmckl_exit_code rc;
if (ctx->det.type == 'G') {
rc = qmckl_compute_det_inv_matrix_alpha(context,
ctx->det.det_num_alpha,
@ -1621,7 +1621,7 @@ qmckl_exit_code qmckl_provide_det_inv_matrix_alpha(qmckl_context context) {
QMCKL_FAILURE,
"compute_det_inv_matrix_alpha",
"Not yet implemented");
}
}
if (rc != QMCKL_SUCCESS) {
return rc;
}
@ -1649,14 +1649,14 @@ qmckl_exit_code qmckl_provide_det_inv_matrix_beta(qmckl_context context) {
"qmckl_electron",
NULL);
}
if(!(ctx->electron.provided)) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_electron",
NULL);
}
if (!ctx->ao_basis.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
@ -1685,7 +1685,7 @@ qmckl_exit_code qmckl_provide_det_inv_matrix_beta(qmckl_context context) {
if (ctx->det.det_inv_matrix_beta == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->det.walk_num * ctx->det.det_num_beta *
mem_info.size = ctx->det.walk_num * ctx->det.det_num_beta *
ctx->electron.down_num * ctx->electron.down_num * sizeof(double);
double* det_inv_matrix_beta = (double*) qmckl_malloc(context, mem_info);
@ -1701,7 +1701,7 @@ qmckl_exit_code qmckl_provide_det_inv_matrix_beta(qmckl_context context) {
if (ctx->det.det_adj_matrix_beta == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->det.walk_num * ctx->det.det_num_beta *
mem_info.size = ctx->det.walk_num * ctx->det.det_num_beta *
ctx->electron.down_num * ctx->electron.down_num * sizeof(double);
double* det_adj_matrix_beta = (double*) qmckl_malloc(context, mem_info);
@ -1729,7 +1729,7 @@ qmckl_exit_code qmckl_provide_det_inv_matrix_beta(qmckl_context context) {
ctx->det.det_value_beta = det_value_beta;
}
qmckl_exit_code rc;
qmckl_exit_code rc;
if (ctx->det.type == 'G') {
rc = qmckl_compute_det_inv_matrix_beta(context,
ctx->det.det_num_beta,
@ -1744,7 +1744,7 @@ qmckl_exit_code qmckl_provide_det_inv_matrix_beta(qmckl_context context) {
QMCKL_FAILURE,
"compute_det_inv_matrix_beta",
"Not yet implemented");
}
}
if (rc != QMCKL_SUCCESS) {
return rc;
}
@ -1757,7 +1757,7 @@ qmckl_exit_code qmckl_provide_det_inv_matrix_beta(qmckl_context context) {
return QMCKL_SUCCESS;
}
#+end_src
*** Compute alpha
:PROPERTIES:
:Name: qmckl_compute_det_inv_matrix_alpha
@ -1776,7 +1776,7 @@ qmckl_exit_code qmckl_provide_det_inv_matrix_beta(qmckl_context context) {
| ~det_value_alpha~ | ~double[det_num_alpha][walk_num]~ | out | value of determinant matrix |
| ~det_adj_matrix_alpha~ | ~double[det_num_alpha][walk_num][alpha_num][alpha_num]~ | out | adjoint of determinant matrix |
| ~det_inv_matrix_alpha~ | ~double[det_num_alpha][walk_num][alpha_num][alpha_num]~ | out | inverse of determinant matrix |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_det_inv_matrix_alpha_f(context, &
det_num_alpha, walk_num, alpha_num, det_vgl_alpha, det_value_alpha, det_adj_matrix_alpha, det_inv_matrix_alpha) &
@ -1828,13 +1828,13 @@ integer function qmckl_compute_det_inv_matrix_alpha_f(context, &
res = qmckl_adjugate(context, &
alpha_num, matA, LDA, &
det_adj_matrix_alpha(1, 1, iwalk, idet), &
det_adj_matrix_alpha(1:alpha_num, 1:alpha_num, iwalk, idet), &
int(size(det_adj_matrix_alpha,1),8), &
det_l)
det_inv_matrix_alpha(1:alpha_num, 1:alpha_num, iwalk, idet) = &
(1.d0/det_l) * &
det_adj_matrix_alpha(1:alpha_num, 1:alpha_num, iwalk, idet)
det_adj_matrix_alpha(1:alpha_num, 1:alpha_num, iwalk, idet)
det_value_alpha(iwalk, idet) = det_l
end do
@ -1856,7 +1856,7 @@ end function qmckl_compute_det_inv_matrix_alpha_f
const double* det_vgl_alpha,
double* const det_value_alpha,
double* const det_adj_matrix_alpha,
double* const det_inv_matrix_alpha );
double* const det_inv_matrix_alpha );
#+end_src
#+CALL: generate_c_interface(table=qmckl_det_inv_matrix_alpha_args,rettyp=get_value("CRetType"),fname="qmckl_compute_det_inv_matrix_alpha"))
@ -1918,7 +1918,7 @@ end function qmckl_compute_det_inv_matrix_alpha_f
| ~det_value_beta~ | ~double[det_num_beta][walk_num]~ | out | value of determinant matrix |
| ~det_adj_matrix_beta~ | ~double[det_num_beta][walk_num][beta_num][beta_num]~ | out | adjoint of determinant matrix |
| ~det_inv_matrix_beta~ | ~double[det_num_beta][walk_num][beta_num][beta_num]~ | out | inverse of determinant matrix |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_det_inv_matrix_beta_f(context, &
det_num_beta, walk_num, beta_num, det_vgl_beta, det_value_beta, det_adj_matrix_beta, det_inv_matrix_beta) &
@ -1976,7 +1976,7 @@ integer function qmckl_compute_det_inv_matrix_beta_f(context, &
det_inv_matrix_beta(1:beta_num, 1:beta_num, iwalk, idet) = &
(1.d0/det_l) * &
det_adj_matrix_beta(1:beta_num, 1:beta_num, iwalk, idet)
det_adj_matrix_beta(1:beta_num, 1:beta_num, iwalk, idet)
det_value_beta(iwalk, idet) = det_l
end do
@ -1999,7 +1999,7 @@ end function qmckl_compute_det_inv_matrix_beta_f
const double* det_vgl_beta,
double* const det_value_beta,
double* const det_adj_matrix_beta,
double* const det_inv_matrix_beta );
double* const det_inv_matrix_beta );
#+end_src
#+CALL: generate_c_interface(table=qmckl_det_inv_matrix_beta_args,rettyp=get_value("CRetType"),fname="qmckl_compute_det_inv_matrix_beta"))

169
org/qmckl_electron.org
View File

@ -394,7 +394,7 @@ qmckl_get_electron_rescale_factor_en (const qmckl_context context, double* const
*** Electron coordinates
Returns the current electron coordinates. The pointer is assumed
to point on a memory block of size ~3 * elec_num * walk_num~.
to point on a memory block of size ~size_max~ \ge ~3 * elec_num * walk_num~.
The order of the indices is:
| | Normal | Transposed |
@ -479,7 +479,7 @@ qmckl_get_electron_coord (const qmckl_context context, const char transp, double
#+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code qmckl_set_electron_num (qmckl_context context, const int64_t up_num, const int64_t down_num);
qmckl_exit_code qmckl_set_electron_walk_num (qmckl_context context, const int64_t walk_num);
qmckl_exit_code qmckl_set_electron_coord (qmckl_context context, const char transp, const double* coord);
qmckl_exit_code qmckl_set_electron_coord (qmckl_context context, const char transp, const double* coord, const int64_t size_max);
qmckl_exit_code qmckl_set_electron_rescale_factor_ee (qmckl_context context, const double kappa_ee);
qmckl_exit_code qmckl_set_electron_rescale_factor_en (qmckl_context context, const double kappa_en);
@ -664,7 +664,11 @@ end interface
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_set_electron_coord(qmckl_context context, const char transp, const double* coord) {
qmckl_set_electron_coord(qmckl_context context,
const char transp,
const double* coord,
const int64_t size_max)
{
<<pre2>>
@ -705,6 +709,13 @@ qmckl_set_electron_coord(qmckl_context context, const char transp, const double*
"walk_num is not set");
}
if (size_max < walk_num*3*elec_num) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_4,
"qmckl_set_electron_coord",
"destination array is too small");
}
/* If num and walk_num are set, the arrays should be allocated */
assert (ctx->electron.coord_old != NULL);
assert (ctx->electron.coord_new != NULL);
@ -742,7 +753,7 @@ qmckl_set_electron_coord(qmckl_context context, const char transp, const double*
#+begin_src f90 :comments org :tangle (eval fh_func) :noweb yes
interface
integer(c_int32_t) function qmckl_set_electron_coord(context, transp, coord) bind(C)
integer(c_int32_t) function qmckl_set_electron_coord(context, transp, coord, size_max) bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
@ -750,6 +761,7 @@ interface
integer (c_int64_t) , intent(in) , value :: context
character , intent(in) , value :: transp
double precision , intent(in) :: coord(*)
integer (c_int64_t) , intent(in) , value :: size_max
end function
end interface
#+end_src
@ -848,7 +860,7 @@ assert(w == walk_num);
assert(qmckl_electron_provided(context));
rc = qmckl_set_electron_coord (context, 'N', elec_coord);
rc = qmckl_set_electron_coord (context, 'N', elec_coord, walk_num*elec_num*3);
assert(rc == QMCKL_SUCCESS);
double elec_coord2[walk_num*3*elec_num];
@ -1219,14 +1231,14 @@ qmckl_exit_code qmckl_provide_ee_distance_rescaled(qmckl_context context)
:END:
#+NAME: qmckl_ee_distance_rescaled_args
| Variable | Type | In/Out | Description |
|----------------------------------------+---------------------------+--------+--------------------------------------|
| ~context~ | ~qmckl_context~ | in | Global state |
| ~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 |
| ~ee_distance~ | ~double[walk_num][elec_num][elec_num]~ | out | Electron-electron rescaled distances |
| Variable | Type | In/Out | Description |
|---------------------------+----------------------------------------+--------+--------------------------------------|
| ~context~ | ~qmckl_context~ | in | Global state |
| ~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 |
| ~ee_distance~ | ~double[walk_num][elec_num][elec_num]~ | out | Electron-electron rescaled distances |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_ee_distance_rescaled_f(context, elec_num, rescale_factor_kappa_ee, walk_num, &
@ -2077,7 +2089,7 @@ assert(fabs(en_distance[1][0][1] - 3.1804527583077356) < 1.e-12);
** Electron-nucleus rescaled distances
~en_distance_rescaled~ stores the matrix of the rescaled distances between
electrons and nucleii.
electrons and nuclei.
\[
C_{ij} = \left( 1 - \exp{-\kappa C_{ij}}\right)/\kappa
@ -2746,7 +2758,7 @@ qmckl_exit_code qmckl_provide_en_potential(qmckl_context context)
|---------------+----------------------------------------+--------+--------------------------------------|
| ~context~ | ~qmckl_context~ | in | Global state |
| ~elec_num~ | ~int64_t~ | in | Number of electrons |
| ~nucl_num~ | ~int64_t~ | in | Number of nucleii |
| ~nucl_num~ | ~int64_t~ | in | Number of nuclei |
| ~walk_num~ | ~int64_t~ | in | Number of walkers |
| ~charge~ | ~double[nucl_num]~ | in | charge of nucleus |
| ~en_distance~ | ~double[walk_num][nucl_num][elec_num]~ | in | Electron-electron rescaled distances |
@ -2847,6 +2859,133 @@ rc = qmckl_get_electron_en_potential(context, &(en_pot[0]));
assert (rc == QMCKL_SUCCESS);
#+end_src
** Generate initial coordinates
*** Compute :noexport:
# begin_src f90 :comments org :tangle (eval f) :noweb yes
subroutine draw_init_points
implicit none
BEGIN_DOC
! Place randomly electrons around nuclei
END_DOC
integer :: iwalk
logical, allocatable :: do_elec(:)
integer :: acc_num
real, allocatable :: xmin(:,:)
integer :: i, j, k, l, kk
real :: norm
allocate (do_elec(elec_num), xmin(3,elec_num))
xmin = -huge(1.)
norm = 0.
do i=1,elec_alpha_num
do j=1,ao_num
norm += mo_coef_transp(i,j)*mo_coef_transp(i,j)
enddo
enddo
norm = sqrt(norm/float(elec_alpha_num))
call rinfo( irp_here, 'Norm : ', norm )
call rinfo( irp_here, 'mo_scale: ' , mo_scale )
mo_coef_transp = mo_coef_transp/norm
double precision :: qmc_ranf
real :: mo_max
do i=1,elec_alpha_num
l=1
xmin(1,i) = mo_coef_transp(i,1)*mo_coef_transp(i,1) - 0.001*qmc_ranf()
do j=2,ao_num
xmin(2,i) = mo_coef_transp(i,j)*mo_coef_transp(i,j) - 0.001*qmc_ranf()
if (xmin(2,i) > xmin(1,i) ) then
xmin(1,i) = xmin(2,i)
l = ao_nucl(j)
endif
enddo
xmin(1,i) = nucl_coord(l,1)
xmin(2,i) = nucl_coord(l,2)
xmin(3,i) = nucl_coord(l,3)
enddo
call iinfo(irp_here, 'Det num = ', det_num )
do k=1,elec_beta_num
i = k+elec_alpha_num
l=1
xmin(1,i) = mo_coef_transp(k,1)*mo_coef_transp(k,1) - 0.001*qmc_ranf()
do j=2,ao_num
xmin(2,i) = mo_coef_transp(k,j)*mo_coef_transp(k,j) - 0.001*qmc_ranf()
if (xmin(2,i) > xmin(1,i) ) then
xmin(1,i) = xmin(2,i)
l = ao_nucl(j)
endif
enddo
xmin(1,i) = nucl_coord(l,1)
xmin(2,i) = nucl_coord(l,2)
xmin(3,i) = nucl_coord(l,3)
enddo
call rinfo( irp_here, 'time step =', time_step )
do iwalk=1,walk_num
print *, 'Generating initial positions for walker', iwalk
acc_num = 0
do_elec = .True.
integer :: iter
do iter = 1,10000
if (acc_num >= elec_num) then
exit
endif
double precision :: gauss
real :: re_compute
re_compute = 0.
do while (re_compute < 1.e-6)
do i=1,elec_num
if (do_elec(i)) then
do l=1,3
elec_coord(i,l) = xmin(l,i) + 1.5*(0.5-qmc_ranf())
enddo
endif
enddo
TOUCH elec_coord
re_compute = minval(nucl_elec_dist(1:nucl_num,1:elec_num))
enddo
do i=1,elec_alpha_num
if (do_elec(i)) then
if ( mo_value_transp(i,i)**2 >= qmc_ranf()) then
acc_num += 1
do_elec(i) = .False.
endif
endif
enddo
do i=1,elec_beta_num
if (do_elec(i+elec_alpha_num)) then
if ( mo_value_transp(i,i+elec_alpha_num)**2 >= qmc_ranf()) then
acc_num += 1
do_elec(i+elec_alpha_num) = .False.
endif
endif
enddo
enddo
do l=1,3
do i=1,elec_num+1
elec_coord_full(i,l,iwalk) = elec_coord(i,l)
enddo
enddo
enddo
if (.not.is_worker) then
call ezfio_set_electrons_elec_coord_pool_size(walk_num)
call ezfio_set_electrons_elec_coord_pool(elec_coord_full)
endif
SOFT_TOUCH elec_coord elec_coord_full
deallocate (do_elec, xmin)
end
# end_src
* End of files :noexport:
#+begin_src c :tangle (eval h_private_type)

42
org/qmckl_jastrow.org
View File

@ -103,26 +103,26 @@ int main() {
computed data:
| Variable | Type | In/Out | Description |
|------------+-----------------------------------------------------------------------+---------------------------------------------------------------------------------------------------------+-------------|
| ~dim_cord_vect~ | ~int64_t~ | Number of unique C coefficients | |
| ~dim_cord_vect_date~ | ~uint64_t~ | Number of unique C coefficients | |
| ~asymp_jasb~ | ~double[2]~ | Asymptotic component | |
| ~asymp_jasb_date~ | ~uint64_t~ | Asymptotic component | |
| ~cord_vect_full~ | ~double[dim_cord_vect][nucl_num]~ | vector of non-zero coefficients | |
| ~cord_vect_full_date~ | ~uint64_t~ | Keep track of changes here | |
| ~lkpm_combined_index~ | ~int64_t[4][dim_cord_vect]~ | Transform l,k,p, and m into consecutive indices | |
| ~lkpm_combined_index_date~ | ~uint64_t~ | Transform l,k,p, and m into consecutive indices | |
| ~tmp_c~ | ~double[elec_num][nucl_num][ncord + 1][ncord][walk_num]~ | vector of non-zero coefficients | |
| ~dtmp_c~ | ~double[elec_num][4][nucl_num][ncord + 1][ncord][walk_num]~ | vector of non-zero coefficients | |
| ~een_rescaled_e~ | ~double[walk_num][elec_num][elec_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord | |
| ~een_rescaled_e_date~ | ~uint64_t~ | Keep track of the date of creation | |
| ~een_rescaled_n~ | ~double[walk_num][elec_num][nucl_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord | |
| ~een_rescaled_n_date~ | ~uint64_t~ | Keep track of the date of creation | |
| ~een_rescaled_e_deriv_e~ | ~double[walk_num][elec_num][4][elec_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons | |
| ~een_rescaled_e_deriv_e_date~ | ~uint64_t~ | Keep track of the date of creation | |
| ~een_rescaled_n_deriv_e~ | ~double[walk_num][elec_num][4][nucl_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons | |
| ~een_rescaled_n_deriv_e_date~ | ~uint64_t~ | Keep track of the date of creation | |
| Variable | Type | In/Out | Description |
|-------------------------------+-------------------------------------------------------------+---------------------------------------------------------------------------------------------------------+-------------|
| ~dim_cord_vect~ | ~int64_t~ | Number of unique C coefficients | |
| ~dim_cord_vect_date~ | ~uint64_t~ | Number of unique C coefficients | |
| ~asymp_jasb~ | ~double[2]~ | Asymptotic component | |
| ~asymp_jasb_date~ | ~uint64_t~ | Asymptotic component | |
| ~cord_vect_full~ | ~double[dim_cord_vect][nucl_num]~ | vector of non-zero coefficients | |
| ~cord_vect_full_date~ | ~uint64_t~ | Keep track of changes here | |
| ~lkpm_combined_index~ | ~int64_t[4][dim_cord_vect]~ | Transform l,k,p, and m into consecutive indices | |
| ~lkpm_combined_index_date~ | ~uint64_t~ | Transform l,k,p, and m into consecutive indices | |
| ~tmp_c~ | ~double[elec_num][nucl_num][ncord + 1][ncord][walk_num]~ | vector of non-zero coefficients | |
| ~dtmp_c~ | ~double[elec_num][4][nucl_num][ncord + 1][ncord][walk_num]~ | vector of non-zero coefficients | |
| ~een_rescaled_e~ | ~double[walk_num][elec_num][elec_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord | |
| ~een_rescaled_e_date~ | ~uint64_t~ | Keep track of the date of creation | |
| ~een_rescaled_n~ | ~double[walk_num][elec_num][nucl_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord | |
| ~een_rescaled_n_date~ | ~uint64_t~ | Keep track of the date of creation | |
| ~een_rescaled_e_deriv_e~ | ~double[walk_num][elec_num][4][elec_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons | |
| ~een_rescaled_e_deriv_e_date~ | ~uint64_t~ | Keep track of the date of creation | |
| ~een_rescaled_n_deriv_e~ | ~double[walk_num][elec_num][4][nucl_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons | |
| ~een_rescaled_n_deriv_e_date~ | ~uint64_t~ | Keep track of the date of creation | |
For H2O we have the following data:
@ -1088,7 +1088,7 @@ assert(w == walk_num);
assert(qmckl_electron_provided(context));
rc = qmckl_set_electron_coord (context, 'N', elec_coord);
rc = qmckl_set_electron_coord (context, 'N', elec_coord, walk_num*3*elec_num);
assert(rc == QMCKL_SUCCESS);
double elec_coord2[walk_num*3*elec_num];

106
org/qmckl_local_energy.org
View File

@ -13,7 +13,7 @@ E_L = KE + PE
Where the kinetic energy is given as:
\[
KE = -\frac{1}{2}\frac{\bigtriangleup \Psi}{\Psi}
KE = -\frac{1}{2}\frac{\bigtriangleup \Psi}{\Psi}
\]
The laplacian of the wavefunction in the single-determinant
@ -59,7 +59,7 @@ F(x) = 2\frac{\nabla \Psi(r)}{\Psi(r)}
y = x + D F(x) \delta t + \chi
\]
Where \[\chi\] is a random number with gaussian distribution centered at 0
Where \[\chi\] is a random number with gaussian distribution centered at 0
and width of \[2D\delta t\]. Here \[D\] is the drift parameter.
3. Acceptance probability - \[min\left[1, q(y,x)\right]\]
@ -102,7 +102,7 @@ int main() {
qmckl_exit_code rc;
#+end_src
#+begin_src c :tangle (eval c)
#ifdef HAVE_CONFIG_H
#include "config.h"
@ -138,7 +138,7 @@ int main() {
| |
Computed data:
|--------------+---------------------------+----------------------------|
| ~e_kin~ | ~[walk_num]~ | total kinetic energy |
| ~e_pot~ | ~[walk_num]~ | total potential energy |
@ -147,9 +147,9 @@ int main() {
| ~y_move~ | ~[3][walk_num]~ | The diffusion move |
| ~accep_prob~ | ~[walk_num]~ | The acceptance probability |
|--------------+---------------------------+----------------------------|
** Data structure
#+begin_src c :comments org :tangle (eval h_private_type)
typedef struct qmckl_local_energy_struct {
double * e_kin;
@ -188,7 +188,7 @@ typedef struct qmckl_local_energy_struct {
Where the kinetic energy is given as:
\[
KE = -\frac{1}{2}\frac{\bigtriangleup \Psi}{\Psi}
KE = -\frac{1}{2}\frac{\bigtriangleup \Psi}{\Psi}
\]
The laplacian of the wavefunction in the single-determinant
@ -199,14 +199,14 @@ case is given as follows:
\]
*** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_kinetic_energy(qmckl_context context, double* const kinetic_energy);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_kinetic_energy(qmckl_context context, double * const kinetic_energy) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
@ -245,7 +245,7 @@ qmckl_exit_code qmckl_provide_kinetic_energy(qmckl_context context);
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_provide_kinetic_energy(qmckl_context context) {
qmckl_exit_code rc;
qmckl_exit_code rc;
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
@ -259,14 +259,14 @@ qmckl_exit_code qmckl_provide_kinetic_energy(qmckl_context context) {
"qmckl_electron",
NULL);
}
if(!(ctx->electron.provided)) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_electron",
NULL);
}
if (!ctx->ao_basis.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
@ -344,7 +344,7 @@ qmckl_exit_code qmckl_provide_kinetic_energy(qmckl_context context) {
QMCKL_FAILURE,
"compute_kinetic_energy",
"Not yet implemented");
}
}
if (rc != QMCKL_SUCCESS) {
return rc;
}
@ -380,7 +380,7 @@ qmckl_exit_code qmckl_provide_kinetic_energy(qmckl_context context) {
| ~double~ | ~det_inv_matrix_alpha[det_num_alpha][walk_num][alpha_num][alpha_num]~ | in | Value, gradients and Laplacian of the Det |
| ~double~ | ~det_inv_matrix_beta[det_num_beta][walk_num][beta_num][beta_num]~ | in | Value, gradients and Laplacian of the Det |
| ~double~ | ~e_kin[walk_num]~ | out | Kinetic energy |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_kinetic_energy_f(context, walk_num, &
det_num_alpha, det_num_beta, alpha_num, beta_num, elec_num, mo_index_alpha, mo_index_beta, &
@ -438,34 +438,20 @@ integer function qmckl_compute_kinetic_energy_f(context, walk_num, &
do idet = 1, det_num_alpha
do iwalk = 1, walk_num
! Alpha part
tmp_e = 0.0d0
do imo = 1, alpha_num
do ielec = 1, alpha_num
mo_id = mo_index_alpha(imo, iwalk, idet)
e_kin(iwalk) = e_kin(iwalk) - 0.5d0 * det_inv_matrix_alpha(imo, ielec, iwalk, idet) * &
mo_vgl(mo_id, ielec, 5)
!print *,"det alpha = ",det_inv_matrix_alpha(imo,ielec,iwalk,idet)
!print *,mo_vgl(mo_id,ielec,5)
!!print *," det val = ",det_value_alpha(iwalk,idet)
!tmp_e = tmp_e - 0.5d0 * det_inv_matrix_alpha(imo, ielec, iwalk, idet) * &
! mo_vgl(mo_id, ielec, 5)
end do
!print *,"e_kin = ",tmp_e
end do
! Beta part
tmp_e = 0.0d0
do imo = 1, beta_num
do ielec = 1, beta_num
mo_id = mo_index_beta(imo, iwalk, idet)
e_kin(iwalk) = e_kin(iwalk) - 0.5d0 * det_inv_matrix_beta(imo, ielec, iwalk, idet) * &
mo_vgl(mo_id, alpha_num + ielec, 5)
!print *,"det beta = ",det_inv_matrix_beta(imo,ielec,iwalk,idet)
!print *,mo_vgl(mo_id,alpha_num+ielec,5)
!!print *," det val = ",det_value_alpha(iwalk,idet)
!tmp_e = tmp_e - 0.5d0 * det_inv_matrix_beta(imo, ielec, iwalk, idet) * &
! mo_vgl(mo_id, alpha_num + ielec, 5)
end do
!print *,"e_kin = ",tmp_e
end do
end do
end do
@ -493,7 +479,7 @@ end function qmckl_compute_kinetic_energy_f
const double* det_value_beta,
const double* det_inv_matrix_alpha,
const double* det_inv_matrix_beta,
double* const e_kin );
double* const e_kin );
#+end_src
#+CALL: generate_c_interface(table=qmckl_compute_kinetic_energy_args,rettyp=get_value("CRetType"),fname="qmckl_compute_kinetic_energy"))
@ -560,7 +546,7 @@ end function qmckl_compute_kinetic_energy_f
end function qmckl_compute_kinetic_energy
#+end_src
*** Test
#+begin_src c :tangle (eval c_test) :exports none
@ -576,7 +562,7 @@ const int64_t nucl_num = chbrclf_nucl_num;
const double* nucl_charge = chbrclf_charge;
const double* nucl_coord = &(chbrclf_nucl_coord[0][0]);
rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num);
rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num);
assert (rc == QMCKL_SUCCESS);
rc = qmckl_set_electron_walk_num (context, walk_num);
@ -584,7 +570,7 @@ assert (rc == QMCKL_SUCCESS);
assert(qmckl_electron_provided(context));
rc = qmckl_set_electron_coord (context, 'N', elec_coord);
rc = qmckl_set_electron_coord (context, 'N', elec_coord, walk_num*elec_num*3);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_num (context, nucl_num);
@ -756,7 +742,7 @@ rc = qmckl_get_kinetic_energy(context, &(kinetic_energy[0]));
assert (rc == QMCKL_SUCCESS);
#+end_src
** Potential energy
:PROPERTIES:
:Name: qmckl_compute_potential_energy
@ -786,14 +772,14 @@ contributions.
\]
*** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_potential_energy(qmckl_context context, double* const potential_energy);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_potential_energy(qmckl_context context, double * const potential_energy) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
@ -839,21 +825,21 @@ qmckl_exit_code qmckl_provide_potential_energy(qmckl_context context) {
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
qmckl_exit_code rc;
qmckl_exit_code rc;
if(!(ctx->nucleus.provided)) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_electron",
NULL);
}
if(!(ctx->electron.provided)) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_electron",
NULL);
}
if (!ctx->ao_basis.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
@ -932,7 +918,7 @@ qmckl_exit_code qmckl_provide_potential_energy(qmckl_context context) {
QMCKL_FAILURE,
"compute_potential_energy",
"Not yet implemented");
}
}
if (rc != QMCKL_SUCCESS) {
return rc;
}
@ -960,7 +946,7 @@ qmckl_exit_code qmckl_provide_potential_energy(qmckl_context context) {
| ~double~ | ~en_pot[walk_num]~ | in | en potential |
| ~double~ | ~repulsion~ | in | en potential |
| ~double~ | ~e_pot[walk_num]~ | out | Potential energy |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_potential_energy_f(context, walk_num, &
elec_num, nucl_num, ee_pot, en_pot, repulsion, e_pot) &
@ -1014,7 +1000,7 @@ end function qmckl_compute_potential_energy_f
const double* ee_pot,
const double* en_pot,
const double repulsion,
double* const e_pot );
double* const e_pot );
#+end_src
#+CALL: generate_c_interface(table=qmckl_compute_potential_energy_args,rettyp=get_value("CRetType"),fname="qmckl_compute_potential_energy"))
@ -1043,7 +1029,7 @@ end function qmckl_compute_potential_energy_f
end function qmckl_compute_potential_energy
#+end_src
*** Test
#+begin_src c :tangle (eval c_test) :exports none
// Calculate the Potential energy
@ -1070,14 +1056,14 @@ 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);
#+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) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
@ -1129,14 +1115,14 @@ qmckl_exit_code qmckl_provide_local_energy(qmckl_context context) {
"qmckl_electron",
NULL);
}
if(!(ctx->electron.provided)) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_electron",
NULL);
}
if (!ctx->ao_basis.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
@ -1205,7 +1191,7 @@ qmckl_exit_code qmckl_provide_local_energy(qmckl_context context) {
QMCKL_FAILURE,
"compute_local_energy",
"Not yet implemented");
}
}
if (rc != QMCKL_SUCCESS) {
return rc;
}
@ -1230,7 +1216,7 @@ qmckl_exit_code qmckl_provide_local_energy(qmckl_context context) {
| ~double~ | ~e_kin[walk_num]~ | in | e kinetic |
| ~double~ | ~e_pot[walk_num]~ | in | e potential |
| ~double~ | ~e_local[walk_num]~ | out | local energy |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_local_energy_f(context, walk_num, &
e_kin, e_pot, e_local) &
@ -1273,7 +1259,7 @@ end function qmckl_compute_local_energy_f
const int64_t walk_num,
const double* e_kin,
const double* e_pot,
double* const e_local );
double* const e_local );
#+end_src
#+CALL: generate_c_interface(table=qmckl_compute_local_energy_args,rettyp=get_value("CRetType"),fname="qmckl_compute_local_energy"))
@ -1299,7 +1285,7 @@ end function qmckl_compute_local_energy_f
end function qmckl_compute_local_energy
#+end_src
*** Test
#+begin_src c :tangle (eval c_test) :exports none
// Calculate the Local energy
@ -1318,7 +1304,7 @@ assert (rc == QMCKL_SUCCESS);
:FRetType: qmckl_exit_code
:END:
The drift vector is calculated as the ration of the gradient
The drift vector is calculated as the ration of the gradient
with the determinant of the wavefunction.
\[
@ -1326,14 +1312,14 @@ with the determinant of the wavefunction.
\]
*** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_drift_vector(qmckl_context context, double* const drift_vector);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_drift_vector(qmckl_context context, double * const drift_vector) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
@ -1385,14 +1371,14 @@ qmckl_exit_code qmckl_provide_drift_vector(qmckl_context context) {
"qmckl_electron",
NULL);
}
if(!(ctx->electron.provided)) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_electron",
NULL);
}
if (!ctx->ao_basis.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
@ -1433,7 +1419,7 @@ qmckl_exit_code qmckl_provide_drift_vector(qmckl_context context) {
ctx->local_energy.r_drift = r_drift;
}
qmckl_exit_code rc;
qmckl_exit_code rc;
if (ctx->det.type == 'G') {
rc = qmckl_compute_drift_vector(context,
ctx->det.walk_num,
@ -1454,7 +1440,7 @@ qmckl_exit_code qmckl_provide_drift_vector(qmckl_context context) {
QMCKL_FAILURE,
"compute_drift_vector",
"Not yet implemented");
}
}
if (rc != QMCKL_SUCCESS) {
return rc;
}
@ -1488,7 +1474,7 @@ qmckl_exit_code qmckl_provide_drift_vector(qmckl_context context) {
| ~double~ | ~det_inv_matrix_alpha[det_num_alpha][walk_num][alpha_num][alpha_num]~ | in | Value, gradients and Laplacian of the Det |
| ~double~ | ~det_inv_matrix_beta[det_num_beta][walk_num][beta_num][beta_num]~ | in | Value, gradients and Laplacian of the Det |
| ~double~ | ~r_drift[walk_num][elec_num][3]~ | out | Kinetic energy |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_drift_vector_f(context, walk_num, &
det_num_alpha, det_num_beta, alpha_num, beta_num, elec_num, mo_index_alpha, mo_index_beta, &
@ -1593,7 +1579,7 @@ end function qmckl_compute_drift_vector_f
const double* mo_vgl,
const double* det_inv_matrix_alpha,
const double* det_inv_matrix_beta,
double* const r_drift );
double* const r_drift );
#+end_src
#+CALL: generate_c_interface(table=qmckl_compute_drift_vector_args,rettyp=get_value("CRetType"),fname="qmckl_compute_drift_vector"))
@ -1654,7 +1640,7 @@ end function qmckl_compute_drift_vector_f
end function qmckl_compute_drift_vector
#+end_src
*** Test
#+begin_src c :tangle (eval c_test) :exports none
// Calculate the Drift vector

50
org/qmckl_mo.org
View File

@ -2,7 +2,7 @@
#+SETUPFILE: ../tools/theme.setup
#+INCLUDE: ../tools/lib.org
The molecular orbitals (MOs) are defined in the basis of AOs along with a AO to MO
The molecular orbitals (MOs) are defined in the basis of AOs along with a AO to MO
coefficient matrix \[C\]. Using these coefficients (e.g. from Hartree Fock SCF method)
the MOs are defined as follows:
@ -89,7 +89,7 @@ int main() {
| ~coefficient~ | ~[mo_num][ao_num]~ | Orbital coefficients |
Computed data:
|---------------+-------------------------+----------------------------------------------------------------------------------------|
|---------------+-------------------------+----------------------------------------------------------------------------------------|
| ~mo_vgl~ | ~[5][elec_num][mo_num]~ | Value, gradients, Laplacian of the MOs at electron positions |
@ -100,7 +100,7 @@ int main() {
#+begin_src c :comments org :tangle (eval h_private_type)
typedef struct qmckl_mo_basis_struct {
int64_t mo_num;
int64_t mo_num;
double * coefficient;
double * mo_vgl;
@ -118,10 +118,10 @@ typedef struct qmckl_mo_basis_struct {
Some values are initialized by default, and are not concerned by
this mechanism.
#+begin_src c :comments org :tangle (eval h_private_func)
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_exit_code qmckl_init_mo_basis(qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c)
qmckl_exit_code qmckl_init_mo_basis(qmckl_context context) {
@ -137,7 +137,7 @@ qmckl_exit_code qmckl_init_mo_basis(qmckl_context context) {
return QMCKL_SUCCESS;
}
#+end_src
** Access functions
#+begin_src c :comments org :tangle (eval h_func) :exports none
@ -344,6 +344,7 @@ qmckl_exit_code qmckl_finalize_mo_basis(qmckl_context context);
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_finalize_mo_basis(qmckl_context context) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT,
@ -354,14 +355,7 @@ qmckl_exit_code qmckl_finalize_mo_basis(qmckl_context context) {
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
qmckl_exit_code rc;
rc = qmckl_provide_mo_vgl(context);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_finalize_mo_basis",
NULL);
}
qmckl_exit_code rc = QMCKL_SUCCESS;
return rc;
}
#+end_src
@ -378,7 +372,7 @@ qmckl_exit_code qmckl_get_mo_basis_vgl(qmckl_context context, double* const mo_v
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_mo_basis_vgl(qmckl_context context, double* const mo_vgl) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
@ -425,7 +419,7 @@ qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context);
qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context)
{
qmckl_exit_code rc;
qmckl_exit_code rc;
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
@ -454,7 +448,7 @@ qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context)
"qmckl_electron",
NULL);
}
if (!ctx->mo_basis.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
@ -498,7 +492,7 @@ qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context)
return QMCKL_SUCCESS;
}
#+end_src
*** Compute
:PROPERTIES:
:Name: qmckl_compute_mo_basis_vgl
@ -515,7 +509,7 @@ qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context)
| ~double~ | ~ao_vgl[5][elec_num][ao_num]~ | in | Value, gradients and Laplacian of the AOs |
| ~double~ | ~mo_vgl[5][elec_num][mo_num]~ | out | Value, gradients and Laplacian of the MOs |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_mo_basis_vgl_f(context, &
ao_num, mo_num, elec_num, &
@ -553,7 +547,7 @@ integer function qmckl_compute_mo_basis_vgl_f(context, &
info = QMCKL_SUCCESS
info_qmckl_dgemm_value = QMCKL_SUCCESS
! Don't compute exponentials when the result will be almost zero.
! TODO : Use numerical precision here
cutoff = -dlog(1.d-15)
@ -577,7 +571,7 @@ integer function qmckl_compute_mo_basis_vgl_f(context, &
!end do
ao_vgl_big = reshape(ao_vgl(:, :, :),(/ao_num, elec_num*5_8/))
LDB = size(ao_vgl_big,1)
LDB = size(ao_vgl_big,1)
LDC = size(mo_vgl_big,1)
info = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, &
@ -604,7 +598,7 @@ end function qmckl_compute_mo_basis_vgl_f
const int64_t elec_num,
const double* coef_normalized,
const double* ao_vgl,
double* const mo_vgl );
double* const mo_vgl );
#+end_src
@ -640,7 +634,7 @@ end function qmckl_compute_mo_basis_vgl_f
import numpy as np
def f(a,x,y):
return np.sum( [c * np.exp( -b*(np.linalg.norm(x-y))**2) for b,c in a] )
return np.sum( [c * np.exp( -b*(np.linalg.norm(x-y))**2) for b,c in a] )
def df(a,x,y,n):
h0 = 1.e-6
@ -712,7 +706,7 @@ const int64_t nucl_num = chbrclf_nucl_num;
const double* nucl_charge = chbrclf_charge;
const double* nucl_coord = &(chbrclf_nucl_coord[0][0]);
rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num);
rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num);
assert (rc == QMCKL_SUCCESS);
rc = qmckl_set_electron_walk_num (context, walk_num);
@ -720,7 +714,7 @@ assert (rc == QMCKL_SUCCESS);
assert(qmckl_electron_provided(context));
rc = qmckl_set_electron_coord (context, 'N', elec_coord);
rc = qmckl_set_electron_coord (context, 'N', elec_coord, walk_num*elec_num*3);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_num (context, nucl_num);
@ -854,7 +848,7 @@ assert (rc == QMCKL_SUCCESS);
// elec_coord[0] = point_x[i];
// elec_coord[1] = point_y[j];
// elec_coord[2] = point_z[k];
// rc = qmckl_set_electron_coord (context, 'N', elec_coord);
// rc = qmckl_set_electron_coord (context, 'N', elec_coord);
// assert(rc == QMCKL_SUCCESS);
//
// // Calculate value of MO (1st electron)
@ -882,8 +876,8 @@ printf(" mo_vgl mo_vgl[1][26][223] %25.15e\n", mo_vgl[1][2][3]);
printf(" mo_vgl mo_vgl[0][26][224] %25.15e\n", mo_vgl[0][2][3]);
printf(" mo_vgl mo_vgl[1][26][224] %25.15e\n", mo_vgl[1][2][3]);
printf("\n");
}
}
#+end_src
* End of files :noexport:

590
org/qmckl_point.org Normal file
View File

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

23
org/table_of_contents
View File

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