QMCkl source code documentation

1. Introduction
2. Documentation
3. Acknowledgments

1 Introduction

The ultimate goal of the QMCkl library is to provide a high-performance implementation of the main kernels of QMC. In this particular implementation of the library, we focus on the definition of the API and the tests, and on a pedagogical presentation of the algorithms. We expect the HPC experts to use this repository as a reference for re-writing optimized libraries.

1.1 Literate programming

In a traditional source code, most of the lines of source files of a program are code, scripts, Makefiles, and only a few lines are comments explaining parts of the code that are non-trivial to understand. The documentation of the prorgam is usually written in a separate directory, and is often outdated compared to the code.

Literate programming is a different approach to programming, where the program is considered as a publishable-quality document. Most of the lines of the source files are text, mathematical formulas, tables, figures, etc, and the lines of code are just the translation in a computer language of the ideas and algorithms expressed in the text. More importantly, the "document" is structured like a text document with sections, subsections, a bibliography, a table of contents etc, and the place where pieces of code appear are the places where they should belong for the reader to understand the logic of the program, not the places where the compiler expects to find them. Both the publishable-quality document and the binary executable are produced from the same source files.

Literate programming is particularly well adapted in this context, as the central part of this project is the documentation of an API. The implementation of the algorithms is just an expression of the algorithms in a language that can be compiled, so that the correctness of the algorithms can be tested.

We have chosen to write the source files in org-mode format, as any text editor can be used to edit org-mode files. To produce the documentation, there exists multiple possibilities to convert org-mode files into different formats such as HTML or PDF. The source code is easily extracted from the org-mode files invoking the Emacs text editor from the command-line in the Makefile, and then the produced files are compiled. Moreover, within the Emacs text editor the source code blocks can be executed interactively, in the same spirit as Jupyter notebooks.

1.2 Source code editing

For a tutorial on literate programming with org-mode, follow this link.

Any text editor can be used to edit org-mode files. For a better user experience Emacs is recommended. For users hating Emacs, it is good to know that Emacs can behave like Vim when switched into ``Evil'' mode.

In the tools/init.el file, we provide a minimal Emacs configuration file for vim users. This file should be copied into .emacs.d/init.el.

For users with a preference for Jupyter notebooks, we also provide the tools/nb_to_org.sh script can convert jupyter notebooks into org-mode files.

Note that pandoc can be used to convert multiple markdown formats into org-mode.

1.3 Choice of the programming language

Most of the codes of the 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.

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.

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.

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. Fortran interfaces should also be written in the qmckl_f.f90 file.

For more guidelines on using Fortran to generate a C interface, see this link.

1.4 Design of the library

The proposed API should allow the library to: deal with memory transfers between CPU and accelerators, and to use different levels of floating-point precision. We chose a multi-layered design with low-level and high-level functions (see below).

1.4.1 Naming conventions

To avoid namespace collisions, we use qmckl_ as a prefix for all exported functions and variables. All exported header files should have a file name prefixed with qmckl_.

If the name of the org-mode file is xxx.org, the name of the produced C files should be xxx.c and xxx.h and the name of the produced Fortran 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.

1.4.2 Application programming interface

In the C language, the number of bits used by the integer types can change from one architecture to another one. To circumvent this problem, we choose to use the integer types defined in <stdint.h> where the number of bits used for the integers are fixed.

To ensure that the library will be easily usable in any other language than C, we restrict the data types in the interfaces to the following:

32-bit and 64-bit integers, scalars and and arrays (int32_t and int64_t)
32-bit and 64-bit floats, scalars and and arrays (float and double)
Pointers are always casted into 64-bit integers, even on legacy 32-bit architectures
ASCII strings are represented as a pointers to character arrays and terminated by a '\0' character (C convention).
Complex numbers can be represented by an array of 2 floats.
Boolean variables are stored as integers, 1 for true and 0 for false
Floating point variables should be by default
double unless explicitly mentioned
integers used for counting should always be int64_t

To facilitate the use in other languages than C, we will provide some bindings in other languages in other repositories.

1.4.3 Global state

Global variables should be avoided in the library, because it is possible that one single program needs to use multiple instances of the library. To solve this problem we propose to use a pointer to a context variable, built by the library with the qmckl_context_create function. The context contains the global state of the library, and is used as the first argument of many QMCkl functions.

The internal structure of the context is not specified, to give a maximum of freedom to the different implementations. Modifying the state is done by setters and getters, prefixed by qmckl_context_set_ an qmckl_context_get_. When a context variable is modified by a setter, a copy of the old data structure is made and updated, and the pointer to the new data structure is returned, such that the old contexts can still be accessed. It is also possible to modify the state in an impure fashion, using the qmckl_context_update_ functions. The context and its old versions can be destroyed with qmckl_context_destroy.

1.4.4 Low-level functions

