UP | HOME

Atomic Orbitals

Table of Contents

1 Introduction

The atomic basis set is defined as a list of shells. Each shell \(s\) is centered on a nucleus \(A\), possesses a given angular momentum \(l\) and a radial function \(R_s\). The radial function is a linear combination of primitive functions that can be of type Slater (\(p=1\)) or Gaussian (\(p=2\)):

\[ R_s(\mathbf{r}) = \mathcal{N}_s |\mathbf{r}-\mathbf{R}_A|^{n_s} \sum_{k=1}^{N_{\text{prim}}} a_{ks}\, f_{ks} \exp \left( - \gamma_{ks} | \mathbf{r}-\mathbf{R}_A | ^p \right). \]

In the case of Gaussian functions, \(n_s\) is always zero. The normalization factor \(\mathcal{N}_s\) ensures that all the functions of the shell are normalized (integrate) to unity. Usually, basis sets are given a combination of normalized primitives, so the normalization coefficients of the primitives, \(f_{ks}\), need also to be provided.

Atomic orbitals (AOs) are defined as

\[ \chi_i (\mathbf{r}) = \mathcal{M}_i\, P_{\eta(i)}(\mathbf{r})\, R_{\theta(i)} (\mathbf{r}) \]

where \(\theta(i)\) returns the shell on which the AO is expanded, and \(\eta(i)\) denotes which angular function is chosen and \(P\) are the generating functions of the given angular momentum \(\eta(i)\). Here, the parameter \(\mathcal{M}_i\) is an extra parameter which allows the normalization of the different functions of the same shell to be different, as in GAMESS for example.

In this section we describe first how the basis set is stored in the context, and then we present the kernels used to compute the values, gradients and Laplacian of the atomic basis functions.

2 Context

2.1 Constant data

The following arrays are stored in the context, and need to be set when initializing the library:

Variable Type Description
type char Gaussian ('G') or Slater ('S')
shell_num int64_t Number of shells
prim_num int64_t Total number of primitives
nucleus_index int64_t[nucl_num] Index of the first shell of each nucleus
nucleus_shell_num int64_t[nucl_num] Number of shells per nucleus
shell_ang_mom int32_t[shell_num] Angular momentum of each shell
shell_prim_num int64_t[shell_num] Number of primitives in each shell
shell_prim_index int64_t[shell_num] Address of the first primitive of each shell in the EXPONENT array
shell_factor double[shell_num] Normalization factor for each shell
exponent double[prim_num] Array of exponents
coefficient double[prim_num] Array of coefficients
prim_factor double[prim_num] Normalization factors of the primtives
ao_num int64_t Number of AOs
ao_cartesian bool If true, use polynomials. Otherwise, use spherical harmonics
ao_factor double[ao_num] Normalization factor of the AO

The following data is computed when the basis is finalized:

Variable Type Description
nucleus_prim_index int64_t[nucl_num+1] Index of the first primitive of each nucleus
nucleus_max_ang_mom int32_t[nucl_num] Maximum angular momentum of each nucleus
coefficient_normalized~ double[prim_num] Normalized array of coefficients
ao_ang_mom int32_t[ao_num] Angular momentum of the shell to which the AO belongs
ao_nucl int64_t[ao_num] Nucleus on which the AO is centered
nucleus_range double[nucl_num] Distance beyond which all AOs are zero

For H2 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
ao_num = 30

nucleus_index = [0 , 6]

shell_ang_mom = [0, 0, 0, 1, 1, 2, 0, 0, 0, 1, 1, 2]
shell_factor  = [ 1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.]
shell_prim_num = [5, 1, 1, 1, 1, 1, 5, 1, 1, 1, 1, 1]
shell_prim_index = [0 , 5 , 6 , 7 , 8 , 9 , 10, 15, 16, 17, 18, 19]

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]

prim_factor = [ 1.0006253235944540e+01, 2.4169531573445120e+00, 7.9610924849766440e-01
  3.0734305383061117e-01, 1.2929684417481876e-01, 3.0734305383061117e-01,
  1.2929684417481876e-01, 2.1842769845268308e+00, 4.3649547399719840e-01,
  1.8135965626177861e+00, 1.0006253235944540e+01, 2.4169531573445120e+00,
  7.9610924849766440e-01, 3.0734305383061117e-01, 1.2929684417481876e-01,
  3.0734305383061117e-01, 1.2929684417481876e-01, 2.1842769845268308e+00,
  4.3649547399719840e-01, 1.8135965626177861e+00 ]

A scalar variable $V$ present in this table can be set or get by calling the functions: 

qmckl_exit_code qmckl_set_ao_basis_$V$ ( qmckl_context context,
                                         const $type_of_V$ $V$);

qmckl_exit_code qmckl_get_ao_basis_$V$ ( const qmckl_context context,
                                         $type_of_V$ const $V$);

interface
   integer(c_int32_t) function qmckl_set_ao_basis_$V$ (context, $V$) &
        bind(C)
     use, intrinsic :: iso_c_binding
     import
     implicit none
     integer (c_int64_t) , intent(in)  , value :: context
     $f_type_of_V$       , intent(in)  , value :: $V$
   end function qmckl_set_ao_basis_$V$
end interface

interface
   integer(c_int32_t) function qmckl_get_ao_basis_$V$ (context, $V$) &
        bind(C)
     use, intrinsic :: iso_c_binding
     import
     implicit none
     integer (c_int64_t) , intent(in)  , value :: context
     $f_type_of_V$       , intent(out)         :: $V$
   end function qmckl_get_ao_basis_$V$
end interface

For array variables, use the rule: 

qmckl_exit_code qmckl_set_ao_basis_$V$ ( qmckl_context context,
                                         const $type_of_V$ $V$,
                                         const int64_t size_max);

qmckl_exit_code qmckl_get_ao_basis_$V$ ( const qmckl_context context,
                                         $type_of_V$ const $V$,
                                         const int64_t size_max);

interface
   integer(c_int32_t) function qmckl_set_ao_basis_$V$ (context, &
        $V$, size_max) bind(C)
     use, intrinsic :: iso_c_binding
     import
     implicit none
     integer (c_int64_t) , intent(in)  , value :: context
     $f_type_of_V$       , intent(in)  , value :: $V$
     integer (c_int64_t) , intent(in)  , value :: size_max
   end function qmckl_set_ao_basis_$V$
end interface

interface
   integer(c_int32_t) function qmckl_get_ao_basis_$V$ (context, &
        $V$, size_max) bind(C)
     use, intrinsic :: iso_c_binding
     import
     implicit none
     integer (c_int64_t) , intent(in)  , value :: context
     $f_type_of_V$       , intent(out)         :: $V$
     integer (c_int64_t) , intent(in)  , value :: size_max
   end function qmckl_get_ao_basis_$V$
end interface

2.1.1 Initialization functions

size_max is the dimension of the input array, which should be equal of larger than the value given in the table of section 2.

2.1.1.1 C interface

To set the basis set, all the following functions need to be called. 

qmckl_exit_code
qmckl_set_ao_basis_type (qmckl_context context,
                         const char basis_type);

qmckl_exit_code
qmckl_set_ao_basis_shell_num (qmckl_context context,
                              const int64_t shell_num);

qmckl_exit_code
qmckl_set_ao_basis_prim_num (qmckl_context context,
                             const int64_t prim_num);

qmckl_exit_code
qmckl_set_ao_basis_nucleus_shell_num (qmckl_context context,
                                      const int64_t* nucleus_shell_num,
                                      const int64_t size_max);

qmckl_exit_code
qmckl_set_ao_basis_nucleus_index (qmckl_context context,
                                  const int64_t* nucleus_index,
                                  const int64_t size_max);

qmckl_exit_code
qmckl_set_ao_basis_shell_ang_mom (qmckl_context context,
                                  const int32_t* shell_ang_mom,
                                  const int64_t size_max);

qmckl_exit_code
qmckl_set_ao_basis_shell_prim_num (qmckl_context context,
                                   const int64_t* shell_prim_num,
                                   const int64_t size_max);

qmckl_exit_code
qmckl_set_ao_basis_shell_prim_index (qmckl_context context,
                                     const int64_t* shell_prim_index,
                                     const int64_t size_max);

qmckl_exit_code
qmckl_set_ao_basis_shell_factor (qmckl_context context,
                                 const double* shell_factor,
                                 const int64_t size_max);

qmckl_exit_code
qmckl_set_ao_basis_exponent (qmckl_context context,
                             const double* exponent,
                             const int64_t size_max);

