Atomic Orbitals

Table of Contents

This files contains all the routines for the computation of the values, gradients and Laplacian of the atomic basis functions.

4 files are produced:

1 Polynomials

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

\begin{eqnarray*} \frac{\partial }{\partial x} P_l\left(\mathbf{r},\mathbf{R}_i \right) & = & a (x-X_i)^{a-1} (y-Y_i)^b (z-Z_i)^c \\ \frac{\partial }{\partial y} P_l\left(\mathbf{r},\mathbf{R}_i \right) & = & b (x-X_i)^a (y-Y_i)^{b-1} (z-Z_i)^c \\ \frac{\partial }{\partial z} P_l\left(\mathbf{r},\mathbf{R}_i \right) & = & c (x-X_i)^a (y-Y_i)^b (z-Z_i)^{c-1} \\ \end{eqnarray*} \begin{eqnarray*} \left( \frac{\partial }{\partial x^2} + \frac{\partial }{\partial y^2} + \frac{\partial }{\partial z^2} \right) P_l \left(\mathbf{r},\mathbf{R}_i \right) & = & a(a-1) (x-X_i)^{a-2} (y-Y_i)^b (z-Z_i)^c + \\ && b(b-1) (x-X_i)^a (y-Y_i)^{b-1} (z-Z_i)^c + \\ && c(c-1) (x-X_i)^a (y-Y_i)^b (z-Z_i)^{c-1} \end{eqnarray*}

1.1 qmckl_ao_powers

Computes all the powers of the n input data up to the given maximum value given in input for each of the \(n\) points:

\[ P_{ij} = X_j^i \]

1.1.1 Arguments

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

1.1.2 Requirements

  • context is not 0
  • n > 0
  • X is allocated with at least \(n \times 8\) bytes
  • LMAX is allocated with at least \(n \times 4\) bytes
  • P is allocated with at least \(n \times \max_i \text{LMAX}_i \times 8\) bytes
  • LDP >= \(\max_i\) LMAX[i]

1.1.3 Header

qmckl_exit_code qmckl_ao_powers(const qmckl_context context,
		const int64_t n, 
		const double *X, const int32_t *LMAX,
		const double *P, const int64_t LDP);

1.1.4 Source

integer function qmckl_ao_powers_f(context, n, X, LMAX, P, ldp) result(info)
  implicit none
  integer*8 , intent(in)  :: context
  integer*8 , intent(in)  :: n
  real*8    , intent(in)  :: X(n)
  integer   , intent(in)  :: LMAX(n)
  real*8    , intent(out) :: P(ldp,n)
  integer*8 , intent(in)  :: ldp

  integer*8  :: i,j

  info = 0

  if (context == 0_8) then
     info = -1
     return
  endif

  if (LDP < MAXVAL(LMAX)) then
     info = -2
     return
  endif

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

end function qmckl_ao_powers_f

1.2 qmckl_ao_polynomial_vgl

Computes the values, gradients and Laplacians at a given point of all polynomials with an angular momentum up to lmax.

1.2.1 Arguments

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

1.2.2 Requirements

  • context is not 0
  • n > 0
  • lmax >= 0
  • ldl >= 3
  • ldv >= (=lmax=+1)(=lmax=+2)(=lmax=+3)/6
  • X is allocated with at least \(3 \times 8\) bytes
  • R is allocated with at least \(3 \times 8\) bytes
  • L is allocated with at least \(3 \times n \times 4\) bytes
  • VGL is allocated with at least \(n \times 5 \times 8\) bytes
  • On output, n should be equal to (=lmax=+1)(=lmax=+2)(=lmax=+3)/6

1.2.3 Header

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

1.2.4 Source

integer function qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL, ldv) result(info)
  implicit none
  integer*8 , intent(in)  :: context
  real*8    , intent(in)  :: X(3), R(3)
  integer   , intent(in)  :: lmax
  integer*8 , intent(out) :: n
  integer   , intent(out) :: L(ldl,(lmax+1)*(lmax+2)*(lmax+3)/6)
  integer*8 , intent(in)  :: ldl
  real*8    , intent(out) :: VGL(ldv,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)
  integer, external :: qmckl_ao_powers_f
  double precision  :: xy, yz, xz
  double precision  :: da, db, dc, dd

  info = 0

  if (context == 0_8) then
     info = -1
     return
  endif

  if (ldl < 3) then
     info = -2
     return
  endif

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

  if (lmax <= 0) then
     info = -4
     return
  endif


  do i=1,3
     Y(i) = X(i) - R(i)
  end do
  pows(-2:-1,1:3) = 0.d0
  pows(0,1:3) = 1.d0
  lmax_array(1:3) = lmax
  info = qmckl_ao_powers_f(context, 1_8, Y(1), (/lmax/), pows(1,1), size(pows,1,kind=8)) 
  if (info /= 0) return
  info = qmckl_ao_powers_f(context, 1_8, Y(2), (/lmax/), pows(1,2), size(pows,1,kind=8)) 
  if (info /= 0) return
  info = qmckl_ao_powers_f(context, 1_8, Y(3), (/lmax/), pows(1,3), size(pows,1,kind=8)) 
  if (info /= 0) return


  vgl(1,1) = 1.d0
  vgl(1,2:5) = 0.d0
  l(1:3,1) = 0
  n=1
  dd = 1.d0
  do d=1,lmax
     da = 0.d0
     do a=0,d
	db = 0.d0
	do b=0,d-a
	   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

  if (n /= (lmax+1)*(lmax+2)*(lmax+3)/6) then
    info = -5
    return
  endif

  info = 0

end function qmckl_ao_polynomial_vgl_f

2 Gaussian basis functions

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

2.1.1 Arguments

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

2.1.2 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

2.1.3 Header

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

2.1.4 Source

integer function qmckl_ao_gaussians_vgl_f(context, X, R, n, A, VGL, ldv) result(info)
  implicit none
  integer*8 , intent(in)  :: context
  real*8    , intent(in)  :: X(3), R(3)
  integer*8 , intent(in)  :: n
  real*8    , intent(in)  :: A(n)
  real*8    , intent(out) :: VGL(ldv,5)
  integer*8 , intent(in)  :: ldv

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

  info = 0

  if (context == 0_8) then
     info = -1
     return
  endif

  if (n <= 0) then
     info = -2
     return
  endif

  if (ldv < n) then
     info = -3
     return
  endif


  do i=1,3
     Y(i) = X(i) - R(i)
  end do
  r2 = Y(1)*Y(1) + Y(2)*Y(2) + Y(3)*Y(3)

  do i=1,n
     VGL(i,1) = dexp(-A(i) * r2)
  end do

  do i=1,n
     VGL(i,5) = A(i) * VGL(i,1)
  end do

  t = -2.d0 * ( X(1) - R(1) )
  u = -2.d0 * ( X(2) - R(2) )
  v = -2.d0 * ( X(3) - R(3) )

  do i=1,n
     VGL(i,2) = t * VGL(i,5)
     VGL(i,3) = u * VGL(i,5)
     VGL(i,4) = v * VGL(i,5)
  end do

  t = 4.d0 * r2
  do i=1,n
     VGL(i,5) = (t * A(i) - 6.d0) *  VGL(i,5)
  end do

end function qmckl_ao_gaussians_vgl_f

3 TODO Slater basis functions

Created: 2020-11-04 Wed 23:47

Emacs 25.2.2 (Org mode 8.2.10)

Validate