Low-level functions are very simple functions which are leaves of the function call tree (they don't call any other QMCkl function).

These functions are pure, and unaware of the QMCkl context. They are not allowed to allocate/deallocate memory, and if they need temporary memory it should be provided in input.

1.4.5 High-level functions

High-level functions are at the top of the function call tree. They are able to choose which lower-level function to call depending on the required precision, and do the corresponding type conversions. These functions are also responsible for allocating temporary storage, to simplify the use of accelerators.

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.

1.4.6 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.

1.5 Algorithms

Reducing the scaling of an algorithm usually implies also reducing its arithmetic complexity (number of flops per byte). Therefore, for small sizes \(\mathcal{O}(N^3)\) and \(\mathcal{O}(N^2)\) algorithms are better adapted than linear scaling algorithms. As QMCkl is a general purpose library, multiple algorithms should be implemented adapted to different problem sizes.

1.6 Rules for the API

stdint should be used for integers (int32_t, int64_t)
integers used for counting should always be int64_t
floats should be by default double, unless explicitly mentioned
pointers are converted to int64_t to increase portability

2 Documentation

2.1 `qmckl.h` header file

The qmckl.h header file has to be included in C codes when QMCkl functions are used: #+BEGIN_SRC C :tangle none #include "qmckl.h" #+END_SRC f90

In Fortran programs, the qmckl_f.f90 interface file should be included in the source code using the library, and the Fortran codes should use the qmckl module as #+BEGIN_SRC f90 :tangle none use qmckl #+END_SRC f90

2.2 Context

This file is written in C because it is more natural to express the context in C than in Fortran.

2 files are produced:

a source file : qmckl_context.c
a test file : test_qmckl_context.c

2.2.1 Context

The context variable is a handle for the state of the library, and is stored in the following data structure, which can't be seen outside of the library. To simplify compatibility with other languages, the pointer to the internal data structure is converted into a 64-bit signed integer, defined in the qmckl_context type. A value of 0 for the context is equivalent to a NULL pointer.

typedef int64_t qmckl_context ;

Data for error handling

We define here the the data structure containing the strings necessary for error handling.

#define  QMCKL_MAX_FUN_LEN  256
#define  QMCKL_MAX_MSG_LEN  1024

typedef struct qmckl_error_struct {

  qmckl_exit_code exit_code;
  char function[QMCKL_MAX_FUN_LEN];
  char message [QMCKL_MAX_MSG_LEN];

} qmckl_error_struct;

Basis set data structure

Data structure for the info related to the atomic orbitals basis set.

typedef struct qmckl_ao_basis_struct {

  int64_t   shell_num;
  int64_t   prim_num;
  int64_t * shell_center;
  int32_t * shell_ang_mom;
  double  * shell_factor;
  double  * exponent    ;
  double  * coefficient ;
  int64_t * shell_prim_num;
  char      type;

} qmckl_ao_basis_struct;

Source

The tag is used internally to check if the memory domain pointed by a pointer is a valid context.

typedef struct qmckl_context_struct {

  struct qmckl_context_struct * prev;

  /* Molecular system */
  //  struct qmckl_nucleus_struct     * nucleus;
  //  struct qmckl_electron_struct    * electron;
  struct qmckl_ao_basis_struct      * ao_basis;
  //  struct qmckl_mo_struct          * mo;
  //  struct qmckl_determinant_struct * det;

  /* Numerical precision */
  uint32_t tag;
  int32_t precision;
  int32_t range;

  /* Error handling */
  struct qmckl_error_struct   * error;

} qmckl_context_struct;

#define VALID_TAG   0xBEEFFACE
#define INVALID_TAG 0xDEADBEEF

qmckl_context_update_error

qmckl_exit_code
qmckl_context_update_error(qmckl_context context, const qmckl_exit_code exit_code, const char* function, const char* message);

Source

qmckl_exit_code
qmckl_context_update_error(qmckl_context context, const qmckl_exit_code exit_code, const char* function, const char* message)
{
  assert (context != 0);
  assert (function != NULL);
  assert (message != NULL);
  assert (exit_code > 0);
  assert (exit_code < QMCKL_INVALID_EXIT_CODE);

  qmckl_context_struct* ctx = (qmckl_context_struct*) context;
  if (ctx == NULL) return QMCKL_FAILURE;

  if (ctx->error != NULL) {
    free(ctx->error);
    ctx->error = NULL;
  }

  qmckl_error_struct* error = (qmckl_error_struct*) qmckl_malloc (context, sizeof(qmckl_error_struct));
  error->exit_code = exit_code;
  strcpy(error->function, function);
  strcpy(error->message, message);

  ctx->error = error;

  return QMCKL_SUCCESS;
}

TODO Test

qmckl_context_set_error

qmckl_context
qmckl_context_set_error(qmckl_context context, const qmckl_exit_code exit_code, const char* function, const char* message);

Source

qmckl_context
qmckl_context_set_error(qmckl_context context, const qmckl_exit_code exit_code, const char* function, const char* message)
{
  assert (context != 0);
  assert (function != NULL);
  assert (message != NULL);
  assert (exit_code > 0);
  assert (exit_code < QMCKL_INVALID_EXIT_CODE);

  qmckl_context new_context = qmckl_context_copy(context);
  if (new_context == 0) return context;

  if (qmckl_context_update_error(new_context, exit_code,
                                 function, message) != QMCKL_SUCCESS) {
    return context;
  }

  return new_context;
}

TODO Test

qmckl_context_check

Checks if the domain pointed by the pointer is a valid context. Returns the input qmckl_context if the context is valid, 0 otherwise.

qmckl_context qmckl_context_check(const qmckl_context context) ;

Source

qmckl_context qmckl_context_check(const qmckl_context context) {

  if (context == (qmckl_context) 0) return (qmckl_context) 0;

  const qmckl_context_struct * ctx = (qmckl_context_struct*) context;

  if (ctx->tag != VALID_TAG) return (qmckl_context) 0;

  return context;
}

qmckl_context_create

To create a new context, use qmckl_context_create().

On success, returns a pointer to a context using the qmckl_context type
Returns 0 upon failure to allocate the internal data structure
```
qmckl_context qmckl_context_create();
```

Source

qmckl_context qmckl_context_create() {

  qmckl_context_struct* context =
    (qmckl_context_struct*) qmckl_malloc ((qmckl_context) 0, sizeof(qmckl_context_struct));
  if (context == NULL) {
    return (qmckl_context) 0;
  }

  context->prev      = NULL;
  context->ao_basis  = NULL;
  context->precision = QMCKL_DEFAULT_PRECISION;
  context->range     = QMCKL_DEFAULT_RANGE;
  context->tag       = VALID_TAG;
  context->error     = NULL;

  return (qmckl_context) context;
}

Fortran interface

interface
   integer (c_int64_t) function qmckl_context_create() bind(C)
     use, intrinsic :: iso_c_binding
   end function qmckl_context_create
end interface

qmckl_context_copy

This function makes a shallow copy of the current context.

Copying the 0-valued context returns 0
On success, returns a pointer to the new context using the qmckl_context type
Returns 0 upon failure to allocate the internal data structure for the new context
```
qmckl_context qmckl_context_copy(const qmckl_context context);
```

Source

qmckl_context qmckl_context_copy(const qmckl_context context) {

  const qmckl_context checked_context = qmckl_context_check(context);

  if (checked_context == (qmckl_context) 0) {
    return (qmckl_context) 0;
  }

  qmckl_context_struct* old_context = (qmckl_context_struct*) checked_context;

  qmckl_context_struct* new_context =
    (qmckl_context_struct*) qmckl_malloc (context, sizeof(qmckl_context_struct));
  if (new_context == NULL) {
    return (qmckl_context) 0;
  }

  new_context->prev      = old_context;
  new_context->ao_basis  = old_context->ao_basis;
  new_context->precision = old_context->precision;
  new_context->range     = old_context->range;
  new_context->tag       = VALID_TAG;
  new_context->error     = old_context->error;

  return (qmckl_context) new_context;
}

Fortran interface

interface
   integer (c_int64_t) function qmckl_context_copy(context) bind(C)
     use, intrinsic :: iso_c_binding
     integer (c_int64_t), intent(in), value :: context
   end function qmckl_context_copy
end interface

qmckl_context_previous

Returns the previous context

On success, returns the ancestor of the current context
Returns 0 for the initial context

Returns 0 for the 0-valued context

qmckl_context qmckl_context_previous(const qmckl_context context);

Source

qmckl_context qmckl_context_previous(const qmckl_context context) {

  const qmckl_context checked_context = qmckl_context_check(context);
  if (checked_context == (qmckl_context) 0) {
    return (qmckl_context) 0;
  }

  const qmckl_context_struct* ctx = (qmckl_context_struct*) checked_context;
  return qmckl_context_check((qmckl_context) ctx->prev);
}

Fortran interface

interface
   integer (c_int64_t) function qmckl_context_previous(context) bind(C)
     use, intrinsic :: iso_c_binding
     integer (c_int64_t), intent(in), value :: context
   end function qmckl_context_previous
end interface

qmckl_context_destroy

Destroys the current context, leaving the ancestors untouched.

Succeeds if the current context is properly destroyed
Fails otherwise
Fails if the 0-valued context is given in argument
Fails if the the pointer is not a valid context

qmckl_exit_code qmckl_context_destroy(qmckl_context context);

Source

qmckl_exit_code qmckl_context_destroy(const qmckl_context context) {

  const qmckl_context checked_context = qmckl_context_check(context);
  if (checked_context == (qmckl_context) 0) return QMCKL_FAILURE;

  qmckl_context_struct* ctx = (qmckl_context_struct*) context;
  if (ctx == NULL) return QMCKL_FAILURE;

  ctx->tag = INVALID_TAG;
  return qmckl_free(context,ctx);
}

Fortran interface

interface
   integer (c_int32_t) function qmckl_context_destroy(context) bind(C)
     use, intrinsic :: iso_c_binding
     integer (c_int64_t), intent(in), value :: context
   end function qmckl_context_destroy
end interface

Basis set

For H₂ with the following basis set,

HYDROGEN
S   5
1         3.387000E+01           6.068000E-03
2         5.095000E+00           4.530800E-02
3         1.159000E+00           2.028220E-01
4         3.258000E-01           5.039030E-01
5         1.027000E-01           3.834210E-01
S   1
1         3.258000E-01           1.000000E+00
S   1
1         1.027000E-01           1.000000E+00
P   1
1         1.407000E+00           1.000000E+00
P   1
1         3.880000E-01           1.000000E+00
D   1
1         1.057000E+00           1.0000000

we have:

type = 'G'
shell_num = 12
prim_num = 20
SHELL_CENTER = [1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2]
SHELL_ANG_MOM = ['S', 'S', 'S', 'P', 'P', 'D', 'S', 'S', 'S', 'P', 'P', 'D']
SHELL_PRIM_NUM = [5, 1, 1, 1, 1, 1, 5, 1, 1, 1, 1, 1]
prim_index = [1, 6, 7, 8, 9, 10, 11, 16, 17, 18, 19, 20]
EXPONENT = [ 33.87, 5.095, 1.159, 0.3258, 0.1027, 0.3258, 0.1027,
             1.407, 0.388, 1.057, 33.87, 5.095, 1.159, 0.3258, 0.1027,
             0.3258, 0.1027, 1.407, 0.388, 1.057]
COEFFICIENT = [ 0.006068, 0.045308, 0.202822, 0.503903, 0.383421,
                1.0, 1.0, 1.0, 1.0, 1.0, 0.006068, 0.045308, 0.202822,
                0.503903, 0.383421, 1.0, 1.0, 1.0, 1.0, 1.0]

qmckl_context_update_ao_basis

Updates the data describing the AO basis set into the context.

`type`	Gaussian or Slater
`shell_num`	Number of shells
`prim_num`	Total number of primitives
`SHELL_CENTER(shell_num)`	Id of the nucleus on which the shell is centered
`SHELL_ANG_MOM(shell_num)`	Id of the nucleus on which the shell is centered
`SHELL_FACTOR(shell_num)`	Normalization factor for the shell
`SHELL_PRIM_NUM(shell_num)`	Number of primitives in the shell
`SHELL_PRIM_INDEX(shell_num)`	Address of the first primitive of the shelll in the `EXPONENT` array
`EXPONENT(prim_num)`	Array of exponents
`COEFFICIENT(prim_num)`	Array of coefficients

qmckl_exit_code
qmckl_context_update_ao_basis(qmckl_context   context     , const char      type,
                              const int64_t   shell_num   , const int64_t   prim_num,
                              const int64_t * SHELL_CENTER, const int32_t * SHELL_ANG_MOM,
                              const double  * SHELL_FACTOR, const int64_t * SHELL_PRIM_NUM,
                              const int64_t * SHELL_PRIM_INDEX,
                              const double  * EXPONENT    , const double  * COEFFICIENT);

Source

qmckl_exit_code
qmckl_context_update_ao_basis(qmckl_context   context     , const char      type,
                              const int64_t   shell_num   , const int64_t   prim_num,
                              const int64_t * SHELL_CENTER, const int32_t * SHELL_ANG_MOM,
                              const double  * SHELL_FACTOR, const int64_t * SHELL_PRIM_NUM,
                              const int64_t * SHELL_PRIM_INDEX,
                              const double  * EXPONENT    , const double  * COEFFICIENT)
{

  int64_t i;

  /* Check input */

  if (type != 'G' && type != 'S') return QMCKL_FAILURE;
  if (shell_num <= 0) return QMCKL_FAILURE;
  if (prim_num  <= 0) return QMCKL_FAILURE;
  if (prim_num  < shell_num) return QMCKL_FAILURE;

  for (i=0 ; i<shell_num ; i++) {
    if (SHELL_CENTER[i] <= 0) return QMCKL_FAILURE;
    if (SHELL_PRIM_NUM[i] <= 0) return QMCKL_FAILURE;
    if (SHELL_ANG_MOM[i] < 0) return QMCKL_FAILURE;
    if (SHELL_PRIM_INDEX[i] < 0) return QMCKL_FAILURE;
  }

  for (i=0 ; i<prim_num ; i++) {
    if (EXPONENT[i] <= 0) return QMCKL_FAILURE;
  }

  qmckl_context_struct* ctx = (qmckl_context_struct*) context;
  if (ctx == NULL) return QMCKL_FAILURE;

  qmckl_ao_basis_struct* basis = (qmckl_ao_basis_struct*) malloc (sizeof(qmckl_ao_basis_struct));
  if (basis == NULL) return QMCKL_FAILURE;


  /* Memory allocations */

  basis->shell_center  = (int64_t*) malloc (shell_num * sizeof(int64_t));
  if (basis->shell_center == NULL) {
    qmckl_free(context, basis);
    return QMCKL_FAILURE;
  }

  basis->shell_ang_mom = (int32_t*) malloc (shell_num * sizeof(int32_t));
  if (basis->shell_ang_mom == NULL) {
    qmckl_free(context, basis->shell_center);
    qmckl_free(context, basis);
    return QMCKL_FAILURE;
  }

  basis->shell_prim_num= (int64_t*) malloc (shell_num * sizeof(int64_t));
  if (basis->shell_prim_num == NULL) {
    qmckl_free(context, basis->shell_ang_mom);
    qmckl_free(context, basis->shell_center);
    qmckl_free(context, basis);
    return QMCKL_FAILURE;
  }

  basis->shell_factor  = (double *) malloc (shell_num * sizeof(double ));
  if (basis->shell_factor == NULL) {
    qmckl_free(context, basis->shell_prim_num);
    qmckl_free(context, basis->shell_ang_mom);
    qmckl_free(context, basis->shell_center);
    qmckl_free(context, basis);
    return QMCKL_FAILURE;
  }

  basis->exponent      = (double *) malloc (prim_num  * sizeof(double ));
  if (basis->exponent == NULL) {
    qmckl_free(context, basis->shell_factor);
    qmckl_free(context, basis->shell_prim_num);
    qmckl_free(context, basis->shell_ang_mom);
    qmckl_free(context, basis->shell_center);
    qmckl_free(context, basis);
    return QMCKL_FAILURE;
  }

  basis->coefficient   = (double *) malloc (prim_num  * sizeof(double ));
  if (basis->coefficient == NULL) {
    qmckl_free(context, basis->exponent);
    qmckl_free(context, basis->shell_factor);
    qmckl_free(context, basis->shell_prim_num);
    qmckl_free(context, basis->shell_ang_mom);
    qmckl_free(context, basis->shell_center);
    qmckl_free(context, basis);
    return QMCKL_FAILURE;
  }


  /* Assign data */

  basis->type      = type;
  basis->shell_num = shell_num;
  basis->prim_num  = prim_num;

  for (i=0 ; i<shell_num ; i++) {
    basis->shell_center  [i] = SHELL_CENTER  [i];
    basis->shell_ang_mom [i] = SHELL_ANG_MOM [i];
    basis->shell_prim_num[i] = SHELL_PRIM_NUM[i];
    basis->shell_factor  [i] = SHELL_FACTOR  [i];
  }

  for (i=0 ; i<prim_num ; i++) {
    basis->exponent   [i] = EXPONENT[i];
    basis->coefficient[i] = COEFFICIENT[i];
  }

  ctx->ao_basis = basis;
  return QMCKL_SUCCESS;
}

Fortran interface

interface
   integer (c_int32_t) function qmckl_context_update_ao_basis(context, &
        typ, shell_num, prim_num, SHELL_CENTER, SHELL_ANG_MOM, SHELL_FACTOR, &
        SHELL_PRIM_NUM, SHELL_PRIM_INDEX, EXPONENT, COEFFICIENT) bind(C)
     use, intrinsic :: iso_c_binding
     integer (c_int64_t), intent(in), value :: context
     character(c_char)  , intent(in), value :: typ
     integer (c_int64_t), intent(in), value :: shell_num
     integer (c_int64_t), intent(in), value :: prim_num
     integer (c_int64_t), intent(in)        :: SHELL_CENTER(shell_num)
     integer (c_int32_t), intent(in)        :: SHELL_ANG_MOM(shell_num)
     double precision   , intent(in)        :: SHELL_FACTOR(shell_num)
     integer (c_int64_t), intent(in)        :: SHELL_PRIM_NUM(shell_num)
     integer (c_int64_t), intent(in)        :: SHELL_PRIM_INDEX(shell_num)
     double precision   , intent(in)        :: EXPONENT(prim_num)
     double precision   , intent(in)        :: COEFFICIENT(prim_num)
   end function qmckl_context_update_ao_basis
end interface

TODO Test

qmckl_context_set_ao_basis

Sets the data describing the AO basis set into the context.

`type`	Gaussian or Slater
`shell_num`	Number of shells
`prim_num`	Total number of primitives
`SHELL_CENTER(shell_num)`	Id of the nucleus on which the shell is centered
`SHELL_ANG_MOM(shell_num)`	Id of the nucleus on which the shell is centered
`SHELL_FACTOR(shell_num)`	Normalization factor for the shell
`SHELL_PRIM_NUM(shell_num)`	Number of primitives in the shell
`SHELL_PRIM_INDEX(shell_num)`	Address of the first primitive of the shelll in the `EXPONENT` array
`EXPONENT(prim_num)`	Array of exponents
`COEFFICIENT(prim_num)`	Array of coefficients

qmckl_context
qmckl_context_set_ao_basis(const qmckl_context context     , const char      type,
                           const int64_t       shell_num   , const int64_t   prim_num,
                           const int64_t *     SHELL_CENTER, const int32_t * SHELL_ANG_MOM,
                           const double  *     SHELL_FACTOR, const int64_t * SHELL_PRIM_NUM,
                           const int64_t *     SHELL_PRIM_INDEX,
                           const double  *     EXPONENT    , const double  * COEFFICIENT);

Source

qmckl_context
qmckl_context_set_ao_basis(const qmckl_context context     , const char      type,
                           const int64_t       shell_num   , const int64_t   prim_num,
                           const int64_t *     SHELL_CENTER, const int32_t * SHELL_ANG_MOM,
                           const double  *     SHELL_FACTOR, const int64_t * SHELL_PRIM_NUM,
                           const int64_t *     SHELL_PRIM_INDEX,
                           const double  *     EXPONENT    , const double  * COEFFICIENT)
{

  qmckl_context new_context = qmckl_context_copy(context);
  if (new_context == 0) return 0;

  if (qmckl_context_update_ao_basis(new_context, type, shell_num, prim_num,
                                    SHELL_CENTER, SHELL_ANG_MOM, SHELL_FACTOR,
                                    SHELL_PRIM_NUM, SHELL_PRIM_INDEX, EXPONENT,
                                    COEFFICIENT
                                    ) == QMCKL_FAILURE)
    return 0;

  return new_context;
}

Fortran interface

interface
   integer (c_int64_t) function qmckl_context_set_ao_basis(context, &
        typ, shell_num, prim_num, SHELL_CENTER, SHELL_ANG_MOM, SHELL_FACTOR, &
        SHELL_PRIM_NUM, SHELL_PRIM_INDEX, EXPONENT, COEFFICIENT) bind(C)
     use, intrinsic :: iso_c_binding
     integer (c_int64_t), intent(in), value :: context
     character(c_char)  , intent(in), value :: typ
     integer (c_int64_t), intent(in), value :: shell_num
     integer (c_int64_t), intent(in), value :: prim_num
     integer (c_int64_t), intent(in)        :: SHELL_CENTER(shell_num)
     integer (c_int32_t), intent(in)        :: SHELL_ANG_MOM(shell_num)
     double precision   , intent(in)        :: SHELL_FACTOR(shell_num)
     integer (c_int64_t), intent(in)        :: SHELL_PRIM_NUM(shell_num)
     integer (c_int64_t), intent(in)        :: SHELL_PRIM_INDEX(shell_num)
     double precision   , intent(in)        :: EXPONENT(prim_num)
     double precision   , intent(in)        :: COEFFICIENT(prim_num)
   end function qmckl_context_set_ao_basis
end interface

TODO Test

Precision

The following functions set and get the expected required precision and range. precision should be an integer between 2 and 53, and range should be an integer between 2 and 11.

The setter functions functions return a new context as a 64-bit integer. The getter functions return the value, as a 32-bit integer. The update functions return QMCKL_SUCCESS or QMCKL_FAILURE.

qmckl_context_update_precision

Modifies the parameter for the numerical precision in a given context.

qmckl_exit_code qmckl_context_update_precision(const qmckl_context context, const int precision);

Source

qmckl_exit_code qmckl_context_update_precision(const qmckl_context context, const int precision) {

  if (precision <  2) return QMCKL_FAILURE;
  if (precision > 53) return QMCKL_FAILURE;

  qmckl_context_struct* ctx = (qmckl_context_struct*) context;
  if (ctx == NULL) return QMCKL_FAILURE;

  ctx->precision = precision;
  return QMCKL_SUCCESS;
}

Fortran interface

interface
   integer (c_int32_t) function qmckl_context_update_precision(context, precision) bind(C)
     use, intrinsic :: iso_c_binding
     integer (c_int64_t), intent(in), value :: context
     integer (c_int32_t), intent(in), value :: precision
   end function qmckl_context_update_precision
end interface

qmckl_context_update_range

Modifies the parameter for the numerical range in a given context.

qmckl_exit_code qmckl_context_update_range(const qmckl_context context, const int range);

Source

qmckl_exit_code qmckl_context_update_range(const qmckl_context context, const int range) {

  if (range <  2) return QMCKL_FAILURE;
  if (range > 11) return QMCKL_FAILURE;

  qmckl_context_struct* ctx = (qmckl_context_struct*) context;
  if (ctx == NULL) return QMCKL_FAILURE;

  ctx->range = range;
  return QMCKL_SUCCESS;
}

Fortran interface

interface
   integer (c_int32_t) function qmckl_context_update_range(context, range) bind(C)
     use, intrinsic :: iso_c_binding
     integer (c_int64_t), intent(in), value :: context
     integer (c_int32_t), intent(in), value :: range
   end function qmckl_context_update_range
end interface

qmckl_context_set_precision

Returns a copy of the context with a different precision parameter.

qmckl_context qmckl_context_set_precision(const qmckl_context context, const int precision);

Source

qmckl_context qmckl_context_set_precision(const qmckl_context context, const int precision) {
  qmckl_context new_context = qmckl_context_copy(context);
  if (new_context == 0) return 0;

  if (qmckl_context_update_precision(new_context, precision) == QMCKL_FAILURE) return 0;

  return new_context;
}

Fortran interface

interface
   integer (c_int64_t) function qmckl_context_set_precision(context, precision) bind(C)
     use, intrinsic :: iso_c_binding
     integer (c_int64_t), intent(in), value :: context
     integer (c_int32_t), intent(in), value :: precision
   end function qmckl_context_set_precision
end interface

qmckl_context_set_range

Returns a copy of the context with a different precision parameter.

qmckl_context qmckl_context_set_range(const qmckl_context context, const int range);

Source

qmckl_context qmckl_context_set_range(const qmckl_context context, const int range) {
  qmckl_context new_context = qmckl_context_copy(context);
  if (new_context == 0) return 0;

  if (qmckl_context_update_range(new_context, range) == QMCKL_FAILURE) return 0;

  return new_context;
}

Fortran interface

interface
   integer (c_int64_t) function qmckl_context_set_range(context, range) bind(C)
     use, intrinsic :: iso_c_binding
     integer (c_int64_t), intent(in), value :: context
     integer (c_int32_t), intent(in), value :: range
   end function qmckl_context_set_range
end interface

qmckl_context_get_precision

Returns the value of the numerical precision in the context

int32_t qmckl_context_get_precision(const qmckl_context context);

Source

int qmckl_context_get_precision(const qmckl_context context) {
  const qmckl_context_struct* ctx = (qmckl_context_struct*) context;
  return ctx->precision;
}

Fortran interface

interface
   integer (c_int32_t) function qmckl_context_get_precision(context) bind(C)
     use, intrinsic :: iso_c_binding
     integer (c_int64_t), intent(in), value :: context
   end function qmckl_context_get_precision
end interface

qmckl_context_get_range

Returns the value of the numerical range in the context

int32_t qmckl_context_get_range(const qmckl_context context);

Source

int qmckl_context_get_range(const qmckl_context context) {
  const qmckl_context_struct* ctx = (qmckl_context_struct*) context;
  return ctx->range;
}

Fortran interface

interface
   integer (c_int32_t) function qmckl_context_get_range(context) bind(C)
     use, intrinsic :: iso_c_binding
     integer (c_int64_t), intent(in), value :: context
   end function qmckl_context_get_range
end interface

qmckl_context_get_epsilon

Returns \(\epsilon = 2^{1-n}\) where n is the precision

double qmckl_context_get_epsilon(const qmckl_context context);

Source

double qmckl_context_get_epsilon(const qmckl_context context) {
  const qmckl_context_struct* ctx = (qmckl_context_struct*) context;
  return pow(2.0,(double) 1-ctx->precision);
}

Fortran interface

interface
   real (c_double) function qmckl_context_get_epsilon(context) bind(C)
     use, intrinsic :: iso_c_binding
     integer (c_int64_t), intent(in), value :: context
   end function qmckl_context_get_epsilon
end interface

2.3 Error handling

This file is written in C because it is more natural to express the error handling in C than in Fortran.

2 files are produced:

a source file : qmckl_error.c
a test file : test_qmckl_error.c

2.3.1 Error handling

The library should never make the calling programs abort, nor perform any input/output operations. This decision has to be taken by the developer of the code calling the library.

All the functions return with an exit code, defined as

typedef int32_t qmckl_exit_code;

The exit code returns the completion status of the function to the calling program. When a function call completed successfully, the QMCKL_SUCCESS exit code is returned. If one of the functions of the library fails to complete the requested task, an appropriate error code is returned to the program.

Here is the complete list of exit codes.

`QMCKL_SUCCESS`	0
`QMCKL_INVALID_ARG_1`	1
`QMCKL_INVALID_ARG_2`	2
`QMCKL_INVALID_ARG_3`	3
`QMCKL_INVALID_ARG_4`	4
`QMCKL_INVALID_ARG_5`	5
`QMCKL_INVALID_ARG_6`	6
`QMCKL_INVALID_ARG_7`	7
`QMCKL_INVALID_ARG_8`	8
`QMCKL_INVALID_ARG_9`	9
`QMCKL_INVALID_ARG_10`	10
`QMCKL_NULL_CONTEXT`	101
`QMCKL_FAILURE`	102
`QMCKL_ERRNO`	103
`QMCKL_INVALID_EXIT_CODE`	104

""" This script generates the C and Fortran constants for the error
    codes from the org-mode table.
"""

result = [ "#+BEGIN_SRC C :comments org :tangle qmckl.h" ] 
for (text, code) in table:
    text=text.replace("~","")
    result += [ f"#define  {text:30s} {code:d}" ]
result += [ "#+END_SRC" ]
    
result += [ "" ] 

result += [ "#+BEGIN_SRC f90 :comments org :tangle qmckl_f.f90" ]
for (text, code) in table:
    text=text.replace("~","")
    result += [ f"   integer, parameter :: {text:30s} = {code:d}" ]
result += [ "#+END_SRC" ]

return '\n'.join(result)

#define  QMCKL_SUCCESS                  0
#define  QMCKL_INVALID_ARG_1            1
#define  QMCKL_INVALID_ARG_2            2
#define  QMCKL_INVALID_ARG_3            3
#define  QMCKL_INVALID_ARG_4            4
#define  QMCKL_INVALID_ARG_5            5
#define  QMCKL_INVALID_ARG_6            6
#define  QMCKL_INVALID_ARG_7            7
#define  QMCKL_INVALID_ARG_8            8
#define  QMCKL_INVALID_ARG_9            9
#define  QMCKL_INVALID_ARG_10           10
#define  QMCKL_NULL_CONTEXT             101
#define  QMCKL_FAILURE                  102
#define  QMCKL_ERRNO                    103
#define  QMCKL_INVALID_EXIT_CODE        104

integer, parameter :: QMCKL_SUCCESS                  = 0
integer, parameter :: QMCKL_INVALID_ARG_1            = 1
integer, parameter :: QMCKL_INVALID_ARG_2            = 2
integer, parameter :: QMCKL_INVALID_ARG_3            = 3
integer, parameter :: QMCKL_INVALID_ARG_4            = 4
integer, parameter :: QMCKL_INVALID_ARG_5            = 5
integer, parameter :: QMCKL_INVALID_ARG_6            = 6
integer, parameter :: QMCKL_INVALID_ARG_7            = 7
integer, parameter :: QMCKL_INVALID_ARG_8            = 8
integer, parameter :: QMCKL_INVALID_ARG_9            = 9
integer, parameter :: QMCKL_INVALID_ARG_10           = 10
integer, parameter :: QMCKL_NULL_CONTEXT             = 101
integer, parameter :: QMCKL_FAILURE                  = 102
integer, parameter :: QMCKL_ERRNO                    = 103
integer, parameter :: QMCKL_INVALID_EXIT_CODE        = 104

To make a function fail, the qmckl_failwith function should be called, such that information about the failure is stored in the context. The desired exit code is given as an argument, as well as the name of the function and an error message. The return code of the function is the desired return code.

qmckl_exit_code qmckl_failwith(qmckl_context context,
                               const qmckl_exit_code exit_code,
                               const char* function,
                               const char* message) ;

qmckl_exit_code qmckl_failwith(qmckl_context context,
                               const qmckl_exit_code exit_code,
                               const char* function,
                               const char* message) {
  if (context == 0) return QMCKL_NULL_CONTEXT;
  assert (exit_code > 0);
  assert (exit_code < QMCKL_INVALID_EXIT_CODE);
  assert (function != NULL);
  assert (message  != NULL);
  assert (strlen(function) < QMCKL_MAX_FUN_LEN);
  assert (strlen(message)  < QMCKL_MAX_MSG_LEN);

  context = qmckl_context_set_error(context, exit_code, function, message);
  return exit_code;
}

For example, this function can be used as

if (x < 0) {
  return qmckl_failwith(context,
                        QMCKL_INVALID_ARG_2,
                        "qmckl_function", 
                        "Expected x >= 0");
 }

2.3.2 Multi-precision related constants

Controlling numerical precision enables optimizations. Here, the default parameters determining the target numerical precision and range are defined.

#define QMCKL_DEFAULT_PRECISION 53
#define QMCKL_DEFAULT_RANGE     11

integer, parameter :: QMCKL_DEFAULT_PRECISION = 53
integer, parameter :: QMCKL_DEFAULT_RANGE = 11

2.4 Memory management

We override the allocation functions to enable the possibility of optimized libraries to fine-tune the memory allocation.

2 files are produced:

a source file : qmckl_memory.c
a test file : test_qmckl_memory.c

2.4.1 `qmckl_malloc`

Memory allocation function, letting the library choose how the memory will be allocated, and a pointer is returned to the user. The context is passed to let the library store data related to the allocation inside the context.

void* qmckl_malloc(const qmckl_context ctx, const size_t size);

interface
   type (c_ptr) function qmckl_malloc (context, size) bind(C)
     use, intrinsic :: iso_c_binding
     integer (c_int64_t), intent(in), value :: context
     integer (c_int64_t), intent(in), value :: size
   end function qmckl_malloc
end interface

Source

void* qmckl_malloc(const qmckl_context ctx, const size_t size) {
  if (ctx == (qmckl_context) 0) {}; /* Avoid unused argument warning */
  void * result = malloc( (size_t) size );
  assert (result != NULL) ;
  return result;
}

2.4.2 `qmckl_free`

The context is passed, in case some important information has been stored related to memory allocation and needs to be updated.

qmckl_exit_code qmckl_free(qmckl_context context, void *ptr);

interface
   integer (c_int32_t) function qmckl_free (context, ptr) bind(C)
     use, intrinsic :: iso_c_binding
     integer (c_int64_t), intent(in), value :: context
     type (c_ptr), intent(in), value :: ptr
   end function qmckl_free
end interface

Source

qmckl_exit_code qmckl_free(qmckl_context context, void *ptr) {
  if (context == 0) return QMCKL_INVALID_ARG_1;
  if (ptr == NULL)  return QMCKL_INVALID_ARG_2;
  free(ptr);
  return QMCKL_SUCCESS;
}

2.5 Computation of distances

Function for the computation of distances between particles.

3 files are produced:

a source file : qmckl_distance.f90
a C test file : test_qmckl_distance.c
a Fortran test file : test_qmckl_distance_f.f90

2.5.1 Squared distance

qmckl_distance_sq

Computes the matrix of the squared distances between all pairs of points in two sets, one point within each set: \[ C_{ij} = \sum_{k=1}^3 (A_{k,i}-B_{k,j})^2 \]

Arguments

`context`	input	Global state
`transa`	input	Array `A` is `N`: Normal, `T`: Transposed
`transb`	input	Array `B` is `N`: Normal, `T`: Transposed
`m`	input	Number of points in the first set
`n`	input	Number of points in the second set
`A(lda,3)`	input	Array containing the \(m \times 3\) matrix \(A\)
`lda`	input	Leading dimension of array `A`
`B(ldb,3)`	input	Array containing the \(n \times 3\) matrix \(B\)
`ldb`	input	Leading dimension of array `B`
`C(ldc,n)`	output	Array containing the \(m \times n\) matrix \(C\)
`ldc`	input	Leading dimension of array `C`

Requirements
- context is not 0
- m > 0
- n > 0
- lda >= 3 if transa is N
- lda >= m if transa is T
- ldb >= 3 if transb is N
- ldb >= n if transb is T
- ldc >= m
- A is allocated with at least \(3 \times m \times 8\) bytes
- B is allocated with at least \(3 \times n \times 8\) bytes
- C is allocated with at least \(m \times n \times 8\) bytes

Performance

This function might be more efficient when A and B are transposed.

qmckl_exit_code qmckl_distance_sq(const qmckl_context context,
                                  const char transa, const char transb,
                                  const int64_t m, const int64_t n,
                                  const double *A, const int64_t lda,
                                  const double *B, const int64_t ldb,
                                  const double *C, const int64_t ldc);

Source

integer function qmckl_distance_sq_f(context, transa, transb, m, n, A, LDA, B, LDB, C, LDC) result(info)
  implicit none
  integer*8  , intent(in)  :: context
  character  , intent(in)  :: transa, transb
  integer*8  , intent(in)  :: m, n
  integer*8  , intent(in)  :: lda
  real*8     , intent(in)  :: A(lda,*)
  integer*8  , intent(in)  :: ldb
  real*8     , intent(in)  :: B(ldb,*)
  integer*8  , intent(in)  :: ldc
  real*8     , intent(out) :: C(ldc,*)

  integer*8 :: i,j
  real*8    :: x, y, z
  integer   :: transab

  info = 0

  if (context == 0_8) then
     info = -1
     return
  endif

  if (m <= 0_8) then
     info = -2
     return
  endif

  if (n <= 0_8) then
     info = -3
     return
  endif

  if (transa == 'N' .or. transa == 'n') then
     transab = 0
  else if (transa == 'T' .or. transa == 't') then
     transab = 1
  else
     transab = -100
  endif

  if (transb == 'N' .or. transb == 'n') then
     continue
  else if (transa == 'T' .or. transa == 't') then
     transab = transab + 2
  else
     transab = -100
  endif

  if (transab < 0) then
     info = -4
     return 
  endif

  if (iand(transab,1) == 0 .and. LDA < 3) then
     info = -5
     return
  endif

  if (iand(transab,1) == 1 .and. LDA < m) then
     info = -6
     return
  endif

  if (iand(transab,2) == 0 .and. LDA < 3) then
     info = -6
     return
  endif

  if (iand(transab,2) == 2 .and. LDA < m) then
     info = -7
     return
  endif


  select case (transab)

  case(0)

     do j=1,n
        do i=1,m
           x = A(1,i) - B(1,j)
           y = A(2,i) - B(2,j)
           z = A(3,i) - B(3,j)
           C(i,j) = x*x + y*y + z*z
        end do
     end do

  case(1)

     do j=1,n
        do i=1,m
           x = A(i,1) - B(1,j)
           y = A(i,2) - B(2,j)
           z = A(i,3) - B(3,j)
           C(i,j) = x*x + y*y + z*z
        end do
     end do

  case(2)

     do j=1,n
        do i=1,m
           x = A(1,i) - B(j,1)
           y = A(2,i) - B(j,2)
           z = A(3,i) - B(j,3)
           C(i,j) = x*x + y*y + z*z
        end do
     end do

  case(3)

     do j=1,n
        do i=1,m
           x = A(i,1) - B(j,1)
           y = A(i,2) - B(j,2)
           z = A(i,3) - B(j,3)
           C(i,j) = x*x + y*y + z*z
        end do
     end do

  end select

end function qmckl_distance_sq_f

2.6 Atomic Orbitals

This files contains all the routines for the computation of the values, gradients and Laplacian of the atomic basis functions.

3 files are produced:

a source file : qmckl_ao.f90
a C test file : test_qmckl_ao.c
a Fortran test file : test_qmckl_ao_f.f90

2.6.1 Polynomials

\[ P_l(\mathbf{r},\mathbf{R}_i) = (x-X_i)^a (y-Y_i)^b (z-Z_i)^c \]

\begin{eqnarray*} \frac{\partial }{\partial x} P_l\left(\mathbf{r},\mathbf{R}_i \right) & = & a (x-X_i)^{a-1} (y-Y_i)^b (z-Z_i)^c \\ \frac{\partial }{\partial y} P_l\left(\mathbf{r},\mathbf{R}_i \right) & = & b (x-X_i)^a (y-Y_i)^{b-1} (z-Z_i)^c \\ \frac{\partial }{\partial z} P_l\left(\mathbf{r},\mathbf{R}_i \right) & = & c (x-X_i)^a (y-Y_i)^b (z-Z_i)^{c-1} \\ \end{eqnarray*} \begin{eqnarray*} \left( \frac{\partial }{\partial x^2} + \frac{\partial }{\partial y^2} + \frac{\partial }{\partial z^2} \right) P_l \left(\mathbf{r},\mathbf{R}_i \right) & = & a(a-1) (x-X_i)^{a-2} (y-Y_i)^b (z-Z_i)^c + \\ && b(b-1) (x-X_i)^a (y-Y_i)^{b-1} (z-Z_i)^c + \\ && c(c-1) (x-X_i)^a (y-Y_i)^b (z-Z_i)^{c-1} \end{eqnarray*}

qmckl_ao_power

Computes all the powers of the n input data up to the given maximum value given in input for each of the \(n\) points:

\[ P_{ij} = X_j^i \]

Arguments

`context`	input	Global state
`n`	input	Number of values
`X(n)`	input	Array containing the input values
`LMAX(n)`	input	Array containing the maximum power for each value
`P(LDP,n)`	output	Array containing all the powers of `X`
`LDP`	input	Leading dimension of array `P`

Requirements
- context is not 0
- n > 0
- X is allocated with at least \(n \times 8\) bytes
- LMAX is allocated with at least \(n \times 4\) bytes
- P is allocated with at least \(n \times \max_i \text{LMAX}_i \times 8\) bytes
- LDP >= \(\max_i\) LMAX[i]

Header

qmckl_exit_code qmckl_ao_power(const qmckl_context context,
                const int64_t n, 
                const double *X, const int32_t *LMAX,
                const double *P, const int64_t LDP);

Source

integer function qmckl_ao_power_f(context, n, X, LMAX, P, ldp) result(info)
  implicit none
  integer*8 , intent(in)  :: context
  integer*8 , intent(in)  :: n
  real*8    , intent(in)  :: X(n)
  integer   , intent(in)  :: LMAX(n)
  real*8    , intent(out) :: P(ldp,n)
  integer*8 , intent(in)  :: ldp

  integer*8  :: i,j

  info = 0

  if (context == 0_8) then
     info = -1
     return
  endif

  if (LDP < MAXVAL(LMAX)) then
     info = -2
     return
  endif

  do j=1,n
    P(1,j) = X(j)
    do i=2,LMAX(j)
       P(i,j) = P(i-1,j) * X(j) 
    end do
  end do

end function qmckl_ao_power_f

qmckl_ao_polynomial_vgl

Computes the values, gradients and Laplacians at a given point of all polynomials with an angular momentum up to lmax.

Arguments

`context`	input	Global state
`X(3)`	input	Array containing the coordinates of the points
`R(3)`	input	Array containing the x,y,z coordinates of the center
`lmax`	input	Maximum angular momentum
`n`	output	Number of computed polynomials
`L(ldl,n)`	output	Contains a,b,c for all `n` results
`ldl`	input	Leading dimension of `L`
`VGL(ldv,n)`	output	Value, gradients and Laplacian of the polynomials
`ldv`	input	Leading dimension of array `VGL`

Requirements
- context is not 0
- n > 0
- lmax >= 0
- ldl >= 3
- ldv >= 5
- X is allocated with at least \(3 \times 8\) bytes
- R is allocated with at least \(3 \times 8\) bytes
- n >= (lmax+1)(lmax+2)(lmax+3)/6
- L is allocated with at least \(3 \times n \times 4\) bytes
- VGL is allocated with at least \(5 \times n \times 8\) bytes
- On output, n should be equal to (lmax+1)(lmax+2)(lmax+3)/6
- On output, the powers are given in the following order (l=a+b+c):
  - Increase values of l
  - Within a given value of l, alphabetical order of the string made by a*"x" + b*"y" + c*"z" (in Python notation). For example, with a=0, b=2 and c=1 the string is "yyz"
Error codes

-1 Null context

-2 Inconsistent ldl

-3 Inconsistent ldv

-4 Inconsistent lmax

Header

qmckl_exit_code qmckl_ao_polynomial_vgl(const qmckl_context context,
                const double *X, const double *R,
                const int32_t lmax, const int64_t *n,
                const int32_t *L,   const int64_t ldl,
                const double *VGL,  const int64_t ldv);

Source

integer function qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL, ldv) result(info)
  implicit none
  integer*8 , intent(in)  :: context
  real*8    , intent(in)  :: X(3), R(3)
  integer   , intent(in)  :: lmax
  integer*8 , intent(out) :: n
  integer   , intent(out) :: L(ldl,(lmax+1)*(lmax+2)*(lmax+3)/6)
  integer*8 , intent(in)  :: ldl
  real*8    , intent(out) :: VGL(ldv,(lmax+1)*(lmax+2)*(lmax+3)/6)
  integer*8 , intent(in)  :: ldv

  integer*8         :: i,j
  integer           :: a,b,c,d
  real*8            :: Y(3)
  integer           :: lmax_array(3)
  real*8            :: pows(-2:lmax,3)
  integer, external :: qmckl_ao_power_f
  double precision  :: xy, yz, xz
  double precision  :: da, db, dc, dd

  info = 0

  if (context == 0_8) then
     info = -1
     return
  endif

  if (ldl < 3) then
     info = -2
     return
  endif

  if (ldv < 5) then
     info = -3
     return
  endif

  if (lmax <= 0) then
     info = -4
     return
  endif


  do i=1,3
     Y(i) = X(i) - R(i)
  end do

  lmax_array(1:3) = lmax
  if (lmax == 0) then
    VGL(1,1) = 1.d0
    vgL(2:5,1) = 0.d0
    l(1:3,1) = 0
    n=1
  else if (lmax > 0) then
    pows(-2:0,1:3) = 1.d0
    do i=1,lmax
        pows(i,1) = pows(i-1,1) * Y(1) 
        pows(i,2) = pows(i-1,2) * Y(2) 
        pows(i,3) = pows(i-1,3) * Y(3) 
    end do

    VGL(1:5,1:4) = 0.d0
    l(1:3,1:4) = 0

    VGL(1,1) = 1.d0
    vgl(1:5,2:4) = 0.d0

    l(1,2) = 1
    vgl(1,2) = pows(1,1)
    vgL(2,2) = 1.d0

    l(2,3) = 1
    vgl(1,3) = pows(1,2)
    vgL(3,3) = 1.d0

    l(3,4) = 1
    vgl(1,4) = pows(1,3)
    vgL(4,4) = 1.d0

    n=4
  endif

  ! l>=2
  dd = 2.d0
  do d=2,lmax
     da = dd
     do a=d,0,-1
        db = dd-da
        do b=d-a,0,-1
           c  = d  - a  - b
           dc = dd - da - db
           n = n+1

           l(1,n) = a
           l(2,n) = b
           l(3,n) = c

           xy = pows(a,1) * pows(b,2)
           yz = pows(b,2) * pows(c,3)
           xz = pows(a,1) * pows(c,3)

           vgl(1,n) = xy * pows(c,3)

           xy = dc * xy
           xz = db * xz
           yz = da * yz

           vgl(2,n) = pows(a-1,1) * yz
           vgl(3,n) = pows(b-1,2) * xz
           vgl(4,n) = pows(c-1,3) * xy

           vgl(5,n) = &
                (da-1.d0) * pows(a-2,1) * yz + &
                (db-1.d0) * pows(b-2,2) * xz + &
                (dc-1.d0) * pows(c-2,3) * xy

           db = db - 1.d0
        end do
        da = da - 1.d0
     end do
     dd = dd + 1.d0
  end do

  info = 0

