QMCkl source code documentation

1. Introduction
2. Documentation
3. Acknowledgments

1 Introduction

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

Literate programming is particularly adapted in this context. Source files are written in org-mode format, to provide useful comments and LaTex formulas close to the code. There exists multiple possibilities to convert org-mode files into different formats such as HTML or pdf. For a tutorial on literate programming with org-mode, follow this link.

The code is extracted from the org files using Emacs as a command-line tool in the Makefile, and then the produced files are compiled.

1.1 Language used

Fortran is one of the most common languages used by the community, and is simple enough to make the algorithms readable. Hence we propose in this pedagogical implementation of QMCkl to use Fortran to express the algorithms. For specific internal functions where the C language is more natural, C is used.

As Fortran modules generate compiler-dependent files, the use of modules is restricted to the internal use of the library, otherwise the compliance with C is violated.

The external dependencies should be kept as small as possible, so external libraries should be used only if their used is strongly justified.

1.2 Source code editing

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

For users with a preference for Jupyter notebooks, the following script can convert jupyter notebooks to org-mode files:

#!/bin/bash
# $ nb_to_org.sh notebook.ipynb
# produces the org-mode file notebook.org

set -e

nb=$(basename $1 .ipynb)
jupyter nbconvert --to markdown ${nb}.ipynb --output ${nb}.md
pandoc ${nb}.md -o ${nb}.org
rm ${nb}.md

And pandoc can convert multiple markdown formats into org-mode.

1.3 Writing in Fortran

The Fortran source files should provide a C interface using iso_c_binding. The name of the Fortran source files should end with _f.f90 to be properly handled by the Makefile. The names of the functions defined in fortran should be the same as those exposed in the API suffixed by _f. Fortran interface files should also be written in the qmckl_f.f90 file.

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

1.4 Coding style

To improve readability, we maintain a consistent coding style in the library.

For C source files, we will use (decide on a coding style)
For Fortran source files, we will use (decide on a coding style)

Coding style can be automatically checked with clang-format.

1.5 Design of the library

The proposed API should allow the library to:

deal with memory transfers between CPU and accelerators
use different levels of floating-point precision

We chose a multi-layered design with low-level and high-level functions (see below).

1.5.1 Naming conventions

Use qmckl_ as a prefix for all exported functions and variables. All exported header files should have a filename with the prefix qmckl_.

If the name of the org-mode file is xxx.org, the name of the produced C files should be xxx.c and xxx.h and the name of the produced Fortran files should be xxx.f90

Arrays are in uppercase and scalars are in lowercase.

In the names of the variables and functions, only the singular form is allowed.

1.5.2 Application programming interface

The application programming interface (API) is designed to be compatible with the C programming language (not C++), to ensure that the library will be easily usable in any language. This implies that only the following data types are allowed in the API:

32-bit and 64-bit floats and arrays (real and double)
32-bit and 64-bit integers and arrays (int32_t and int64_t)
Pointers should be represented as 64-bit integers (even on 32-bit architectures)
ASCII strings are represented as a pointers to a character arrays and terminated by a zero character (C convention).

Complex numbers can be represented by an array of 2 floats.

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

1.5.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.5.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.5.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.5.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.6 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.7 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

This file produces the qmckl.h header file, which is to be included when qmckl functions are used.

We also create here the qmckl_f.f90 which is the Fortran interface file.

2.1.1 Constants

Success/failure

These are the codes returned by the functions to indicate success or failure. All such functions should have as a return type qmckl_exit_code.

#define QMCKL_SUCCESS 0
#define QMCKL_FAILURE 1

typedef int32_t qmckl_exit_code;
typedef int64_t qmckl_context ;

integer, parameter :: QMCKL_SUCCESS = 0
integer, parameter :: QMCKL_FAILURE = 0

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.2 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.2.1 `qmckl_malloc`

Memory allocation function, letting the library choose how the memory will be allocated, and a pointer is returned to the user.

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) {
    /* Avoids unused parameter error */
    return malloc( (size_t) size );
  }
  return malloc( (size_t) size );
}

2.2.2 `qmckl_free`

void qmckl_free(void *ptr);

interface
   subroutine qmckl_free (ptr) bind(C)
     use, intrinsic :: iso_c_binding
     type (c_ptr), intent(in), value :: ptr
   end subroutine qmckl_free
end interface

Source

void qmckl_free(void *ptr) {
  free(ptr);
}

2.3 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.3.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.

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;

} qmckl_context_struct;

#define VALID_TAG   0xBEEFFACE
#define INVALID_TAG 0xDEADBEEF

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;

  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;

  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;
  qmckl_free(ctx);
  return QMCKL_SUCCESS;
}

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

2.3.2 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) {
    free(basis);
    return QMCKL_FAILURE;
  }

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

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

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

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

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


  /* Assign data */

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

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

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

  ctx->ao_basis = basis;
  return QMCKL_SUCCESS;
}

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(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

2.3.3 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(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(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.4 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.4.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.5 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.5.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.5.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.5.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 Language used

1.2 Source code editing

1.3 Writing in Fortran

1.4 Coding style

1.5 Design of the library

1.5.1 Naming conventions

1.5.2 Application programming interface

1.5.3 Global state

1.5.4 Low-level functions

1.5.5 High-level functions

1.5.6 Numerical precision

1.6 Algorithms

1.7 Rules for the API

2 Documentation

2.1 qmckl.h header file

2.1.1 Constants

2.2 Memory management

2.2.1 qmckl_malloc

2.2.2 qmckl_free

2.3 Context

2.3.1 Context

2.3.2 Basis set

2.3.3 Precision

2.4 Computation of distances

2.4.1 Squared distance

2.5 Atomic Orbitals

2.5.1 Polynomials

2.5.2 Gaussian basis functions

2.5.3 TODO Slater basis functions

3 Acknowledgments

2.1 `qmckl.h` header file

2.2.1 `qmckl_malloc`

2.2.2 `qmckl_free`