TREX/qmckl

Fork 0

mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-06-01 02:45:43 +02:00

Anthony Scemama dcb392c0af Swap indices 1..5 with points in AOs/MOs

2022-02-11 16:07:25 +01:00

173 KiB

Raw Blame History

Atomic Orbitals

Introduction
Context
- Constant data
  - Initialization functions
    - C interface
    - Fortran interface
  - Access functions
    - C interface
    - Fortran interface
- Computed data
  - After initialization
  - Access functions
Radial part
Polynomial part
- General functions for Powers of $x-X_i$
- General functions for Value, Gradient and Laplacian of a polynomial
Combining radial and polynomial parts

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.

Context

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

For H_2 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 = 38

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

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 /TREX/qmckl/src/commit/dcb392c0afc4bcbef1c259b6fcd40a20b8d5227b/org/Context.

C interface

#+NAME:pre2

#+NAME:post2

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_ao_num (qmckl_context context,
                      const int64_t ao_num);

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_nucleus_shell_num (qmckl_context context,
                                 const int64_t* nucleus_shell_num,
                                 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_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);

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

Access functions

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);

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

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 primitives at current positions
`ao_vgl_date`	`uint64_t`	Last modification date of Value, gradients, Laplacian of the AOs at current positions

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.

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 /TREX/qmckl/src/commit/dcb392c0afc4bcbef1c259b6fcd40a20b8d5227b/org/Computation%20of%20primitives.

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 /TREX/qmckl/src/commit/dcb392c0afc4bcbef1c259b6fcd40a20b8d5227b/org/Computation%20of%20shells.

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 /TREX/qmckl/src/commit/dcb392c0afc4bcbef1c259b6fcd40a20b8d5227b/org/Combining%20radial%20and%20polynomial%20parts.

Radial part

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

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 = -dlog(1.d-15)

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

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
`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 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, 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)
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 = -dlog(1.d-12)

do inucl=1,nucl_num

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

     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

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}

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`	int64_t	in	Number of values
`X`	double[n]	in	Array containing the input values
`LMAX`	int32_t[n]	in	Array containing the maximum power for each value
`P`	double[n][ldp]	out	Array containing all the powers of `X`
`ldp`	int64_t	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

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 );

integer function qmckl_ao_polynomial_vgl_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,(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)
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

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

end function qmckl_ao_polynomial_vgl_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 );

integer function qmckl_ao_polynomial_transp_vgl_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)
integer           :: lmax_array(3)
real*8            :: 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 < (lmax+1)*(lmax+2)*(lmax+3)/6) then
 info = QMCKL_INVALID_ARG_9
 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(1,2:5) = 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

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

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

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

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

 l  (3,4) = 1
 VGL(4,1) = 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(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_f

Combining radial and polynomial parts

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

Unoptimized version

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_f

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

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

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

info = QMCKL_SUCCESS

! Don't compute polynomials when the radial part is zero.
cutoff = -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)
 k=1
 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 + z*z + z*z

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

    ! Compute polynomials 
    info = qmckl_ao_polynomial_vgl_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
       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

HPC version

integer function qmckl_compute_ao_vgl_hpc_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_transp_vgl_f

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

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

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

info = QMCKL_SUCCESS

! Don't compute polynomials when the radial part is zero.
cutoff = -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)
 k=1
 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 + z*z + z*z

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

    ! Compute polynomials 
    info = qmckl_ao_polynomial_transp_vgl_f(context, e_coord, n_coord, &
         nucleus_max_ang_mom(inucl), n_poly, powers, 3_8, &
         poly_vgl, int(ao_num,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
       l = shell_ang_mom(ishell)
       if (shell_vgl(ishell,1,ipoint) /= 0.d0) then
          do il = lstart(l), lstart(l+1)-1
             ! Value
             ao_vgl(k,1,ipoint) = &
                  poly_vgl(il,1) * shell_vgl(ishell,1,ipoint) * ao_factor(k)
             
             ! Grad_x
             ao_vgl(k,2,ipoint) = ( &
                  poly_vgl(il,2) * shell_vgl(ishell,1,ipoint) + &
                  poly_vgl(il,1) * shell_vgl(ishell,2,ipoint) &
                  ) * ao_factor(k)
             
             ! Grad_y
             ao_vgl(k,3,ipoint) = ( &
                  poly_vgl(il,3) * shell_vgl(ishell,1,ipoint) + &
                  poly_vgl(il,1) * shell_vgl(ishell,3,ipoint) &
                  ) * ao_factor(k)
             
             ! Grad_z
             ao_vgl(k,4,ipoint) = ( &
                  poly_vgl(il,4) * shell_vgl(ishell,1,ipoint) + &
                  poly_vgl(il,1) * shell_vgl(ishell,4,ipoint) &
                  ) * ao_factor(k)
             
             ! Lapl_z
             ao_vgl(k,5,ipoint) = ( &
                  poly_vgl(il,5) * shell_vgl(ishell,1,ipoint) + &
                  poly_vgl(il,1) * shell_vgl(ishell,5,ipoint) + &
                  2.d0 * ( &
                  poly_vgl(il,2) * shell_vgl(ishell,2,ipoint) + &
                  poly_vgl(il,3) * shell_vgl(ishell,3,ipoint) + &
                  poly_vgl(il,4) * shell_vgl(ishell,4,ipoint) ) &
                  ) * ao_factor(k)
             k = k+1
          end do
       else
          do il = lstart(l), lstart(l+1)-1
             ao_vgl(k,1,ipoint) = 0.d0
             ao_vgl(k,2,ipoint) = 0.d0
             ao_vgl(k,3,ipoint) = 0.d0
             ao_vgl(k,4,ipoint) = 0.d0
             ao_vgl(k,5,ipoint) = 0.d0
             k = k+1
          end do
       end if
    end do
 end do
end do

deallocate(poly_vgl, powers)
end function qmckl_compute_ao_vgl_hpc_f

Interfaces

qmckl_exit_code qmckl_compute_ao_vgl (
      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 );

173 KiB Raw Blame History

Atomic Orbitals

Introduction

Context

Constant data

Initialization functions

C interface

Fortran interface

Access functions

C interface

Fortran interface

Computed data

After initialization

Access functions

Radial part

General functions for Gaussian basis functions

Computation of primitives

Computation of shells

Polynomial part

General functions for Powers of $x-X_i$

General functions for Value, Gradient and Laplacian of a polynomial

Combining radial and polynomial parts

Unoptimized version

HPC version

Interfaces

173 KiB

Raw Blame History