qmckl/org/qmckl_ao.org

Atomic Orbitals

• Context
  • Data structure
  • Access functions
  • Initialization functions
  • Fortran interfaces
• Polynomial part
  • Powers of $x-X_i$
    • Requirements
    • C Header
    • Source
    • C interface
    • Fortran interface
    • Test
  • Value, Gradient and Laplacian of a polynomial
    • Requirements
    • C Header
    • Source
    • C interface
    • Fortran interface
    • Test
• Radial part
  • Gaussian basis functions
  • Slater basis functions
  • Radial functions on a grid
• Combining radial and polynomial parts

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 \emph{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} \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 to unity. As this normalization requires the ability to compute overlap integrals, it should be written in the file to ensure that the file is self-contained and does not require the client program to have the ability to compute such integrals.

Atomic orbitals (AOs) are defined as

\[ \chi_i (\mathbf{r}) = 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.

In this section we describe the kernels used to compute the values, gradients and Laplacian of the atomic basis functions.

Context

The following arrays are stored in the context:

type		Gaussian ('G') or Slater ('S')
shell_num		Number of shells
prim_num		Total number of primitives
shell_center	[shell_num]	Id of the nucleus on which each shell is centered
shell_ang_mom	[shell_num]	Angular momentum of each shell
shell_prim_num	[shell_num]	Number of primitives in each shell
shell_prim_index	[shell_num]	Address of the first primitive of each shell in the EXPONENT array
shell_factor	[shell_num]	Normalization factor for each shell
exponent	[prim_num]	Array of exponents
coefficient	[prim_num]	Array of coefficients

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
shell_center = [1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2]
shell_ang_mom = ['S', 'S', 'S', 'P', 'P', 'D', 'S', 'S', 'S', 'P', 'P', 'D']
shell_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 = [1, 6, 7, 8, 9, 10, 11, 16, 17, 18, 19, 20]
exponent = [ 33.87, 5.095, 1.159, 0.3258, 0.1027, 0.3258, 0.1027,
          1.407, 0.388, 1.057, 33.87, 5.095, 1.159, 0.3258, 0.1027,
          0.3258, 0.1027, 1.407, 0.388, 1.057]
coefficient = [ 0.006068, 0.045308, 0.202822, 0.503903, 0.383421,
             1.0, 1.0, 1.0, 1.0, 1.0, 0.006068, 0.045308, 0.202822,
             0.503903, 0.383421, 1.0, 1.0, 1.0, 1.0, 1.0]

Data structure

typedef struct qmckl_ao_basis_struct {
int32_t   uninitialized;
int64_t   shell_num;
int64_t   prim_num;
int64_t * shell_center;
char    * shell_ang_mom;
int64_t * shell_prim_num;
int64_t * shell_prim_index;
double  * shell_factor;
double  * exponent    ;
double  * coefficient ;
bool      provided;
char      type;
} qmckl_ao_basis_struct;

The uninitialized integer contains one bit set to one for each initialization function which has not bee called. It becomes equal to zero after all initialization functions have been called. The struct is then initialized and provided == true.

Access functions

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

bool      qmckl_ao_basis_provided           (const qmckl_context context);

#+NAME:post

Initialization functions

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

qmckl_exit_code  qmckl_set_ao_basis_type             (qmckl_context context, const char      t);
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_shell_prim_index (qmckl_context context, const int64_t * shell_prim_index);
qmckl_exit_code  qmckl_set_ao_basis_shell_center     (qmckl_context context, const int64_t * shell_center);
qmckl_exit_code  qmckl_set_ao_basis_shell_ang_mom    (qmckl_context context, const char    * shell_ang_mom);
qmckl_exit_code  qmckl_set_ao_basis_shell_prim_num   (qmckl_context context, const int64_t * shell_prim_num);
qmckl_exit_code  qmckl_set_ao_basis_shell_factor     (qmckl_context context, const double  * shell_factor);
qmckl_exit_code  qmckl_set_ao_basis_exponent         (qmckl_context context, const double  * exponent);
qmckl_exit_code  qmckl_set_ao_basis_coefficient      (qmckl_context context, const double  * coefficient);

#+NAME:pre2

#+NAME:post2

TODO Fortran interfaces

Polynomial part

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 \]

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

