mirror of
https://github.com/TREX-CoE/qmckl.git
synced 2024-11-19 20:42:50 +01:00
128 KiB
128 KiB
Atomic Orbitals
- Context
- Radial part
- Polynomial part
- 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}\, 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
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 | |
nucleus_index |
[nucl_num] |
Index of the first shell of each nucleus |
nucleus_shell_num |
[nucl_num] |
Number of shells per nucleus |
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 |
prim_factor |
[prim_num] |
Normalization factors of the primtives |
ao_num |
Number of AOs | |
ao_cartesian |
If true, use polynomials. Otherwise, use spherical harmonics | |
ao_factor |
[ao_num] |
Normalization factor of the AO |
ao_shell |
[ao_num] |
For each AO, specify to which shell it belongs |
Computed data:
coefficient_normalized |
[prim_num] |
Normalized primitive coefficients |
|---|---|---|
nucleus_prim_index |
[nucl_num] |
Index of the first primitive for each nucleus |
nucleus_max_ang_mom |
[nucl_num] |
Maximum angular momentum for each nucleus |
nucleus_range |
[nucl_num] |
Distance beyond which all the AOs are zero |
primitive_vgl |
[5][walk_num][elec_num][prim_num] |
Value, gradients, Laplacian of the primitives at electron positions |
primitive_vgl_date |
uint64_t |
Late modification date of Value, gradients, Laplacian of the primitives at electron positions |
shell_vgl |
[5][walk_num][elec_num][shell_num] |
Value, gradients, Laplacian of the primitives at electron positions |
shell_vgl_date |
uint64_t |
Late modification date of Value, gradients, Laplacian of the AOs at electron positions |
ao_vgl |
[5][walk_num][elec_num][ao_num] |
Value, gradients, Laplacian of the primitives at electron positions |
ao_vgl_date |
uint64_t |
Late modification date of Value, gradients, Laplacian of the AOs at electron positions |
nucl_shell_index |
[nucl_num] |
Index of the first shell for each nucleus |
exponent_sorted |
[prim_num] |
Array of exponents for sorted primitives |
coeff_norm_sorted |
[prim_num] |
Array of normalized coefficients for sorted primitives |
prim_factor_sorted |
[prim_num] |
Normalization factors of the sorted primtives |
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 ]
Data structure
typedef struct qmckl_ao_basis_struct {
int64_t shell_num;
int64_t prim_num;
int64_t ao_num;
int64_t * nucleus_index;
int64_t * nucleus_shell_num;
int32_t * shell_ang_mom;
int64_t * shell_prim_num;
int64_t * shell_prim_index;
double * shell_factor;
double * exponent;
double * coefficient;
double * prim_factor;
double * ao_factor;
int64_t * nucleus_prim_index;
double * coefficient_normalized;
int32_t * nucleus_max_ang_mom;
double * nucleus_range;
double * primitive_vgl;
int64_t primitive_vgl_date;
double * shell_vgl;
int64_t shell_vgl_date;
double * ao_vgl;
int64_t ao_vgl_date;
int32_t uninitialized;
bool provided;
bool ao_cartesian;
char type;
} qmckl_ao_basis_struct;
The uninitialized integer contains one bit set to one for each
initialization function which has not been called. It becomes equal
to zero after all initialization functions have been called. The
struct is then initialized and provided == true.
Some values are initialized by default, and are not concerned by
this mechanism.
qmckl_exit_code qmckl_init_ao_basis(qmckl_context context);qmckl_exit_code qmckl_init_ao_basis(qmckl_context context) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return false;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
ctx->ao_basis.uninitialized = (1 << 14) - 1;
/* Default values */
ctx->ao_basis.ao_cartesian = true;
return QMCKL_SUCCESS;
}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.
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_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);
qmckl_exit_code qmckl_set_ao_basis_nucleus_shell_num(qmckl_context context, const int64_t * nucleus_shell_num);
qmckl_exit_code qmckl_set_ao_basis_shell_ang_mom (qmckl_context context, const int32_t * 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_prim_index (qmckl_context context, const int64_t * shell_prim_index);
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);
qmckl_exit_code qmckl_set_ao_basis_prim_factor (qmckl_context context, const double * prim_factor);
qmckl_exit_code qmckl_set_ao_basis_ao_factor (qmckl_context context, const double * ao_factor);
qmckl_exit_code qmckl_set_ao_basis_cartesian (qmckl_context context, const bool cartesian);#+NAME:pre2
#+NAME:post2
When the basis set is completely entered, other data structures are 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.
Fortran interfaces
Radial part
TODO Helper functions to accelerate calculations
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 \]
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
contextis not 0n> 0ldv>= 5A(i)> 0 for alliXis allocated with at least $3 \times 8$ bytesRis allocated with at least $3 \times 8$ bytesAis allocated with at least $n \times 8$ bytesVGLis 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_finteger(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_vglTODO General functions for Slater basis functions
TODO General functions for Radial functions on a grid
Computation of primitives
Get
qmckl_exit_code qmckl_get_ao_basis_primitive_vgl(qmckl_context context, double* const primitive_vgl);Provide
Compute
| qmckl_context | context | in | Global state |
| int64_t | prim_num | in | Number of primitives |
| int64_t | elec_num | in | Number of electrons |
| int64_t | nucl_num | in | Number of nuclei |
| int64_t | walk_num | in | Number of walkers |
| int64_t | nucleus_prim_index[nucl_num] | in | Index of the 1st primitive of each nucleus |
| double | elec_coord[walk_num][3][elec_num] | in | Electron coordinates |
| double | nucl_coord[3][elec_num] | in | Nuclear coordinates |
| double | expo[prim_num] | in | Exponents of the primitives |
| double | primitive_vgl[5][walk_num][elec_num][prim_num] | out | Value, gradients and Laplacian of the primitives |
integer function qmckl_compute_ao_basis_primitive_gaussian_vgl_f(context, &
prim_num, elec_num, nucl_num, walk_num, &
nucleus_prim_index, elec_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) :: elec_num
integer*8 , intent(in) :: walk_num
integer*8 , intent(in) :: nucleus_prim_index(nucl_num+1)
double precision , intent(in) :: elec_coord(elec_num,3,walk_num)
double precision , intent(in) :: nucl_coord(nucl_num,3)
double precision , intent(in) :: expo(prim_num)
double precision , intent(out) :: primitive_vgl(prim_num,elec_num,walk_num,5)
integer*8 :: inucl, iprim, iwalk, ielec
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 iwalk = 1, walk_num
do ielec = 1, elec_num
x = elec_coord(ielec,1,iwalk) - nucl_coord(inucl,1)
y = elec_coord(ielec,2,iwalk) - nucl_coord(inucl,2)
z = elec_coord(ielec,3,iwalk) - 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, ielec, iwalk, 1) = v
primitive_vgl(iprim, ielec, iwalk, 2) = two_a * x
primitive_vgl(iprim, ielec, iwalk, 3) = two_a * y
primitive_vgl(iprim, ielec, iwalk, 4) = two_a * z
primitive_vgl(iprim, ielec, iwalk, 5) = two_a * (3.d0 - 2.d0*ar2)
end do
end do
end do
end do
end function qmckl_compute_ao_basis_primitive_gaussian_vgl_fTest
Ideas for improvement
// m : walkers
// j : electrons
// l : primitives
k=0;
for (m=0 ; m<walk_num ; ++m) {
for (j=0 ; j<elec_num ; ++j) {
for (i=0 ; i<nucl_num ; ++i) {
r2 = nucl_elec_dist[i][j];
if (r2 < nucl_radius2[i]) {
for (l=0 ; l<prim_num ; ++l) {
tmp[k].i = i;
tmp[k].j = j;
tmp[k].m = m;
tmp[k].ar2 = -expo[l] *r2;
++k;
}
}
}
}
}
// sort(tmp) in increasing ar2;
// Identify first ar2 above numerical accuracy threshold
// Compute vectorized exponentials on significant values
Computation of shells
Get
qmckl_exit_code qmckl_get_ao_basis_shell_vgl(qmckl_context context, double* const shell_vgl);Provide
Compute
qmckl_context |
context |
in | Global state |
int64_t |
prim_num |
in | Number of primitives |
int64_t |
shell_num |
in | Number of shells |
int64_t |
elec_num |
in | Number of electrons |
int64_t |
nucl_num |
in | Number of nuclei |
int64_t |
walk_num |
in | Number of walkers |
int64_t |
nucleus_shell_num[nucl_num] |
in | Number of shells for each nucleus |
int64_t |
nucleus_index[nucl_num] |
in | Index of the 1st shell of each nucleus |
int64_t |
shell_prim_index[shell_num] |
in | Index of the 1st primitive of each shell |
int64_t |
shell_prim_num[shell_num] |
in | Number of primitives per shell |
double |
elec_coord[walk_num][3][elec_num] |
in | Electron coordinates |
double |
nucl_coord[3][elec_num] |
in | Nuclear coordinates |
double |
expo[prim_num] |
in | Exponents of the primitives |
double |
coef_normalized[prim_num] |
in | Coefficients of the primitives |
double |
shell_vgl[5][walk_num][elec_num][shell_num] |
out | Value, gradients and Laplacian of the shells |
integer function qmckl_compute_ao_basis_shell_gaussian_vgl_f(context, &
prim_num, shell_num, elec_num, nucl_num, walk_num, &
nucleus_shell_num, nucleus_index, shell_prim_index, shell_prim_num, &
elec_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) :: elec_num
integer*8 , intent(in) :: walk_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) :: elec_coord(elec_num,3,walk_num)
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,elec_num,walk_num,5)
integer*8 :: inucl, iprim, iwalk, ielec, ishell
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-15)
do inucl=1,nucl_num
do iwalk = 1, walk_num
do ielec = 1, elec_num
x = elec_coord(ielec,1,iwalk) - nucl_coord(inucl,1)
y = elec_coord(ielec,2,iwalk) - nucl_coord(inucl,2)
z = elec_coord(ielec,3,iwalk) - nucl_coord(inucl,3)
r2 = x*x + y*y + z*z
do ishell=nucleus_index(inucl)+1, nucleus_index(inucl)+nucleus_shell_num(inucl)
! C is zero-based, so shift bounds by one
shell_vgl(ishell, ielec, iwalk, 1) = 0.d0
shell_vgl(ishell, ielec, iwalk, 2) = 0.d0
shell_vgl(ishell, ielec, iwalk, 3) = 0.d0
shell_vgl(ishell, ielec, iwalk, 4) = 0.d0
shell_vgl(ishell, ielec, iwalk, 5) = 0.d0
do iprim = shell_prim_index(ishell)+1, shell_prim_index(ishell)+shell_prim_num(ishell)
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, ielec, iwalk, 1) = &
shell_vgl(ishell, ielec, iwalk, 1) + v
shell_vgl(ishell, ielec, iwalk, 2) = &
shell_vgl(ishell, ielec, iwalk, 2) + two_a * x
shell_vgl(ishell, ielec, iwalk, 3) = &
shell_vgl(ishell, ielec, iwalk, 3) + two_a * y
shell_vgl(ishell, ielec, iwalk, 4) = &
shell_vgl(ishell, ielec, iwalk, 4) + two_a * z
shell_vgl(ishell, ielec, iwalk, 5) = &
shell_vgl(ishell, ielec, iwalk, 5) + two_a * (3.d0 - 2.d0*ar2)
end do
end do
end do
end do
end do
end function qmckl_compute_ao_basis_shell_gaussian_vgl_fqmckl_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 elec_num,
const int64_t nucl_num,
const int64_t walk_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* elec_coord,
const double* nucl_coord,
const double* expo,
const double* coef_normalized,
double* const shell_vgl );Test
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 \]
| 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
contextis notQMCKL_NULL_CONTEXTn> 0Xis allocated with at least $n \times 8$ bytesLMAXis allocated with at least $n \times 4$ bytesPis allocated with at least $n \times \max_i \text{LMAX}_i \times 8$ bytesLDP>= $\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_fC 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_powerGeneral 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.
| 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
contextis notQMCKL_NULL_CONTEXTn> 0lmax>= 0ldl>= 3ldv>= 5Xis allocated with at least $3 \times 8$ bytesRis allocated with at least $3 \times 8$ bytesn>=(lmax+1)(lmax+2)(lmax+3)/6Lis allocated with at least $3 \times n \times 4$ bytesVGLis allocated with at least $5 \times n \times 8$ bytes- On output,
nshould 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"
- Increasing values of
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)
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_fC 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_vglint test_qmckl_ao_polynomial_vgl(qmckl_context context);
assert(0 == test_qmckl_ao_polynomial_vgl(context));Combining radial and polynomial parts
Get
qmckl_exit_code qmckl_get_ao_vgl(qmckl_context context, double* const ao_vgl);Provide
Compute
qmckl_context |
context |
in | Global state |
int64_t |
ao_num |
in | Number of AOs |
int64_t |
shell_num |
in | Number of shells |
int64_t |
elec_num |
in | Number of electrons |
int64_t |
nucl_num |
in | Number of nuclei |
int64_t |
walk_num |
in | Number of walkers |
double |
elec_coord[walk_num][3][elec_num] |
in | Electron coordinates |
double |
nucl_coord[3][nucl_num] |
in | Nuclear coordinates |
int64_t |
nucleus_index[nucl_num] |
in | Index of the 1st shell of each nucleus |
int64_t |
nucleus_shell_num[nucl_num] |
in | Number of shells per nucleus |
double |
nucleus_range[nucl_num] |
in | Range beyond which all is zero |
int32_t |
nucleus_max_ang_mom[nucl_num] |
in | Maximum angular momentum per nucleus |
int32_t |
shell_ang_mom[shell_num] |
in | Angular momentum of each shell |
double |
ao_factor[ao_num] |
in | Normalization factor of the AOs |
double |
shell_vgl[5][walk_num][elec_num][shell_num] |
in | Value, gradients and Laplacian of the shells |
double |
ao_vgl[5][walk_num][elec_num][ao_num] |
out | Value, gradients and Laplacian of the AOs |
integer function qmckl_compute_ao_vgl_f(context, &
ao_num, shell_num, elec_num, nucl_num, walk_num, &
elec_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) :: elec_num
integer*8 , intent(in) :: nucl_num
integer*8 , intent(in) :: walk_num
double precision , intent(in) :: elec_coord(elec_num,3,walk_num)
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,elec_num,walk_num,5)
double precision , intent(out) :: ao_vgl(ao_num,elec_num,walk_num,5)
double precision :: e_coord(3), n_coord(3)
integer*8 :: n_poly
integer :: l, il, k
integer*8 :: ielec, inucl, ishell, iwalk
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.
! TODO : Use numerical precision here
cutoff = -dlog(1.d-15)
do iwalk = 1,walk_num
do ielec = 1, elec_num
e_coord(1) = elec_coord(ielec,1,iwalk)
e_coord(2) = elec_coord(ielec,2,iwalk)
e_coord(3) = elec_coord(ielec,3,iwalk)
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 electron 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
do ishell = nucleus_index(inucl)+1, nucleus_index(inucl)+nucleus_shell_num(inucl)
l = shell_ang_mom(ishell)
do il = lstart(l), lstart(l+1)-1
! Value
ao_vgl(k,ielec,iwalk,1) = &
poly_vgl(1,il) * shell_vgl(ishell,ielec,iwalk,1) * ao_factor(k)
! Grad_x
ao_vgl(k,ielec,iwalk,2) = ( &
poly_vgl(2,il) * shell_vgl(ishell,ielec,iwalk,1) + &
poly_vgl(1,il) * shell_vgl(ishell,ielec,iwalk,2) &
) * ao_factor(k)
! Grad_y
ao_vgl(k,ielec,iwalk,3) = ( &
poly_vgl(3,il) * shell_vgl(ishell,ielec,iwalk,1) + &
poly_vgl(1,il) * shell_vgl(ishell,ielec,iwalk,3) &
) * ao_factor(k)
! Grad_z
ao_vgl(k,ielec,iwalk,4) = ( &
poly_vgl(4,il) * shell_vgl(ishell,ielec,iwalk,1) + &
poly_vgl(1,il) * shell_vgl(ishell,ielec,iwalk,4) &
) * ao_factor(k)
! Lapl_z
ao_vgl(k,ielec,iwalk,5) = ( &
poly_vgl(5,il) * shell_vgl(ishell,ielec,iwalk,1) + &
poly_vgl(1,il) * shell_vgl(ishell,ielec,iwalk,5) + &
2.d0 * ( &
poly_vgl(2,il) * shell_vgl(ishell,ielec,iwalk,2) + &
poly_vgl(3,il) * shell_vgl(ishell,ielec,iwalk,3) + &
poly_vgl(4,il) * shell_vgl(ishell,ielec,iwalk,4) ) &
) * ao_factor(k)
k = k+1
end do
end do
end do
end do
end do
deallocate(poly_vgl, powers)
end function qmckl_compute_ao_vgl_fqmckl_exit_code qmckl_compute_ao_vgl (
const qmckl_context context,
const int64_t ao_num,
const int64_t shell_num,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t walk_num,
const double* elec_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 );