end function qmckl_ao_polynomial_vgl_f

2.6.2 Gaussian basis functions

qmckl_ao_gaussian_vgl

Computes the values, gradients and Laplacians at a given point of n Gaussian functions centered at the same point:

\[ v_i = \exp(-a_i |X-R|^2) \] \[ \nabla_x v_i = -2 a_i (X_x - R_x) v_i \] \[ \nabla_y v_i = -2 a_i (X_y - R_y) v_i \] \[ \nabla_z v_i = -2 a_i (X_z - R_z) v_i \] \[ \Delta v_i = a_i (4 |X-R|^2 a_i - 6) v_i \]

Arguments

`context`	input	Global state
`X(3)`	input	Array containing the coordinates of the points
`R(3)`	input	Array containing the x,y,z coordinates of the center
`n`	input	Number of computed gaussians
`A(n)`	input	Exponents of the Gaussians
`VGL(ldv,5)`	output	Value, gradients and Laplacian of the Gaussians
`ldv`	input	Leading dimension of array `VGL`

Requirements
- context is not 0
- n > 0
- ldv >= 5
- A(i) > 0 for all i
- X is allocated with at least \(3 \times 8\) bytes
- R is allocated with at least \(3 \times 8\) bytes
- A is allocated with at least \(n \times 8\) bytes
- VGL is allocated with at least \(n \times 5 \times 8\) bytes