C Header

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

Source

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

C interface

Fortran interface

Test

integer(c_int32_t) function test_qmckl_ao_power(context) bind(C)
use qmckl
implicit none

integer(qmckl_context), intent(in), value :: context

integer*8                     :: n, LDP
integer, allocatable          :: LMAX(:)
double precision, allocatable :: X(:), P(:,:)
integer*8                     :: i,j
double precision              :: epsilon

epsilon = qmckl_get_numprec_epsilon(context)

n = 100;
LDP = 10;

allocate(X(n), P(LDP,n), LMAX(n))

do j=1,n
 X(j) = -5.d0 + 0.1d0 * dble(j)
 LMAX(j) = 1 + int(mod(j, 5),4)
end do

test_qmckl_ao_power = qmckl_ao_power(context, n, X, LMAX, P, LDP)
if (test_qmckl_ao_power /= QMCKL_SUCCESS) return

test_qmckl_ao_power = QMCKL_FAILURE

do j=1,n
 do i=1,LMAX(j)
    if ( X(j)**i == 0.d0 ) then
       if ( P(i,j) /= 0.d0) return
    else
       if ( dabs(1.d0 - P(i,j) / (X(j)**i)) > epsilon ) return
    end if
 end do
end do

test_qmckl_ao_power = QMCKL_SUCCESS
deallocate(X,P,LMAX)
end function test_qmckl_ao_power

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.

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

Requirements

• context is not 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"

C Header

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

Source

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)
integer, external :: qmckl_ao_power_f
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

C interface

Fortran interface

Test

integer(c_int32_t) function test_qmckl_ao_polynomial_vgl(context) bind(C)
use qmckl
implicit none

integer(c_int64_t), intent(in), value :: context

integer                       :: lmax, d, i
integer, allocatable          :: L(:,:)
integer*8                     :: n, ldl, ldv, j
double precision              :: X(3), R(3), Y(3)
double precision, allocatable :: VGL(:,:)
double precision              :: w
double precision              :: epsilon

epsilon = qmckl_get_numprec_epsilon(context)

X = (/ 1.1 , 2.2 ,  3.3 /)
R = (/ 0.1 , 1.2 , -2.3 /)
Y(:) = X(:) - R(:)

lmax = 4;
ldl = 3;
ldv = 100;

d = (lmax+1)*(lmax+2)*(lmax+3)/6

allocate (L(ldl,d), VGL(ldv,d))

test_qmckl_ao_polynomial_vgl = &
   qmckl_ao_polynomial_vgl(context, X, R, lmax, n, L, ldl, VGL, ldv)

if (test_qmckl_ao_polynomial_vgl /= QMCKL_SUCCESS) return
if (n /= d) return