qmckl_exit_code
qmckl_set_ao_basis_coefficient (qmckl_context context,
                                const double* coefficient,
                                const int64_t size_max);

qmckl_exit_code
qmckl_set_ao_basis_prim_factor (qmckl_context context,
                                const double* prim_factor,
                                const int64_t size_max);

qmckl_exit_code
qmckl_set_ao_basis_ao_num (qmckl_context context,
                           const int64_t ao_num);

qmckl_exit_code
qmckl_set_ao_basis_ao_factor (qmckl_context context,
                              const double* ao_factor,
                              const int64_t size_max);

qmckl_exit_code
qmckl_set_ao_basis_cartesian (qmckl_context context,
                              const bool cartesian);
2.1.1.2 Fortran interface
interface
  integer(c_int32_t) function qmckl_set_ao_basis_type (context, &
       basis_type) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    character(c_char)   , intent(in)  , value :: basis_type
  end function qmckl_set_ao_basis_type
end interface

interface
  integer(c_int32_t) function qmckl_set_ao_basis_shell_num(context, &
       num) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    integer (c_int64_t) , intent(in)  , value :: num
  end function qmckl_set_ao_basis_shell_num
end interface

interface
  integer(c_int32_t) function qmckl_set_ao_basis_prim_num(context, &
       num) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    integer (c_int64_t) , intent(in)  , value :: num
  end function qmckl_set_ao_basis_prim_num
end interface

