Atomic Orbitals

1. Polynomial part
- 1.1. Powers of \(x-X_i\)
- 1.2. Value, Gradient and Laplacian of a polynomial
2. Gaussian basis functions
3. TODO Slater basis functions

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.

1 Polynomial part

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

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

Requirements:

context is not 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,
               const double *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

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

  integer(c_int64_t), 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_context_get_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 /= 0) return

  test_qmckl_ao_power = -1

  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 = 0
  deallocate(X,P,LMAX)
end function test_qmckl_ao_power

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

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

Requirements:

context is not 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,
                        const int64_t *n,
                        const int32_t *L,
                        const int64_t ldl,
                        const double *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)
  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

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_context_get_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);
munit_assert_int(0, ==, test_qmckl_ao_polynomial_vgl(context));

2 Gaussian basis functions

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

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

`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_context_get_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

Atomic Orbitals

Table of Contents

1 Polynomial part

1.1 Powers of \(x-X_i\)

1.2 Value, Gradient and Laplacian of a polynomial

2 Gaussian basis functions

3 TODO Slater basis functions