Header

qmckl_exit_code qmckl_ao_gaussian_vgl(const qmckl_context context,
                const double *X, const double *R,
                const int64_t *n, const int64_t *A,
                const double *VGL,  const int64_t ldv);

Source

integer function qmckl_ao_gaussian_vgl_f(context, X, R, n, A, VGL, ldv) result(info)
  implicit none
  integer*8 , intent(in)  :: context
  real*8    , intent(in)  :: X(3), R(3)
  integer*8 , intent(in)  :: n
  real*8    , intent(in)  :: A(n)
  real*8    , intent(out) :: VGL(ldv,5)
  integer*8 , intent(in)  :: ldv

  integer*8         :: i,j
  real*8            :: Y(3), r2, t, u, v

  info = 0

  if (context == 0_8) then
     info = -1
     return
  endif

  if (n <= 0) then
     info = -2
     return
  endif

  if (ldv < n) then
     info = -3
     return
  endif


  do i=1,3
     Y(i) = X(i) - R(i)
  end do
  r2 = Y(1)*Y(1) + Y(2)*Y(2) + Y(3)*Y(3)

  do i=1,n
     VGL(i,1) = dexp(-A(i) * r2)
  end do

  do i=1,n
     VGL(i,5) = A(i) * VGL(i,1)
  end do

  t = -2.d0 * ( X(1) - R(1) )
  u = -2.d0 * ( X(2) - R(2) )
  v = -2.d0 * ( X(3) - R(3) )

  do i=1,n
     VGL(i,2) = t * VGL(i,5)
     VGL(i,3) = u * VGL(i,5)
     VGL(i,4) = v * VGL(i,5)
  end do

  t = 4.d0 * r2
  do i=1,n
     VGL(i,5) = (t * A(i) - 6.d0) *  VGL(i,5)
  end do