interface
  integer(c_int32_t) function qmckl_set_ao_basis_nucleus_index(context, &
       idx, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    integer (c_int64_t) , intent(in)          :: idx(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_set_ao_basis_nucleus_index
end interface

interface
  integer(c_int32_t) function qmckl_set_ao_basis_nucleus_shell_num(context, &
       shell_num, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    integer (c_int64_t) , intent(in)          :: shell_num(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_set_ao_basis_nucleus_shell_num
end interface

interface
  integer(c_int32_t) function qmckl_set_ao_basis_shell_ang_mom(context, &
       shell_ang_mom, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    integer (c_int32_t) , intent(in)          :: shell_ang_mom(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_set_ao_basis_shell_ang_mom
end interface

interface
  integer(c_int32_t) function qmckl_set_ao_basis_shell_prim_num(context, &
       shell_prim_num, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    integer (c_int64_t) , intent(in)          :: shell_prim_num(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_set_ao_basis_shell_prim_num
end interface

interface
  integer(c_int32_t) function qmckl_set_ao_basis_shell_prim_index(context, &
       shell_prim_index, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    integer (c_int64_t) , intent(in)          :: shell_prim_index(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_set_ao_basis_shell_prim_index
end interface

interface
  integer(c_int32_t) function qmckl_set_ao_basis_shell_factor(context, &
       shell_factor, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    real    (c_double)  , intent(in)          :: shell_factor(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_set_ao_basis_shell_factor
end interface

interface
  integer(c_int32_t) function qmckl_set_ao_basis_exponent(context, &
       exponent, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    real    (c_double)  , intent(in)          :: exponent(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_set_ao_basis_exponent
end interface

interface
  integer(c_int32_t) function qmckl_set_ao_basis_coefficient(context, &
       coefficient, size_max)  bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    real    (c_double)  , intent(in)          :: coefficient(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_set_ao_basis_coefficient
end interface

interface
  integer(c_int32_t) function qmckl_set_ao_basis_prim_factor(context, &
       prim_factor, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    real    (c_double)  , intent(in)          :: prim_factor(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_set_ao_basis_prim_factor
end interface

interface
  integer(c_int32_t) function qmckl_set_ao_basis_ao_num(context, &
       num) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    integer (c_int64_t) , intent(in)  , value :: num
  end function qmckl_set_ao_basis_ao_num
end interface

interface
  integer(c_int32_t) function qmckl_set_ao_basis_cartesian(context, &
       cartesian) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    logical (c_bool)    , intent(in)  , value :: cartesian
  end function qmckl_set_ao_basis_cartesian
end interface

interface
  integer(c_int32_t) function qmckl_set_ao_basis_ao_factor(context, &
       ao_factor, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    real    (c_double)  , intent(in)          :: ao_factor(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_set_ao_basis_ao_factor
end interface

2.1.2 Access functions

size_max is the dimension of the input array, which should be equal of larger than the value given in the table of section 2.

2.1.2.1 C interface
qmckl_exit_code
qmckl_get_ao_basis_type (const qmckl_context context,
                         char* const basis_type);

qmckl_exit_code
qmckl_get_ao_basis_shell_num (const qmckl_context context,
                              int64_t* const shell_num);

qmckl_exit_code
qmckl_get_ao_basis_prim_num (const qmckl_context context,
                             int64_t* const prim_num);

qmckl_exit_code
qmckl_get_ao_basis_nucleus_shell_num (const qmckl_context context,
                                      int64_t* const nucleus_shell_num,
                                      const int64_t size_max);

qmckl_exit_code
qmckl_get_ao_basis_nucleus_index (const qmckl_context context,
                                  int64_t* const nucleus_index,
                                  const int64_t size_max);

qmckl_exit_code
qmckl_get_ao_basis_shell_ang_mom (const qmckl_context context,
                                  int32_t* const shell_ang_mom,
                                  const int64_t size_max);

qmckl_exit_code
qmckl_get_ao_basis_shell_prim_num (const qmckl_context context,
                                   int64_t* const shell_prim_num,
                                   const int64_t size_max);

qmckl_exit_code
qmckl_get_ao_basis_shell_prim_index (const qmckl_context context,
                                     int64_t* const shell_prim_index,
                                     const int64_t size_max);

qmckl_exit_code
qmckl_get_ao_basis_shell_factor (const qmckl_context context,
                                 double*  const shell_factor,
                                 const int64_t size_max);

qmckl_exit_code
qmckl_get_ao_basis_exponent (const qmckl_context context,
                             double*  const exponent,
                             const int64_t size_max);

qmckl_exit_code
qmckl_get_ao_basis_coefficient (const qmckl_context context,
                                double*  const coefficient,
                                const int64_t size_max);

qmckl_exit_code
qmckl_get_ao_basis_prim_factor (const qmckl_context context,
                                double*  const prim_factor,
                                const int64_t size_max);

qmckl_exit_code
qmckl_get_ao_basis_ao_num (const qmckl_context context,
                           int64_t* const ao_num);

qmckl_exit_code
qmckl_get_ao_basis_ao_factor (const qmckl_context context,
                              double* const ao_factor,
                              const int64_t size_max);

When all the data for the AOs have been provided, the following function returns true. 

bool qmckl_ao_basis_provided (const qmckl_context context);
2.1.2.2 Fortran interface
interface
  integer(c_int32_t) function qmckl_get_ao_basis_type (context, &
       basis_type) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    character(c_char)   , intent(out)         :: basis_type
  end function qmckl_get_ao_basis_type
end interface

interface
  integer(c_int32_t) function qmckl_get_ao_basis_shell_num(context, &
       num) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    integer (c_int64_t) , intent(out)         :: num
  end function qmckl_get_ao_basis_shell_num
end interface

interface
  integer(c_int32_t) function qmckl_get_ao_basis_prim_num(context, &
       num) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    integer (c_int64_t) , intent(out)         :: num
  end function qmckl_get_ao_basis_prim_num
end interface

interface
  integer(c_int32_t) function qmckl_get_ao_basis_nucleus_shell_num(context, &
       shell_num, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    integer (c_int64_t) , intent(out)         :: shell_num(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_get_ao_basis_nucleus_shell_num
end interface

interface
  integer(c_int32_t) function qmckl_get_ao_basis_nucleus_index(context, &
       idx, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    integer (c_int64_t) , intent(out)         :: idx(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_get_ao_basis_nucleus_index
end interface

interface
  integer(c_int32_t) function qmckl_get_ao_basis_shell_ang_mom(context, &
       shell_ang_mom, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    integer (c_int32_t) , intent(out)         :: shell_ang_mom(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_get_ao_basis_shell_ang_mom
end interface

interface
  integer(c_int32_t) function qmckl_get_ao_basis_shell_prim_num(context, &
       shell_prim_num, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    integer (c_int64_t) , intent(out)         :: shell_prim_num(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_get_ao_basis_shell_prim_num
end interface

interface
  integer(c_int32_t) function qmckl_get_ao_basis_shell_prim_index(context, &
       shell_prim_index, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    integer (c_int64_t) , intent(out)         :: shell_prim_index(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_get_ao_basis_shell_prim_index
end interface

interface
  integer(c_int32_t) function qmckl_get_ao_basis_shell_factor(context, &
       shell_factor, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    real    (c_double)  , intent(out)         :: shell_factor(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_get_ao_basis_shell_factor
end interface

interface
  integer(c_int32_t) function qmckl_get_ao_basis_exponent(context, &
       exponent, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    real    (c_double)  , intent(out)         :: exponent(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_get_ao_basis_exponent
end interface

interface
  integer(c_int32_t) function qmckl_get_ao_basis_coefficient(context, &
       coefficient, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    real    (c_double)  , intent(out)         :: coefficient(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_get_ao_basis_coefficient
end interface

interface
  integer(c_int32_t) function qmckl_get_ao_basis_prim_factor(context, &
       prim_factor, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    real    (c_double)  , intent(out)         :: prim_factor(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_get_ao_basis_prim_factor
end interface

interface
  integer(c_int32_t) function qmckl_get_ao_basis_ao_num(context, &
       num) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    integer (c_int64_t) , intent(out)         :: num
  end function qmckl_get_ao_basis_ao_num
end interface

interface
  integer(c_int32_t) function qmckl_get_ao_basis_cartesian(context, &
       cartesian) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    logical (c_bool)    , intent(out)         :: cartesian
  end function qmckl_get_ao_basis_cartesian
end interface

interface
  integer(c_int32_t) function qmckl_get_ao_basis_ao_factor(context, &
       ao_factor, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    real    (c_double)  , intent(out)         :: ao_factor(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_get_ao_basis_ao_factor
end interface

2.2 Computed data

The following data is computed as described in the next sections:

Variable Type Description
primitive_vgl double[point_num][5][prim_num] Value, gradients, Laplacian of the primitives at current positions
primitive_vgl_date uint64_t Last modification date of Value, gradients, Laplacian of the primitives at current positions
shell_vgl double[point_num][5][shell_num] Value, gradients, Laplacian of the primitives at current positions
shell_vgl_date uint64_t Last modification date of Value, gradients, Laplacian of the AOs at current positions
ao_vgl double[point_num][5][ao_num] Value, gradients, Laplacian of the AOs at current positions
ao_vgl_date uint64_t Last modification date of Value, gradients, Laplacian of the AOs at current positions
ao_value double[point_num][ao_num] Values of the the AOs at current positions
ao_value_date uint64_t Last modification date of the values of the AOs at current positions

2.2.1 After initialization

When the basis set is completely entered, extra data structures may be computed to accelerate the calculations. The primitives within each contraction are sorted in ascending order of their exponents, such that as soon as a primitive is zero all the following functions vanish. Also, it is possible to compute a nuclear radius beyond which all the primitives are zero up to the numerical accuracy defined in the context.

2.2.2 TODO HPC-specific data structures

For faster access, we provide extra arrays for the shell information as:

  • \(C_{psa}\) = coef_per_nucleus[inucl][ishell][iprim]
  • \(\gamma_{pa}\) expo_per_nucleus[inucl][iprim]

such that the computation of the radial parts can be vectorized.

Exponents are sorted in increasing order to exit loops quickly, and only unique exponents are kept. This also allows to pack the exponents to enable vectorization of exponentials.

The computation of the radial part is made as \[ R_{sa} = \sum_p C_{psa} \gamma_{pa} \] which is a matrix-vector product.

2.2.3 Access functions

qmckl_exit_code
qmckl_get_ao_basis_primitive_vgl (qmckl_context context,
                                  double* const primitive_vgl,
                                  const int64_t size_max);

Returns the array of values, gradients an Laplacian of primitive basis functions evaluated at the current coordinates. See section 3.2. 

qmckl_exit_code
qmckl_get_ao_basis_shell_vgl (qmckl_context context,
                              double* const shell_vgl,
                              const int64_t size_max);

Returns the array of values, gradients an Laplacian of contracted shells evaluated at the current coordinates. See section 3.3. 

qmckl_exit_code
qmckl_get_ao_basis_ao_vgl (qmckl_context context,
                           double* const ao_vgl,
                           const int64_t size_max);

Returns the array of values, gradients an Laplacian of the atomic orbitals evaluated at the current coordinates. See section 5.

Uses the given array to compute the VGL. 

qmckl_exit_code
qmckl_get_ao_basis_ao_vgl_inplace (qmckl_context context,
                                   double* const ao_vgl,
                                   const int64_t size_max);

qmckl_exit_code
qmckl_get_ao_basis_ao_value (qmckl_context context,
                             double* const ao_value,
                             const int64_t size_max);

Returns the array of values of the atomic orbitals evaluated at the current coordinates. See section 5.

Uses the given array to compute the value. 

qmckl_exit_code
qmckl_get_ao_basis_ao_value_inplace (qmckl_context context,
                                     double* const ao_value,
                                     const int64_t size_max);

3 Radial part

3.1 General functions for 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 \]

Variable Type Description
context qmckl_context Global state
X(3) double[3] Array containing the coordinates of the points
R(3) double[3] Array containing the x,y,z coordinates of the center
n int64_t Number of computed Gaussians
A(n) double[n] Exponents of the Gaussians
VGL(ldv,5) double[5][ldv] Value, gradients and Laplacian of the Gaussians
ldv int64_t Leading dimension of array VGL

Requirements:

  • context ≠ 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
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);

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

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

  info = QMCKL_SUCCESS

  if (context == QMCKL_NULL_CONTEXT) then
     info = QMCKL_INVALID_CONTEXT
     return
  endif

  if (n <= 0) then
     info = QMCKL_INVALID_ARG_4
     return
  endif

  if (ldv < n) then
     info = QMCKL_INVALID_ARG_7
     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

3.2 Computation of primitives

Variable Type In/Out Description
context qmckl_context in Global state
prim_num int64_t in Number of primitives
point_num int64_t in Number of points considered
nucl_num int64_t in Number of nuclei
nucleus_prim_index int64_t[nucl_num] in Index of the 1st primitive of each nucleus
coord double[3][point_num] in Coordinates
nucl_coord double[3][nucl_num] in Nuclear coordinates
expo double[prim_num] in Exponents of the primitives
primitive_vgl double[point_num][5][prim_num] out Value, gradients and Laplacian of the primitives 
qmckl_exit_code qmckl_compute_ao_basis_primitive_gaussian_vgl (
      const qmckl_context context,
      const int64_t prim_num,
      const int64_t point_num,
      const int64_t nucl_num,
      const int64_t* nucleus_prim_index,
      const double* coord,
      const double* nucl_coord,
      const double* expo,
      double* const primitive_vgl ); 

integer function qmckl_compute_ao_basis_primitive_gaussian_vgl_f( &
     context, prim_num, point_num, nucl_num,                       &
     nucleus_prim_index, coord, nucl_coord,                  &
     expo, primitive_vgl)                                         &
     result(info)

  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer*8             , intent(in)  :: prim_num
  integer*8             , intent(in)  :: nucl_num
  integer*8             , intent(in)  :: point_num
  integer*8             , intent(in)  :: nucleus_prim_index(nucl_num+1)
  double precision      , intent(in)  :: coord(point_num,3)
  double precision      , intent(in)  :: nucl_coord(nucl_num,3)
  double precision      , intent(in)  :: expo(prim_num)
  double precision      , intent(out) :: primitive_vgl(prim_num,5,point_num)

  integer*8 :: inucl, iprim, ipoint
  double precision :: x, y, z, two_a, ar2, r2, v, cutoff

  info = QMCKL_SUCCESS

  ! Don't compute exponentials when the result will be almost zero.
  cutoff = 27.631021115928547d0  ! -dlog(1.d-12)

  do inucl=1,nucl_num
     ! C is zero-based, so shift bounds by one
     do iprim = nucleus_prim_index(inucl)+1, nucleus_prim_index(inucl+1)
        do ipoint = 1, point_num
           x = coord(ipoint,1) - nucl_coord(inucl,1)
           y = coord(ipoint,2) - nucl_coord(inucl,2)
           z = coord(ipoint,3) - nucl_coord(inucl,3)

           r2 = x*x + y*y + z*z
           ar2 = expo(iprim)*r2
           if (ar2 > cutoff) cycle

           v = dexp(-ar2)
           two_a = -2.d0 * expo(iprim) * v

           primitive_vgl(iprim, 1, ipoint) = v
           primitive_vgl(iprim, 2, ipoint) = two_a * x
           primitive_vgl(iprim, 3, ipoint) = two_a * y
           primitive_vgl(iprim, 4, ipoint) = two_a * z
           primitive_vgl(iprim, 5, ipoint) = two_a * (3.d0 - 2.d0*ar2)

        end do
     end do
  end do

end function qmckl_compute_ao_basis_primitive_gaussian_vgl_f

3.3 Computation of shells

Variable Type In/Out Description
context qmckl_context in Global state
prim_num int64_t in Number of primitives
shell_num int64_t in Number of shells
point_num int64_t in Number of points
nucl_num int64_t in Number of nuclei
nucleus_shell_num int64_t[nucl_num] in Number of shells for each nucleus
nucleus_index int64_t[nucl_num] in Index of the 1st shell of each nucleus
nucleus_range double[nucl_num] in Range of the nucleus
shell_prim_index int64_t[shell_num] in Index of the 1st primitive of each shell
shell_prim_num int64_t[shell_num] in Number of primitives per shell
coord double[3][point_num] in Coordinates
nucl_coord double[3][nucl_num] in Nuclear coordinates
expo double[prim_num] in Exponents of the primitives
coef_normalized double[prim_num] in Coefficients of the primitives
shell_vgl double[point_num][5][shell_num] out Value, gradients and Laplacian of the shells 
qmckl_exit_code qmckl_compute_ao_basis_shell_gaussian_vgl (
      const qmckl_context context,
      const int64_t prim_num,
      const int64_t shell_num,
      const int64_t point_num,
      const int64_t nucl_num,
      const int64_t* nucleus_shell_num,
      const int64_t* nucleus_index,
      const double* nucleus_range,
      const int64_t* shell_prim_index,
      const int64_t* shell_prim_num,
      const double* coord,
      const double* nucl_coord,
      const double* expo,
      const double* coef_normalized,
      double* const shell_vgl ); 

integer function qmckl_compute_ao_basis_shell_gaussian_vgl_f( &
     context, prim_num, shell_num, point_num, nucl_num,       &
     nucleus_shell_num, nucleus_index, nucleus_range,         &
     shell_prim_index, shell_prim_num, coord, nucl_coord,     &
     expo, coef_normalized, shell_vgl)                        &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer*8             , intent(in)  :: prim_num
  integer*8             , intent(in)  :: shell_num
  integer*8             , intent(in)  :: nucl_num
  integer*8             , intent(in)  :: point_num
  integer*8             , intent(in)  :: nucleus_shell_num(nucl_num)
  integer*8             , intent(in)  :: nucleus_index(nucl_num)
  double precision      , intent(in)  :: nucleus_range(nucl_num)
  integer*8             , intent(in)  :: shell_prim_index(shell_num)
  integer*8             , intent(in)  :: shell_prim_num(shell_num)
  double precision      , intent(in)  :: coord(point_num,3)
  double precision      , intent(in)  :: nucl_coord(nucl_num,3)
  double precision      , intent(in)  :: expo(prim_num)
  double precision      , intent(in)  :: coef_normalized(prim_num)
  double precision      , intent(out) :: shell_vgl(shell_num,5,point_num)

  integer*8 :: inucl, iprim, ipoint, ishell
  integer*8 :: ishell_start, ishell_end
  integer*8 :: iprim_start , iprim_end
  double precision :: x, y, z, two_a, ar2, r2, v, cutoff

  info = QMCKL_SUCCESS

  ! Don't compute exponentials when the result will be almost zero.
  ! TODO : Use numerical precision here
  cutoff = 27.631021115928547d0  !-dlog(1.d-12)

  do ipoint = 1, point_num

     do inucl=1,nucl_num

        x = coord(ipoint,1) - nucl_coord(inucl,1)
        y = coord(ipoint,2) - nucl_coord(inucl,2)
        z = coord(ipoint,3) - nucl_coord(inucl,3)

        r2 = x*x + y*y + z*z

        if (r2 > cutoff*nucleus_range(inucl)) then
           cycle
        end if

        ! C is zero-based, so shift bounds by one
        ishell_start = nucleus_index(inucl) + 1
        ishell_end   = nucleus_index(inucl) + nucleus_shell_num(inucl)

        do ishell=ishell_start, ishell_end

           shell_vgl(ishell, 1, ipoint) = 0.d0
           shell_vgl(ishell, 2, ipoint) = 0.d0
           shell_vgl(ishell, 3, ipoint) = 0.d0
           shell_vgl(ishell, 4, ipoint) = 0.d0
           shell_vgl(ishell, 5, ipoint) = 0.d0

           iprim_start = shell_prim_index(ishell) + 1
           iprim_end   = shell_prim_index(ishell) + shell_prim_num(ishell)

           do iprim = iprim_start, iprim_end

              ar2 = expo(iprim)*r2
              if (ar2 > cutoff) then
                 cycle
              end if

              v = coef_normalized(iprim) * dexp(-ar2)
              two_a = -2.d0 * expo(iprim) * v

              shell_vgl(ishell, 1, ipoint) = &
                   shell_vgl(ishell, 1, ipoint) + v

              shell_vgl(ishell, 2, ipoint) = &
                   shell_vgl(ishell, 2, ipoint) + two_a * x

              shell_vgl(ishell, 3, ipoint) = &
                   shell_vgl(ishell, 3, ipoint) + two_a * y

              shell_vgl(ishell, 4, ipoint) = &
                   shell_vgl(ishell, 4, ipoint) + two_a * z

              shell_vgl(ishell, 5, ipoint) = &
                   shell_vgl(ishell, 5, ipoint) + two_a * (3.d0 - 2.d0*ar2)

           end do

        end do
     end do

  end do

end function qmckl_compute_ao_basis_shell_gaussian_vgl_f

4 Polynomial part

Going from the atomic basis set to AOs implies a systematic construction of all the angular functions of each shell. We consider two cases for the angular functions: the real-valued spherical harmonics, and the polynomials in Cartesian coordinates. In the case of spherical harmonics, the AOs are ordered in increasing magnetic quantum number (\(-l \le m \le l\)), and in the case of polynomials we choose the canonical ordering, i.e

\begin{eqnarray} p & : & p_x, p_y, p_z \nonumber \\ d & : & d_{xx}, d_{xy}, d_{xz}, d_{yy}, d_{yz}, d_{zz} \nonumber \\ f & : & f_{xxx}, f_{xxy}, f_{xxz}, f_{xyy}, f_{xyz}, f_{xzz}, f_{yyy}, f_{yyz}, f_{yzz}, f_{zzz} \nonumber \\ {\rm etc.} \nonumber \end{eqnarray}

4.1 General functions for Powers of \(x-X_i\)

The qmckl_ao_power function 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_{ik} = X_i^k \]

Variable Type In/Out Description
context qmckl_context in Global state
n int64t in Number of values
X double[n] in Array containing the input values
LMAX int32t[n] in Array containing the maximum power for each value
P double[n][ldp] out Array containing all the powers of X
ldp int64t in Leading dimension of array P

Requirements:

  • context is not QMCKL_NULL_CONTEXT
  • 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]
qmckl_exit_code qmckl_ao_power (
      const qmckl_context context,
      const int64_t n,
      const double* X,
      const int32_t* LMAX,
      double* const P,
      const int64_t ldp ); 

integer function qmckl_ao_power_f(context, n, X, LMAX, P, ldp) result(info)
  use qmckl
  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,k

  info = QMCKL_SUCCESS

  if (context == QMCKL_NULL_CONTEXT) then
     info = QMCKL_INVALID_CONTEXT
     return
  endif

  if (n <= ldp) then
     info = QMCKL_INVALID_ARG_2
     return
  endif

  k = MAXVAL(LMAX)
  if (LDP < k) then
     info = QMCKL_INVALID_ARG_6
     return
  endif

  if (k <= 0) then
     info = QMCKL_INVALID_ARG_4
     return
  endif

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

end function qmckl_ao_power_f

4.2 General functions for Value, Gradient and Laplacian of a polynomial

A polynomial is centered on a nucleus \(\mathbf{R}_i\)

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

The gradients with respect to electron coordinates are

\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*}

and the Laplacian is

\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_polynomial_vgl computes the values, gradients and Laplacians at a given point in space, of all polynomials with an angular momentum up to lmax.

Variable Type In/Out Description
context qmckl_context in Global state
X double[3] in Array containing the coordinates of the points
R double[3] in Array containing the x,y,z coordinates of the center
lmax int32_t in Maximum angular momentum
n int64_t inout Number of computed polynomials
L int32_t[n][ldl] out Contains a,b,c for all n results
ldl int64_t in Leading dimension of L
VGL double[n][ldv] out Value, gradients and Laplacian of the polynomials
ldv int64_t in Leading dimension of array VGL

Requirements:

  • context ≠ QMCKL_NULL_CONTEXT
  • 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):
    • Increasing 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"
qmckl_exit_code qmckl_ao_polynomial_vgl (
      const qmckl_context context,
      const double* X,
      const double* R,
      const int32_t lmax,
      int64_t* n,
      int32_t* const L,
      const int64_t ldl,
      double* const VGL,
      const int64_t ldv ); 

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

qmckl_exit_code
qmckl_ao_polynomial_vgl (const qmckl_context context,
                         const double* X,
                         const double* R,
                         const int32_t lmax,
                         int64_t* n,
                         int32_t* const L,
                         const int64_t ldl,
                         double* const VGL,
                         const int64_t ldv )
{
#ifdef HAVE_HPC
  //return qmckl_ao_polynomial_vgl_hpc (context, X, R, lmax, n, L, ldl, VGL, ldv);
  return qmckl_ao_polynomial_vgl_doc (context, X, R, lmax, n, L, ldl, VGL, ldv);
#else
  return qmckl_ao_polynomial_vgl_doc (context, X, R, lmax, n, L, ldl, VGL, ldv);
#endif
}

integer function qmckl_ao_polynomial_vgl_doc_f (context,    &
     X, R, lmax, n, L, ldl, VGL, ldv) result(info)
  use qmckl
  implicit none
  integer*8,        intent(in)  :: context
  double precision, 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
  double precision, 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
  double precision  :: Y(3)
  double precision  :: pows(-2:lmax,3)
  double precision  :: xy, yz, xz
  double precision  :: da, db, dc, dd

  info = 0

  if (context == QMCKL_NULL_CONTEXT) then
     info = QMCKL_INVALID_CONTEXT
     return
  endif

  if (lmax < 0) then
     info = QMCKL_INVALID_ARG_4
     return
  endif

  if (ldl < 3) then
     info = QMCKL_INVALID_ARG_7
     return
  endif

  if (ldv < 5) then
     info = QMCKL_INVALID_ARG_9
     return
  endif


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

  if (lmax == 0) then
     VGL(1,1) = 1.d0
     VGL(2,1) = 0.d0
     VGL(3,1) = 0.d0
     VGL(4,1) = 0.d0
     VGL(5,1) = 0.d0
     l(1,1) = 0
     l(2,1) = 0
     l(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
     VGL(1,1) = 1.d0
     VGL(1,2) = pows(1,1)
     VGL(2,2) = 1.d0
     VGL(1,3) = pows(1,2)
     VGL(3,3) = 1.d0
     VGL(1,4) = pows(1,3)
     VGL(4,4) = 1.d0

     l  (1:3,1:4) = 0
     l  (1,2) = 1
     l  (2,3) = 1
     l  (3,4) = 1

     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 = QMCKL_SUCCESS

end function qmckl_ao_polynomial_vgl_doc_f

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

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

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

qmckl_exit_code
qmckl_ao_polynomial_transp_vgl (const qmckl_context context,
                                const double* X,
                                const double* R,
                                const int32_t lmax,
                                int64_t* n,
                                int32_t* const L,
                                const int64_t ldl,
                                double* const VGL,
                                const int64_t ldv )
{
#ifdef HAVE_HPC
  return qmckl_ao_polynomial_transp_vgl_hpc (context, X, R, lmax, n, L, ldl, VGL, ldv);
#else
  return qmckl_ao_polynomial_transp_vgl_doc (context, X, R, lmax, n, L, ldl, VGL, ldv);
#endif
}

integer function qmckl_ao_polynomial_transp_vgl_doc_f (context,    &
     X, R, lmax, n, L, ldl, VGL, ldv) result(info)
  use qmckl
  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,5)
  integer*8 , intent(in)  :: ldv

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

  info = 0

  if (context == QMCKL_NULL_CONTEXT) then
     info = QMCKL_INVALID_CONTEXT
     return
  endif

  if (lmax < 0) then
     info = QMCKL_INVALID_ARG_4
     return
  endif

  if (ldl < 3) then
     info = QMCKL_INVALID_ARG_7
     return
  endif

  if (ldv < (lmax+1)*(lmax+2)*(lmax+3)/6) then
     info = QMCKL_INVALID_ARG_9
     return
  endif


  if (lmax > 0) then

     do i=1,3
        Y(i) = X(i) - R(i)
     end do
     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

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

     VGL(1  ,1  ) = 1.d0

     l  (1,2) = 1
     VGL(2,1) = Y(1)
     VGL(2,2) = 1.d0

     l  (2,3) = 1
     VGL(3,1) = Y(2)
     VGL(3,3) = 1.d0

     l  (3,4) = 1
     VGL(4,1) = Y(3)
     VGL(4,4) = 1.d0

     n=4
  else
     VGL(1,1) = 1.d0
     VGL(1,2:5) = 0.d0
     l(1:3,1) = 0
     n=1
     return
  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

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

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

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

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

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

           VGL(n,5) = &
                (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 = QMCKL_SUCCESS

end function qmckl_ao_polynomial_transp_vgl_doc_f

static inline qmckl_exit_code
qmckl_ao_polynomial_transp_vgl_hpc_inline (const qmckl_context context,
                                    const double* restrict X,
                                    const double* restrict R,
                                    const int32_t lmax,
                                    int64_t* n,
                                    int32_t* restrict const L,
                                    const int64_t ldl,
                                    double* restrict const VGL,
                                    const int64_t ldv )
{
  const qmckl_context_struct* ctx = (qmckl_context_struct*) context;
  assert (ctx != NULL && X != NULL && R != NULL && n != NULL && L != NULL && VGL != NULL);
  if (lmax < 0) return QMCKL_INVALID_ARG_4;
  if (ldl < 3) return QMCKL_INVALID_ARG_7;

  int32_t m;

  switch (lmax) {
  case 0:
    {
      if (ldv < 1) return QMCKL_INVALID_ARG_9;
      L[0] = 0; L[1] = 0; L[2] = 0;

      VGL[0         ] = 1.0;
      VGL[ldv       ] = 0.0;
      VGL[ldv<<1    ] = 0.0;
      VGL[(ldv<<1)+ldv] = 0.0;
      VGL[ldv<<2    ] = 0.0;
      m=1;
      break;
    }
  case 1:
    {
      if (ldv < 4) return QMCKL_INVALID_ARG_9;

      if (ldl == 3) {

        const int32_t ltmp[12] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
        for (int i=0 ; i<12 ; ++i)
          L[i] = ltmp[i];

      } else {

        int32_t* restrict const l[4] = {L, L+ldl, L+(ldl<<1), L+ldl+(ldl<<1)};
        l[0][0] = 0; l[0][1] = 0; l[0][2] = 0;
        l[1][0] = 1; l[1][1] = 0; l[1][2] = 0;
        l[2][0] = 0; l[2][1] = 1; l[2][2] = 0;
        l[3][0] = 0; l[3][1] = 0; l[3][2] = 1;
      }


      if (ldv == 4) {

        const double tmp[20] = {1.0, X[0]-R[0], X[1]-R[1], X[2]-R[2],
                                0.0, 1.0, 0.0, 0.0,
                                0.0, 0.0, 1.0, 0.0,
                                0.0, 0.0, 0.0, 1.0,
                                0.0, 0.0, 0.0, 0.0};

        for (int i=0 ; i<20 ; ++i)
          VGL[i] = tmp[i];

      } else {

        double* restrict const vgl1 = VGL;
        double* restrict const vgl2 = VGL + ldv;
        double* restrict const vgl3 = VGL + (ldv << 1);
        double* restrict const vgl4 = VGL + ldv + (ldv << 1);
        double* restrict const vgl5 = VGL + (ldv << 2);

        for (int32_t k=0 ; k<4 ; ++k) {
          vgl2[k] = 0.0;
          vgl3[k] = 0.0;
          vgl4[k] = 0.0;
          vgl5[k] = 0.0;
        }
        vgl1[0] = 1.0;
        vgl1[1] = X[0]-R[0];
        vgl1[2] = X[1]-R[1];
        vgl1[3] = X[2]-R[2];
        vgl2[1] = 1.0;
        vgl3[2] = 1.0;
        vgl4[3] = 1.0;
      }
      m=4;
      break;
    }
  case 2:
    {
      if (ldv < 10) return QMCKL_INVALID_ARG_9;
      if (ldl == 3) {
        const int32_t ltmp[30] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1,
                                  2, 0, 0, 1, 1, 0, 1, 0, 1, 0, 2, 0,
                                  0, 1, 1, 0, 0, 2};
        for (int i=0 ; i<30 ; ++i)
          L[i] = ltmp[i];

      } else {
        int32_t* restrict l[10];
        for (int32_t i=0 ; i<10 ; ++i) {
          l[i] = L + i*ldl;
        }

        l[0][0] = 0; l[0][1] = 0; l[0][2] = 0;
        l[1][0] = 1; l[1][1] = 0; l[1][2] = 0;
        l[2][0] = 0; l[2][1] = 1; l[2][2] = 0;
        l[3][0] = 0; l[3][1] = 0; l[3][2] = 1;
        l[4][0] = 2; l[4][1] = 0; l[4][2] = 0;
        l[5][0] = 1; l[5][1] = 1; l[5][2] = 0;
        l[6][0] = 1; l[6][1] = 0; l[6][2] = 1;
        l[7][0] = 0; l[7][1] = 2; l[7][2] = 0;
        l[8][0] = 0; l[8][1] = 1; l[8][2] = 1;
        l[9][0] = 0; l[9][1] = 0; l[9][2] = 2;
      }

      const double Y[3] = { X[0]-R[0],
                            X[1]-R[1],
                            X[2]-R[2] };

      if (ldv == 50) {
         const double tmp[50] = {
            1.0, Y[0], Y[1], Y[2], Y[0] * Y[0],
            Y[0] * Y[1], Y[0] * Y[2], Y[1] * Y[1], Y[1] * Y[2], Y[2] * Y[2],
            0.0, 1.0, 0.0, 0.0, Y[0] + Y[0],
            Y[1], Y[2], 0.0, 0.0, 0.0,
            0.0, 0.0, 1.0, 0.0, 0.0,
            Y[0], 0.0, Y[1] + Y[1], Y[2], 0.0,
            0.0, 0.0, 0.0, 1.0, 0.0,
            0.0, Y[0], 0.0, Y[1], Y[2] + Y[2],
            0.0, 0.0, 0.0, 0.0, 2.0,
            0.0, 0.0, 2.0, 0., 2.0 };
         for (int i=0 ; i<50 ; ++i)
           VGL[i] = tmp[i];
      } else {
        double* restrict const vgl1 = VGL;
        double* restrict const vgl2 = VGL + ldv;
        double* restrict const vgl3 = VGL + (ldv << 1);
        double* restrict const vgl4 = VGL + 3*ldv;
        double* restrict const vgl5 = VGL + (ldv << 2);

        vgl1[0] = 1.0      ; vgl1[1] = Y[0]     ; vgl1[2] = Y[1];
        vgl1[3] = Y[2]     ; vgl1[4] = Y[0]*Y[0]; vgl1[5] = Y[0]*Y[1];
        vgl1[6] = Y[0]*Y[2]; vgl1[7] = Y[1]*Y[1]; vgl1[8] = Y[1]*Y[2];
        vgl1[9] = Y[2]*Y[2];

        vgl2[0] = 0.0 ; vgl2[1] = 1.0      ; vgl2[2] = 0.0 ;
        vgl2[3] = 0.0 ; vgl2[4] = Y[0]+Y[0]; vgl2[5] = Y[1];
        vgl2[6] = Y[2]; vgl2[7] = 0.0      ; vgl2[8] = 0.0 ;
        vgl2[9] = 0.0 ;

        vgl3[0] = 0.0; vgl3[1] = 0.0      ; vgl3[2] = 1.0 ;
        vgl3[3] = 0.0; vgl3[4] = 0.0      ; vgl3[5] = Y[0];
        vgl3[6] = 0.0; vgl3[7] = Y[1]+Y[1]; vgl3[8] = Y[2];
        vgl3[9] = 0.0;

        vgl4[0] = 0.0      ; vgl4[1] = 0.0; vgl4[2] = 0.0 ;
        vgl4[3] = 1.0      ; vgl4[4] = 0.0; vgl4[5] = 0.0 ;
        vgl4[6] = Y[0]     ; vgl4[7] = 0.0; vgl4[8] = Y[1];
        vgl4[9] = Y[2]+Y[2];

        vgl5[0] = 0.0; vgl5[1] = 0.0; vgl5[2] = 0.0;
        vgl5[3] = 0.0; vgl5[4] = 2.0; vgl5[5] = 0.0;
        vgl5[6] = 0.0; vgl5[7] = 2.0; vgl5[8] = 0.0;
        vgl5[9] = 2.0;
      }
      m=10;
      break;
    }
  default:
    {
      const int32_t size_max = (lmax+1)*(lmax+2)*(lmax+3)/6;
      if (ldv < size_max) return QMCKL_INVALID_ARG_9;

      double* restrict const vgl1 = VGL;
      double* restrict const vgl2 = VGL + ldv;
      double* restrict const vgl3 = VGL + (ldv<<1);
      double* restrict const vgl4 = VGL + ldv + (ldv<<1);
      double* restrict const vgl5 = VGL + (ldv<<2);

      const double Y[3] = { X[0]-R[0],
                            X[1]-R[1],
                            X[2]-R[2] };

      assert(size_max > lmax+3);
      double pows[3][size_max];

      for (int32_t i=0 ; i<3 ; ++i) {
        pows[0][i] = 1.0;
        pows[1][i] = 1.0;
        pows[2][i] = 1.0;
      }

      for (int32_t i=3 ; i<=lmax+2 ; ++i) {
        pows[0][i] = pows[0][i-1] * Y[0];
        pows[1][i] = pows[1][i-1] * Y[1];
        pows[2][i] = pows[2][i-1] * Y[2];
      }

      int32_t* l[size_max];
      for (int32_t i=0 ; i<size_max ; ++i) {
        l[i] = &(L[i*ldl]);
      }

      for (int32_t i=0 ; i<4 ; ++i) {
        l[i][0] = 0;
        l[i][1] = 0;
        l[i][2] = 0;
      }
      l[1][0] = 1;
      l[2][1] = 1;
      l[3][2] = 1;

      for (int32_t k=0 ; k<4 ; ++k) {
        vgl2[k] = 0.0;
        vgl3[k] = 0.0;
        vgl4[k] = 0.0;
        vgl5[k] = 0.0;
      }
      vgl1[0] = 1.0;
      vgl1[1] = Y[0];
      vgl1[2] = Y[1];
      vgl1[3] = Y[2];

      vgl2[1] = 1.0;
      vgl3[2] = 1.0;
      vgl4[3] = 1.0;
      m=4;

      double dd = 2.0;

      for (int32_t d=2 ; d<= lmax ; ++d) {
        double da = dd;

        for (int32_t a=d ; a>=0 ; --a) {
          double db = dd-da;

          for (int32_t b=d-a ; b>=0 ; --b) {
            const int32_t c  = d - a - b;
            const double dc = dd - da - db;

            double xy = pows[0][a+2] * pows[1][b+2];
            double yz = pows[1][b+2] * pows[2][c+2];
            double xz = pows[0][a+2] * pows[2][c+2];

            l[m][0] = a;
            l[m][1] = b;
            l[m][2] = c;

            vgl1[m] = xy * pows[2][c+2];

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

            vgl2[m] = pows[0][a+1] * yz;
            vgl3[m] = pows[1][b+1] * xz;
            vgl4[m] = pows[2][c+1] * xy;

            vgl5[m] = (da-1.) * pows[0][a] * yz +
              (db-1.) * pows[1][b] * xz +
              (dc-1.) * pows[2][c] * xy;

            db -= 1.0;
            ++m;
          }
          da -= 1.0;
        }
        dd += 1.0;
      }
    }
  }
  *n = m;
  return QMCKL_SUCCESS;
}

qmckl_exit_code
qmckl_ao_polynomial_transp_vgl_hpc (const qmckl_context context,
                                    const double* restrict X,
                                    const double* restrict R,
                                    const int32_t lmax,
                                    int64_t* n,
                                    int32_t* restrict const L,
                                    const int64_t ldl,
                                    double* restrict const VGL,
                                    const int64_t ldv )
{
  return qmckl_ao_polynomial_transp_vgl_hpc_inline (context,
                                    X, R, lmax, n, L, ldl, VGL, ldv );
  }

5 Combining radial and polynomial parts

5.1 Values only

5.1.1 Unoptimized version

Variable Type In/Out Description
context qmckl_context in Global state
ao_num int64_t in Number of AOs
shell_num int64_t in Number of shells
point_num int64_t in Number of points
nucl_num int64_t in Number of nuclei
coord double[3][point_num] in Coordinates
nucl_coord double[3][nucl_num] in Nuclear coordinates
nucleus_index int64_t[nucl_num] in Index of the 1st shell of each nucleus
nucleus_shell_num int64_t[nucl_num] in Number of shells per nucleus
nucleus_range double[nucl_num] in Range beyond which all is zero
nucleus_max_ang_mom int32_t[nucl_num] in Maximum angular momentum per nucleus
shell_ang_mom int32_t[shell_num] in Angular momentum of each shell
ao_factor double[ao_num] in Normalization factor of the AOs
shell_vgl double[point_num][5][shell_num] in Value, gradients and Laplacian of the shells
ao_value double[point_num][ao_num] out Values of the AOs 
integer function qmckl_compute_ao_value_doc_f(context, &
     ao_num, shell_num, point_num, nucl_num, &
     coord, nucl_coord, nucleus_index, nucleus_shell_num, &
     nucleus_range, nucleus_max_ang_mom, shell_ang_mom, &
     ao_factor, shell_vgl, ao_value) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer*8             , intent(in)  :: ao_num
  integer*8             , intent(in)  :: shell_num
  integer*8             , intent(in)  :: point_num
  integer*8             , intent(in)  :: nucl_num
  double precision      , intent(in)  :: coord(point_num,3)
  double precision      , intent(in)  :: nucl_coord(nucl_num,3)
  integer*8             , intent(in)  :: nucleus_index(nucl_num)
  integer*8             , intent(in)  :: nucleus_shell_num(nucl_num)
  double precision      , intent(in)  :: nucleus_range(nucl_num)
  integer               , intent(in)  :: nucleus_max_ang_mom(nucl_num)
  integer               , intent(in)  :: shell_ang_mom(shell_num)
  double precision      , intent(in)  :: ao_factor(ao_num)
  double precision      , intent(in)  :: shell_vgl(shell_num,5,point_num)
  double precision      , intent(out) :: ao_value(ao_num,point_num)

  double precision  :: e_coord(3), n_coord(3)
  integer*8         :: n_poly
  integer           :: l, il, k
  integer*8         :: ipoint, inucl, ishell
  integer*8         :: ishell_start, ishell_end
  integer           :: lstart(0:20)
  double precision  :: x, y, z, r2
  double precision  :: cutoff
  integer, external :: qmckl_ao_polynomial_vgl_doc_f

  double precision, allocatable  :: poly_vgl(:,:)
  integer         , allocatable  :: powers(:,:), ao_index(:)

  allocate(poly_vgl(5,ao_num), powers(3,ao_num), ao_index(ao_num))

  ! Pre-computed data
  do l=0,20
     lstart(l) = l*(l+1)*(l+2)/6 +1
  end do

  k=1
  do inucl=1,nucl_num
     ishell_start = nucleus_index(inucl) + 1
     ishell_end   = nucleus_index(inucl) + nucleus_shell_num(inucl)
     do ishell = ishell_start, ishell_end
        l = shell_ang_mom(ishell)
        ao_index(ishell) = k
        k = k + lstart(l+1) - lstart(l)
     end do
  end do
  info = QMCKL_SUCCESS

  ! Don't compute polynomials when the radial part is zero.
  cutoff = 27.631021115928547d0  !-dlog(1.d-12)

  do ipoint = 1, point_num
     e_coord(1) = coord(ipoint,1)
     e_coord(2) = coord(ipoint,2)
     e_coord(3) = coord(ipoint,3)
     ao_value(:,ipoint) = 0.d0
     do inucl=1,nucl_num
        n_coord(1) = nucl_coord(inucl,1)
        n_coord(2) = nucl_coord(inucl,2)
        n_coord(3) = nucl_coord(inucl,3)

        ! Test if the point is in the range of the nucleus
        x = e_coord(1) - n_coord(1)
        y = e_coord(2) - n_coord(2)
        z = e_coord(3) - n_coord(3)

        r2 = x*x + y*y + z*z

        if (r2 > cutoff*nucleus_range(inucl)) then
           cycle
        end if

        ! Compute polynomials
        info = qmckl_ao_polynomial_vgl_doc_f(context, e_coord, n_coord, &
             nucleus_max_ang_mom(inucl), n_poly, powers, 3_8, &
             poly_vgl, 5_8)

        ! Loop over shells
        ishell_start = nucleus_index(inucl) + 1
        ishell_end   = nucleus_index(inucl) + nucleus_shell_num(inucl)
        do ishell = ishell_start, ishell_end
           k = ao_index(ishell)
           l = shell_ang_mom(ishell)
           do il = lstart(l), lstart(l+1)-1
              ! Value
              ao_value(k,ipoint) = &
                   poly_vgl(1,il) * shell_vgl(ishell,1,ipoint) * ao_factor(k)
              k = k+1
           end do
        end do
     end do
  end do

  deallocate(poly_vgl, powers)
end function qmckl_compute_ao_value_doc_f

5.1.2 HPC version

Variable Type In/Out Description
context qmckl_context in Global state
ao_num int64_t in Number of AOs
shell_num int64_t in Number of shells
prim_num int64_t in Number of primitives
point_num int64_t in Number of points
nucl_num int64_t in Number of nuclei
coord double[3][point_num] in Coordinates
nucl_coord double[3][nucl_num] in Nuclear coordinates
nucleus_index int64_t[nucl_num] in Index of the 1st shell of each nucleus
nucleus_shell_num int64_t[nucl_num] in Number of shells per nucleus
nucleus_range double[nucl_num] in Range beyond which all is zero
nucleus_max_ang_mom int32_t[nucl_num] in Maximum angular momentum per nucleus
shell_ang_mom int32_t[shell_num] in Angular momentum of each shell
shell_prim_index int64_t[shell_num] in Index of the 1st primitive of each shell
shell_prim_num int64_t[shell_num] in Number of primitives per shell
ao_factor double[ao_num] in Normalization factor of the AOs
ao_expo double[prim_num] in Value, gradients and Laplacian of the shells
coef_normalized double[prim_num] in Value, gradients and Laplacian of the shells
ao_value double[point_num][ao_num] out Values of the AOs

5.1.3 Interfaces

qmckl_exit_code qmckl_compute_ao_value_doc (
      const qmckl_context context,
      const int64_t ao_num,
      const int64_t shell_num,
      const int64_t point_num,
      const int64_t nucl_num,
      const double* coord,
      const double* nucl_coord,
      const int64_t* nucleus_index,
      const int64_t* nucleus_shell_num,
      const double* nucleus_range,
      const int32_t* nucleus_max_ang_mom,
      const int32_t* shell_ang_mom,
      const double* ao_factor,
      const double* shell_vgl,
      double* const ao_value );

#ifdef HAVE_HPC
    qmckl_exit_code qmckl_compute_ao_value_hpc_gaussian (
          const qmckl_context context,
          const int64_t ao_num,
          const int64_t shell_num,
          const int32_t* prim_num_per_nucleus,
          const int64_t point_num,
          const int64_t nucl_num,
          const double* coord,
          const double* nucl_coord,
          const int64_t* nucleus_index,
          const int64_t* nucleus_shell_num,
          const double* nucleus_range,
          const int32_t* nucleus_max_ang_mom,
          const int32_t* shell_ang_mom,
          const double* ao_factor,
          const qmckl_matrix expo_per_nucleus,
          const qmckl_tensor coef_per_nucleus,
          double* const ao_value );
#endif

5.2 Value, gradients, Laplacian

5.2.1 Reference version

Variable Type In/Out Description
context qmckl_context in Global state
ao_num int64_t in Number of AOs
shell_num int64_t in Number of shells
point_num int64_t in Number of points
nucl_num int64_t in Number of nuclei
coord double[3][point_num] in Coordinates
nucl_coord double[3][nucl_num] in Nuclear coordinates
nucleus_index int64_t[nucl_num] in Index of the 1st shell of each nucleus
nucleus_shell_num int64_t[nucl_num] in Number of shells per nucleus
nucleus_range double[nucl_num] in Range beyond which all is zero
nucleus_max_ang_mom int32_t[nucl_num] in Maximum angular momentum per nucleus
shell_ang_mom int32_t[shell_num] in Angular momentum of each shell
ao_factor double[ao_num] in Normalization factor of the AOs
shell_vgl double[point_num][5][shell_num] in Value, gradients and Laplacian of the shells
ao_vgl double[point_num][5][ao_num] out Value, gradients and Laplacian of the AOs 
integer function qmckl_compute_ao_vgl_doc_f(context, &
     ao_num, shell_num, point_num, nucl_num, &
     coord, nucl_coord, nucleus_index, nucleus_shell_num, &
     nucleus_range, nucleus_max_ang_mom, shell_ang_mom, &
     ao_factor, shell_vgl, ao_vgl) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer*8             , intent(in)  :: ao_num
  integer*8             , intent(in)  :: shell_num
  integer*8             , intent(in)  :: point_num
  integer*8             , intent(in)  :: nucl_num
  double precision      , intent(in)  :: coord(point_num,3)
  double precision      , intent(in)  :: nucl_coord(nucl_num,3)
  integer*8             , intent(in)  :: nucleus_index(nucl_num)
  integer*8             , intent(in)  :: nucleus_shell_num(nucl_num)
  double precision      , intent(in)  :: nucleus_range(nucl_num)
  integer               , intent(in)  :: nucleus_max_ang_mom(nucl_num)
  integer               , intent(in)  :: shell_ang_mom(shell_num)
  double precision      , intent(in)  :: ao_factor(ao_num)
  double precision      , intent(in)  :: shell_vgl(shell_num,5,point_num)
  double precision      , intent(out) :: ao_vgl(ao_num,5,point_num)

  double precision  :: e_coord(3), n_coord(3)
  integer*8         :: n_poly
  integer           :: l, il, k
  integer*8         :: ipoint, inucl, ishell
  integer*8         :: ishell_start, ishell_end
  integer           :: lstart(0:20)
  double precision  :: x, y, z, r2
  double precision  :: cutoff
  integer, external :: qmckl_ao_polynomial_vgl_doc_f

  double precision, allocatable  :: poly_vgl(:,:)
  integer         , allocatable  :: powers(:,:), ao_index(:)

  allocate(poly_vgl(5,ao_num), powers(3,ao_num), ao_index(ao_num))

  ! Pre-computed data
  do l=0,20
     lstart(l) = l*(l+1)*(l+2)/6 +1
  end do

  k=1
  do inucl=1,nucl_num
     ishell_start = nucleus_index(inucl) + 1
     ishell_end   = nucleus_index(inucl) + nucleus_shell_num(inucl)
     do ishell = ishell_start, ishell_end
        l = shell_ang_mom(ishell)
        ao_index(ishell) = k
        k = k + lstart(l+1) - lstart(l)
     end do
  end do
  info = QMCKL_SUCCESS

  ! Don't compute polynomials when the radial part is zero.
  cutoff = 27.631021115928547d0  ! -dlog(1.d-12)

  do ipoint = 1, point_num
     e_coord(1) = coord(ipoint,1)
     e_coord(2) = coord(ipoint,2)
     e_coord(3) = coord(ipoint,3)
     ao_vgl(:,:,ipoint) = 0.d0
     do inucl=1,nucl_num
        n_coord(1) = nucl_coord(inucl,1)
        n_coord(2) = nucl_coord(inucl,2)
        n_coord(3) = nucl_coord(inucl,3)

        ! Test if the point is in the range of the nucleus
        x = e_coord(1) - n_coord(1)
        y = e_coord(2) - n_coord(2)
        z = e_coord(3) - n_coord(3)

        r2 = x*x + y*y + z*z

        if (r2 > cutoff*nucleus_range(inucl)) then
           cycle
        end if

        ! Compute polynomials
        info = qmckl_ao_polynomial_vgl_doc_f(context, e_coord, n_coord, &
             nucleus_max_ang_mom(inucl), n_poly, powers, 3_8, &
             poly_vgl, 5_8)

        ! Loop over shells
        ishell_start = nucleus_index(inucl) + 1
        ishell_end   = nucleus_index(inucl) + nucleus_shell_num(inucl)
        do ishell = ishell_start, ishell_end
           k = ao_index(ishell)
           l = shell_ang_mom(ishell)
           do il = lstart(l), lstart(l+1)-1
              ! Value
              ao_vgl(k,1,ipoint) = &
                   poly_vgl(1,il) * shell_vgl(ishell,1,ipoint) * ao_factor(k)

              ! Grad_x
              ao_vgl(k,2,ipoint) = ( &
                   poly_vgl(2,il) * shell_vgl(ishell,1,ipoint) + &
                   poly_vgl(1,il) * shell_vgl(ishell,2,ipoint) &
                   ) * ao_factor(k)

              ! Grad_y
              ao_vgl(k,3,ipoint) = ( &
                   poly_vgl(3,il) * shell_vgl(ishell,1,ipoint) + &
                   poly_vgl(1,il) * shell_vgl(ishell,3,ipoint) &
                   ) * ao_factor(k)

              ! Grad_z
              ao_vgl(k,4,ipoint) = ( &
                   poly_vgl(4,il) * shell_vgl(ishell,1,ipoint) + &
                   poly_vgl(1,il) * shell_vgl(ishell,4,ipoint) &
                   ) * ao_factor(k)

              ! Lapl_z
              ao_vgl(k,5,ipoint) = ( &
                   poly_vgl(5,il) * shell_vgl(ishell,1,ipoint) + &
                   poly_vgl(1,il) * shell_vgl(ishell,5,ipoint) + &
                   2.d0 * ( &
                   poly_vgl(2,il) * shell_vgl(ishell,2,ipoint) + &
                   poly_vgl(3,il) * shell_vgl(ishell,3,ipoint) + &
                   poly_vgl(4,il) * shell_vgl(ishell,4,ipoint) ) &
                   ) * ao_factor(k)

              k = k+1
           end do
        end do
     end do
  end do

  deallocate(poly_vgl, powers)
end function qmckl_compute_ao_vgl_doc_f

5.2.2 HPC version

Variable Type In/Out Description
context qmckl_context in Global state
ao_num int64_t in Number of AOs
shell_num int64_t in Number of shells
prim_num int64_t in Number of primitives
point_num int64_t in Number of points
nucl_num int64_t in Number of nuclei
coord double[3][point_num] in Coordinates
nucl_coord double[3][nucl_num] in Nuclear coordinates
nucleus_index int64_t[nucl_num] in Index of the 1st shell of each nucleus
nucleus_shell_num int64_t[nucl_num] in Number of shells per nucleus
nucleus_range double[nucl_num] in Range beyond which all is zero
nucleus_max_ang_mom int32_t[nucl_num] in Maximum angular momentum per nucleus
shell_ang_mom int32_t[shell_num] in Angular momentum of each shell
shell_prim_index int64_t[shell_num] in Index of the 1st primitive of each shell
shell_prim_num int64_t[shell_num] in Number of primitives per shell
ao_factor double[ao_num] in Normalization factor of the AOs
ao_expo double[prim_num] in Value, gradients and Laplacian of the shells
coef_normalized double[prim_num] in Value, gradients and Laplacian of the shells
ao_vgl double[point_num][5][ao_num] out Value, gradients and Laplacian of the AOs

5.2.3 Interfaces

qmckl_exit_code qmckl_compute_ao_vgl_doc (
      const qmckl_context context,
      const int64_t ao_num,
      const int64_t shell_num,
      const int64_t point_num,
      const int64_t nucl_num,
      const double* coord,
      const double* nucl_coord,
      const int64_t* nucleus_index,
      const int64_t* nucleus_shell_num,
      const double* nucleus_range,
      const int32_t* nucleus_max_ang_mom,
      const int32_t* shell_ang_mom,
      const double* ao_factor,
      const double* shell_vgl,
      double* const ao_vgl );

#ifdef HAVE_HPC
    qmckl_exit_code qmckl_compute_ao_vgl_hpc_gaussian (
          const qmckl_context context,
          const int64_t ao_num,
          const int64_t shell_num,
          const int32_t* prim_num_per_nucleus,
          const int64_t point_num,
          const int64_t nucl_num,
          const double* coord,
          const double* nucl_coord,
          const int64_t* nucleus_index,
          const int64_t* nucleus_shell_num,
          const double* nucleus_range,
          const int32_t* nucleus_max_ang_mom,
          const int32_t* shell_ang_mom,
          const double* ao_factor,
          const qmckl_matrix expo_per_nucleus,
          const qmckl_tensor coef_per_nucleus,
          double* const ao_vgl );
#endif

Author: TREX CoE

Created: 2023-09-15 Fri 08:25

Validate