do j=1,n
 test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE
 do i=1,3
    if (L(i,j) < 0) return
 end do
 test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE
 if (dabs(1.d0 - VGL(1,j) / (&
      Y(1)**L(1,j) * Y(2)**L(2,j) * Y(3)**L(3,j)  &
      )) > epsilon ) return

 test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE
 if (L(1,j) < 1) then
    if (VGL(2,j) /= 0.d0) return
 else
    if (dabs(1.d0 - VGL(2,j) / (&
         L(1,j) * Y(1)**(L(1,j)-1) * Y(2)**L(2,j) * Y(3)**L(3,j) &
         )) > epsilon ) return
 end if

 test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE
 if (L(2,j) < 1) then
    if (VGL(3,j) /= 0.d0) return
 else
    if (dabs(1.d0 - VGL(3,j) / (&
         L(2,j) * Y(1)**L(1,j) * Y(2)**(L(2,j)-1) * Y(3)**L(3,j) &
         )) > epsilon ) return
 end if

 test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE
 if (L(3,j) < 1) then
    if (VGL(4,j) /= 0.d0) return
 else
    if (dabs(1.d0 - VGL(4,j) / (&
         L(3,j) * Y(1)**L(1,j) * Y(2)**L(2,j) * Y(3)**(L(3,j)-1) &
         )) > epsilon ) return
 end if

 test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE
 w = 0.d0
 if (L(1,j) > 1) then
    w = w + L(1,j) * (L(1,j)-1) * Y(1)**(L(1,j)-2) * Y(2)**L(2,j) * Y(3)**L(3,j)
 end if
 if (L(2,j) > 1) then
    w = w + L(2,j) * (L(2,j)-1) * Y(1)**L(1,j) * Y(2)**(L(2,j)-2) * Y(3)**L(3,j)
 end if
 if (L(3,j) > 1) then
    w = w + L(3,j) * (L(3,j)-1) * Y(1)**L(1,j) * Y(2)**L(2,j) * Y(3)**(L(3,j)-2)
 end if
 if (dabs(1.d0 - VGL(5,j) / w) > epsilon ) return
end do

test_qmckl_ao_polynomial_vgl = QMCKL_SUCCESS

deallocate(L,VGL)
end function test_qmckl_ao_polynomial_vgl

int  test_qmckl_ao_polynomial_vgl(qmckl_context context);
assert(0 == test_qmckl_ao_polynomial_vgl(context));

Radial part

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 \]

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

Requirements

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

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

integer(c_int32_t) function test_qmckl_ao_gaussian_vgl(context) bind(C)
use qmckl
implicit none

integer(c_int64_t), intent(in), value :: context

integer*8                     :: n, ldv, j, i
double precision              :: X(3), R(3), Y(3), r2
double precision, allocatable :: VGL(:,:), A(:)
double precision              :: epsilon

epsilon = qmckl_get_numprec_epsilon(context)

X = (/ 1.1 , 2.2 ,  3.3 /)
R = (/ 0.1 , 1.2 , -2.3 /)
Y(:) = X(:) - R(:)
r2 = Y(1)**2 + Y(2)**2 + Y(3)**2

n = 10;
ldv = 100;

allocate (A(n), VGL(ldv,5))
do i=1,n
  A(i) = 0.0013 * dble(ishft(1,i))
end do


test_qmckl_ao_gaussian_vgl = &
    qmckl_ao_gaussian_vgl(context, X, R, n, A, VGL, ldv)
if (test_qmckl_ao_gaussian_vgl /= 0) return

test_qmckl_ao_gaussian_vgl = -1

do i=1,n
  test_qmckl_ao_gaussian_vgl = -11
  if (dabs(1.d0 - VGL(i,1) / (&
       dexp(-A(i) * r2) &
       )) > epsilon ) return

  test_qmckl_ao_gaussian_vgl = -12
  if (dabs(1.d0 - VGL(i,2) / (&
       -2.d0 * A(i) * Y(1) * dexp(-A(i) * r2) &
       )) > epsilon ) return

  test_qmckl_ao_gaussian_vgl = -13
  if (dabs(1.d0 - VGL(i,3) / (&
       -2.d0 * A(i) * Y(2) * dexp(-A(i) * r2) &
       )) > epsilon ) return

  test_qmckl_ao_gaussian_vgl = -14
  if (dabs(1.d0 - VGL(i,4) / (&
       -2.d0 * A(i) * Y(3) * dexp(-A(i) * r2) &
       )) > epsilon ) return

  test_qmckl_ao_gaussian_vgl = -15
  if (dabs(1.d0 - VGL(i,5) / (&
       A(i) * (4.d0*r2*A(i) - 6.d0) * dexp(-A(i) * r2) &
       )) > epsilon ) return
end do

test_qmckl_ao_gaussian_vgl = 0

deallocate(VGL)
end function test_qmckl_ao_gaussian_vgl

TODO Slater basis functions

TODO Radial functions on a grid

Combining radial and polynomial parts