end function qmckl_ao_gaussian_vgl_f

2.6.3 TODO Slater basis functions

3 Acknowledgments

TREX: Targeting Real Chemical Accuracy at the Exascale project has received funding from the European Union’s Horizon 2020 - Research and Innovation program - under grant agreement no. 952165. The content of this document does not represent the opinion of the European Union, and the European Union is not responsible for any use that might be made of such content.

-1	Null context
-2	Inconsistent `ldl`
-3	Inconsistent `ldv`
-4	Inconsistent `lmax`

QMCkl source code documentation

Table of Contents

1 Introduction

1.1 Literate programming

1.2 Source code editing

1.3 Choice of the programming language

1.4 Design of the library

1.4.1 Naming conventions

1.4.2 Application programming interface

1.4.3 Global state

1.4.4 Low-level functions

1.4.5 High-level functions

1.4.6 Numerical precision

1.5 Algorithms

1.6 Rules for the API

2 Documentation

2.1 qmckl.h header file

2.2 Context

2.2.1 Context

2.3 Error handling

2.3.1 Error handling

2.3.2 Multi-precision related constants

2.4 Memory management

2.4.1 qmckl_malloc

2.4.2 qmckl_free

2.5 Computation of distances

2.5.1 Squared distance

2.6 Atomic Orbitals

2.6.1 Polynomials

2.6.2 Gaussian basis functions

2.6.3 TODO Slater basis functions

3 Acknowledgments

2.1 `qmckl.h` header file

2.4.1 `qmckl_malloc`

2.4.2 `qmckl_free`