1
0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-12-23 04:44:03 +01:00

Merge branch 'master' into ao_mo_vj

This commit is contained in:
vijay 2021-09-29 00:33:07 +02:00 committed by GitHub
commit 5557187b0d
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
2 changed files with 59322 additions and 389 deletions

View File

@ -37,6 +37,9 @@ 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, context, and then we present the kernels used to compute the values,
gradients and Laplacian of the atomic basis functions. gradients and Laplacian of the atomic basis functions.
* TODO [0/1] Missing features :noexport:
- [ ] Error messages to tell what is missing when initializing
* Headers :noexport: * Headers :noexport:
#+begin_src elisp :noexport :results none #+begin_src elisp :noexport :results none
(org-babel-lob-ingest "../tools/lib.org") (org-babel-lob-ingest "../tools/lib.org")
@ -115,19 +118,19 @@ int main() {
Computed data: Computed data:
|--------------------------+--------------------------------------+-----------------------------------------------------------------------------------------------| |--------------------------+----------------------------+-----------------------------------------------------------------------------------------------|
| ~coefficient_normalized~ | ~[prim_num]~ | Normalized primitive coefficients | | ~coefficient_normalized~ | ~[prim_num]~ | Normalized primitive coefficients |
| ~nucleus_prim_index~ | ~[nucl_num]~ | Index of the first primitive for each nucleus | | ~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_max_ang_mom~ | ~[nucl_num]~ | Maximum angular momentum for each nucleus |
| ~nucleus_range~ | ~[nucl_num]~ | Distance beyond which all the AOs are zero | | ~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~ | ~[5][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 | | ~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~ | ~[5][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 | | ~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~ | ~[5][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 | | ~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 | | ~nucl_shell_index~ | ~[nucl_num]~ | Index of the first shell for each nucleus |
| ~exponent_sorted~ | ~[prim_num]~ | Array of exponents for sorted primitives | | ~exponent_sorted~ | ~[prim_num]~ | Array of exponents for sorted primitives |
| ~coeff_norm_sorted~ | ~[prim_num]~ | Array of normalized coefficients for sorted primitives | | ~coeff_norm_sorted~ | ~[prim_num]~ | Array of normalized coefficients for sorted primitives |
@ -1836,7 +1839,7 @@ qmckl_exit_code qmckl_get_ao_basis_primitive_vgl(qmckl_context context, double*
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL); assert (ctx != NULL);
size_t sze = ctx->ao_basis.prim_num * 5 * ctx->electron.num * ctx->electron.walk_num; size_t sze = ctx->ao_basis.prim_num * 5 * ctx->electron.num;
memcpy(primitive_vgl, ctx->ao_basis.primitive_vgl, sze * sizeof(double)); memcpy(primitive_vgl, ctx->ao_basis.primitive_vgl, sze * sizeof(double));
return QMCKL_SUCCESS; return QMCKL_SUCCESS;
@ -1875,7 +1878,7 @@ qmckl_exit_code qmckl_provide_ao_basis_primitive_vgl(qmckl_context context)
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->ao_basis.prim_num * 5 * ctx->electron.num * mem_info.size = ctx->ao_basis.prim_num * 5 * ctx->electron.num *
ctx->electron.walk_num * sizeof(double); sizeof(double);
double* primitive_vgl = (double*) qmckl_malloc(context, mem_info); double* primitive_vgl = (double*) qmckl_malloc(context, mem_info);
if (primitive_vgl == NULL) { if (primitive_vgl == NULL) {
@ -1893,7 +1896,6 @@ qmckl_exit_code qmckl_provide_ao_basis_primitive_vgl(qmckl_context context)
ctx->ao_basis.prim_num, ctx->ao_basis.prim_num,
ctx->electron.num, ctx->electron.num,
ctx->nucleus.num, ctx->nucleus.num,
ctx->electron.walk_num,
ctx->ao_basis.nucleus_prim_index, ctx->ao_basis.nucleus_prim_index,
ctx->electron.coord_new, ctx->electron.coord_new,
ctx->nucleus.coord, ctx->nucleus.coord,
@ -1928,16 +1930,15 @@ qmckl_exit_code qmckl_provide_ao_basis_primitive_vgl(qmckl_context context)
| int64_t | prim_num | in | Number of primitives | | int64_t | prim_num | in | Number of primitives |
| int64_t | elec_num | in | Number of electrons | | int64_t | elec_num | in | Number of electrons |
| int64_t | nucl_num | in | Number of nuclei | | 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 | | 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 | elec_coord[3][elec_num] | in | Electron coordinates |
| double | nucl_coord[3][elec_num] | in | Nuclear coordinates | | double | nucl_coord[3][elec_num] | in | Nuclear coordinates |
| double | expo[prim_num] | in | Exponents of the primitives | | 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 | | double | primitive_vgl[5][elec_num][prim_num] | out | Value, gradients and Laplacian of the primitives |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes #+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_ao_basis_primitive_gaussian_vgl_f(context, & integer function qmckl_compute_ao_basis_primitive_gaussian_vgl_f(context, &
prim_num, elec_num, nucl_num, walk_num, & prim_num, elec_num, nucl_num, &
nucleus_prim_index, elec_coord, nucl_coord, expo, primitive_vgl) & nucleus_prim_index, elec_coord, nucl_coord, expo, primitive_vgl) &
result(info) result(info)
use qmckl use qmckl
@ -1946,14 +1947,13 @@ integer function qmckl_compute_ao_basis_primitive_gaussian_vgl_f(context, &
integer*8 , intent(in) :: prim_num integer*8 , intent(in) :: prim_num
integer*8 , intent(in) :: nucl_num integer*8 , intent(in) :: nucl_num
integer*8 , intent(in) :: elec_num integer*8 , intent(in) :: elec_num
integer*8 , intent(in) :: walk_num
integer*8 , intent(in) :: nucleus_prim_index(nucl_num+1) 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) :: elec_coord(elec_num,3)
double precision , intent(in) :: nucl_coord(nucl_num,3) double precision , intent(in) :: nucl_coord(nucl_num,3)
double precision , intent(in) :: expo(prim_num) double precision , intent(in) :: expo(prim_num)
double precision , intent(out) :: primitive_vgl(prim_num,elec_num,walk_num,5) double precision , intent(out) :: primitive_vgl(prim_num,elec_num,5)
integer*8 :: inucl, iprim, iwalk, ielec integer*8 :: inucl, iprim, ielec
double precision :: x, y, z, two_a, ar2, r2, v, cutoff double precision :: x, y, z, two_a, ar2, r2, v, cutoff
info = QMCKL_SUCCESS info = QMCKL_SUCCESS
@ -1964,11 +1964,10 @@ integer function qmckl_compute_ao_basis_primitive_gaussian_vgl_f(context, &
do inucl=1,nucl_num do inucl=1,nucl_num
! C is zero-based, so shift bounds by one ! C is zero-based, so shift bounds by one
do iprim = nucleus_prim_index(inucl)+1, nucleus_prim_index(inucl+1) do iprim = nucleus_prim_index(inucl)+1, nucleus_prim_index(inucl+1)
do iwalk = 1, walk_num
do ielec = 1, elec_num do ielec = 1, elec_num
x = elec_coord(ielec,1,iwalk) - nucl_coord(inucl,1) x = elec_coord(ielec,1) - nucl_coord(inucl,1)
y = elec_coord(ielec,2,iwalk) - nucl_coord(inucl,2) y = elec_coord(ielec,2) - nucl_coord(inucl,2)
z = elec_coord(ielec,3,iwalk) - nucl_coord(inucl,3) z = elec_coord(ielec,3) - nucl_coord(inucl,3)
r2 = x*x + y*y + z*z r2 = x*x + y*y + z*z
ar2 = expo(iprim)*r2 ar2 = expo(iprim)*r2
@ -1977,16 +1976,15 @@ integer function qmckl_compute_ao_basis_primitive_gaussian_vgl_f(context, &
v = dexp(-ar2) v = dexp(-ar2)
two_a = -2.d0 * expo(iprim) * v two_a = -2.d0 * expo(iprim) * v
primitive_vgl(iprim, ielec, iwalk, 1) = v primitive_vgl(iprim, ielec, 1) = v
primitive_vgl(iprim, ielec, iwalk, 2) = two_a * x primitive_vgl(iprim, ielec, 2) = two_a * x
primitive_vgl(iprim, ielec, iwalk, 3) = two_a * y primitive_vgl(iprim, ielec, 3) = two_a * y
primitive_vgl(iprim, ielec, iwalk, 4) = two_a * z primitive_vgl(iprim, ielec, 4) = two_a * z
primitive_vgl(iprim, ielec, iwalk, 5) = two_a * (3.d0 - 2.d0*ar2) primitive_vgl(iprim, ielec, 5) = two_a * (3.d0 - 2.d0*ar2)
end do end do
end do end do
end do end do
end do
end function qmckl_compute_ao_basis_primitive_gaussian_vgl_f end function qmckl_compute_ao_basis_primitive_gaussian_vgl_f
#+end_src #+end_src
@ -1997,7 +1995,6 @@ qmckl_exit_code qmckl_compute_ao_basis_primitive_gaussian_vgl(
const int64_t prim_num, const int64_t prim_num,
const int64_t elec_num, const int64_t elec_num,
const int64_t nucl_num, const int64_t nucl_num,
const int64_t walk_num,
const int64_t* nucleus_prim_index, const int64_t* nucleus_prim_index,
const double* elec_coord, const double* elec_coord,
const double* nucl_coord, const double* nucl_coord,
@ -2010,7 +2007,15 @@ qmckl_exit_code qmckl_compute_ao_basis_primitive_gaussian_vgl(
#+RESULTS: #+RESULTS:
#+begin_src f90 :tangle (eval f) :comments org :exports none #+begin_src f90 :tangle (eval f) :comments org :exports none
integer(c_int32_t) function qmckl_compute_ao_basis_primitive_gaussian_vgl & integer(c_int32_t) function qmckl_compute_ao_basis_primitive_gaussian_vgl &
(context, prim_num, elec_num, nucl_num, walk_num, nucleus_prim_index, elec_coord, nucl_coord, expo, primitive_vgl) & (context, &
prim_num, &
elec_num, &
nucl_num, &
nucleus_prim_index, &
elec_coord, &
nucl_coord, &
expo, &
primitive_vgl) &
bind(C) result(info) bind(C) result(info)
use, intrinsic :: iso_c_binding use, intrinsic :: iso_c_binding
@ -2020,16 +2025,23 @@ qmckl_exit_code qmckl_compute_ao_basis_primitive_gaussian_vgl(
integer (c_int64_t) , intent(in) , value :: prim_num integer (c_int64_t) , intent(in) , value :: prim_num
integer (c_int64_t) , intent(in) , value :: elec_num integer (c_int64_t) , intent(in) , value :: elec_num
integer (c_int64_t) , intent(in) , value :: nucl_num integer (c_int64_t) , intent(in) , value :: nucl_num
integer (c_int64_t) , intent(in) , value :: walk_num
integer (c_int64_t) , intent(in) :: nucleus_prim_index(nucl_num) integer (c_int64_t) , intent(in) :: nucleus_prim_index(nucl_num)
real (c_double ) , intent(in) :: elec_coord(elec_num,3,walk_num) real (c_double ) , intent(in) :: elec_coord(elec_num,3)
real (c_double ) , intent(in) :: nucl_coord(elec_num,3) real (c_double ) , intent(in) :: nucl_coord(elec_num,3)
real (c_double ) , intent(in) :: expo(prim_num) real (c_double ) , intent(in) :: expo(prim_num)
real (c_double ) , intent(out) :: primitive_vgl(elec_num,walk_num,5,prim_num) real (c_double ) , intent(out) :: primitive_vgl(prim_num,elec_num,5)
integer(c_int32_t), external :: qmckl_compute_ao_basis_primitive_gaussian_vgl_f integer(c_int32_t), external :: qmckl_compute_ao_basis_primitive_gaussian_vgl_f
info = qmckl_compute_ao_basis_primitive_gaussian_vgl_f & info = 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) (context, &
prim_num, &
elec_num, &
nucl_num, &
nucleus_prim_index, &
elec_coord, &
nucl_coord, &
expo, &
primitive_vgl)
end function qmckl_compute_ao_basis_primitive_gaussian_vgl end function qmckl_compute_ao_basis_primitive_gaussian_vgl
#+end_src #+end_src
@ -2058,40 +2070,25 @@ def lf(a,x,y):
return d2f(a,x,y,1) + d2f(a,x,y,2) + d2f(a,x,y,3) return d2f(a,x,y,1) + d2f(a,x,y,2) + d2f(a,x,y,3)
elec_26_w1 = np.array( [ 1.49050402641, 2.90106987953, -1.05920815468 ] ) elec_26_w1 = np.array( [ 1.49050402641, 2.90106987953, -1.05920815468 ] )
elec_15_w2 = np.array( [ -2.20180344582,-1.9113150239, 2.2193744778600002 ] )
nucl_1 = np.array( [ 1.096243353458458e+00, 8.907054016973815e-01, 7.777092280258892e-01 ] ) nucl_1 = np.array( [ 1.096243353458458e+00, 8.907054016973815e-01, 7.777092280258892e-01 ] )
nucl_2 = np.array( [ 1.168459237342663e+00, 1.125660720053393e+00, 2.833370314829343e+00 ] ) nucl_2 = np.array( [ 1.168459237342663e+00, 1.125660720053393e+00, 2.833370314829343e+00 ] )
#double prim_vgl[prim_num][5][walk_num][elec_num]; #double prim_vgl[prim_num][5][elec_num];
a = 0.9059; x = elec_26_w1 ; y = nucl_1 a = 0.9059; x = elec_26_w1 ; y = nucl_1
print ( "[7][0][0][26] : %e"% f(a,x,y)) print ( "[7][0][26] : %e"% f(a,x,y))
print ( "[7][1][0][26] : %e"% df(a,x,y,1)) print ( "[7][1][26] : %e"% df(a,x,y,1))
print ( "[7][2][0][26] : %e"% df(a,x,y,2)) print ( "[7][2][26] : %e"% df(a,x,y,2))
print ( "[7][3][0][26] : %e"% df(a,x,y,3)) print ( "[7][3][26] : %e"% df(a,x,y,3))
print ( "[7][4][0][26] : %e"% lf(a,x,y)) print ( "[7][4][26] : %e"% lf(a,x,y))
a = 0.32578; x = elec_15_w2 ; y = nucl_2
print ( "[39][0][1][15] : %e"% f(a,x,y))
print ( "[39][1][1][15] : %e"% df(a,x,y,1))
print ( "[39][2][1][15] : %e"% df(a,x,y,2))
print ( "[39][3][1][15] : %e"% df(a,x,y,3))
print ( "[39][4][1][15] : %e"% lf(a,x,y))
#+end_src #+end_src
#+RESULTS: #+RESULTS:
#+begin_example : [7][0][26] : 1.050157e-03
[7][0][0][26] : 1.050157e-03 : [7][1][26] : -7.501497e-04
[7][1][0][26] : -7.501497e-04 : [7][2][26] : -3.825069e-03
[7][2][0][26] : -3.825069e-03 : [7][3][26] : 3.495056e-03
[7][3][0][26] : 3.495056e-03 : [7][4][26] : 2.040013e-02
[7][4][0][26] : 2.040013e-02
[39][0][1][15] : 1.083038e-03
[39][1][1][15] : 2.378275e-03
[39][2][1][15] : 2.143086e-03
[39][3][1][15] : 4.332750e-04
[39][4][1][15] : 7.514605e-03
#+end_example
*** Test *** Test
@ -2118,22 +2115,16 @@ assert(rc == QMCKL_SUCCESS);
double prim_vgl[5][walk_num][elec_num][prim_num]; double prim_vgl[5][elec_num][prim_num];
rc = qmckl_get_ao_basis_primitive_vgl(context, &(prim_vgl[0][0][0][0])); rc = qmckl_get_ao_basis_primitive_vgl(context, &(prim_vgl[0][0][0]));
assert (rc == QMCKL_SUCCESS); assert (rc == QMCKL_SUCCESS);
assert( fabs(prim_vgl[0][0][26][7] - ( 1.0501570432064878E-003)) < 1.e-14 ); assert( fabs(prim_vgl[0][26][7] - ( 1.0501570432064878E-003)) < 1.e-14 );
assert( fabs(prim_vgl[1][0][26][7] - (-7.5014974095310560E-004)) < 1.e-14 ); assert( fabs(prim_vgl[1][26][7] - (-7.5014974095310560E-004)) < 1.e-14 );
assert( fabs(prim_vgl[2][0][26][7] - (-3.8250692897610380E-003)) < 1.e-14 ); assert( fabs(prim_vgl[2][26][7] - (-3.8250692897610380E-003)) < 1.e-14 );
assert( fabs(prim_vgl[3][0][26][7] - ( 3.4950559194080275E-003)) < 1.e-14 ); assert( fabs(prim_vgl[3][26][7] - ( 3.4950559194080275E-003)) < 1.e-14 );
assert( fabs(prim_vgl[4][0][26][7] - ( 2.0392163767356572E-002)) < 1.e-14 ); assert( fabs(prim_vgl[4][26][7] - ( 2.0392163767356572E-002)) < 1.e-14 );
assert( fabs(prim_vgl[0][1][15][39] - ( 1.0825844173157661E-003)) < 1.e-14 );
assert( fabs(prim_vgl[1][1][15][39] - ( 2.3774237611651531E-003)) < 1.e-14 );
assert( fabs(prim_vgl[2][1][15][39] - ( 2.1423191526963063E-003)) < 1.e-14 );
assert( fabs(prim_vgl[3][1][15][39] - ( 4.3312003523048492E-004)) < 1.e-14 );
assert( fabs(prim_vgl[4][1][15][39] - ( 7.5174404780004771E-003)) < 1.e-14 );
} }
@ -2143,13 +2134,11 @@ assert( fabs(prim_vgl[4][1][15][39] - ( 7.5174404780004771E-003)) < 1.e-14 );
*** Ideas for improvement *** Ideas for improvement
#+begin_src c #+begin_src c
// m : walkers
// j : electrons // j : electrons
// l : primitives // l : primitives
k=0; k=0;
for (m=0 ; m<walk_num ; ++m) { for (j=0 ; j<elec_num ; ++j) {
for (j=0 ; j<elec_num ; ++j) {
for (i=0 ; i<nucl_num ; ++i) { for (i=0 ; i<nucl_num ; ++i) {
r2 = nucl_elec_dist[i][j]; r2 = nucl_elec_dist[i][j];
@ -2166,7 +2155,6 @@ for (m=0 ; m<walk_num ; ++m) {
} }
} }
} }
}
// sort(tmp) in increasing ar2; // sort(tmp) in increasing ar2;
// Identify first ar2 above numerical accuracy threshold // Identify first ar2 above numerical accuracy threshold
@ -2196,7 +2184,7 @@ qmckl_exit_code qmckl_get_ao_basis_shell_vgl(qmckl_context context, double* cons
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL); assert (ctx != NULL);
size_t sze = ctx->ao_basis.shell_num * 5 * ctx->electron.num * ctx->electron.walk_num; size_t sze = ctx->ao_basis.shell_num * 5 * ctx->electron.num;
memcpy(shell_vgl, ctx->ao_basis.shell_vgl, sze * sizeof(double)); memcpy(shell_vgl, ctx->ao_basis.shell_vgl, sze * sizeof(double));
return QMCKL_SUCCESS; return QMCKL_SUCCESS;
@ -2255,8 +2243,7 @@ qmckl_exit_code qmckl_provide_ao_basis_shell_vgl(qmckl_context context)
if (ctx->ao_basis.shell_vgl == NULL) { if (ctx->ao_basis.shell_vgl == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->ao_basis.shell_num * 5 * ctx->electron.num * mem_info.size = ctx->ao_basis.shell_num * 5 * ctx->electron.num * sizeof(double);
ctx->electron.walk_num * sizeof(double);
double* shell_vgl = (double*) qmckl_malloc(context, mem_info); double* shell_vgl = (double*) qmckl_malloc(context, mem_info);
if (shell_vgl == NULL) { if (shell_vgl == NULL) {
@ -2275,7 +2262,6 @@ qmckl_exit_code qmckl_provide_ao_basis_shell_vgl(qmckl_context context)
ctx->ao_basis.shell_num, ctx->ao_basis.shell_num,
ctx->electron.num, ctx->electron.num,
ctx->nucleus.num, ctx->nucleus.num,
ctx->electron.walk_num,
ctx->ao_basis.nucleus_shell_num, ctx->ao_basis.nucleus_shell_num,
ctx->ao_basis.nucleus_index, ctx->ao_basis.nucleus_index,
ctx->ao_basis.shell_prim_index, ctx->ao_basis.shell_prim_index,
@ -2315,20 +2301,19 @@ qmckl_exit_code qmckl_provide_ao_basis_shell_vgl(qmckl_context context)
| ~int64_t~ | ~shell_num~ | in | Number of shells | | ~int64_t~ | ~shell_num~ | in | Number of shells |
| ~int64_t~ | ~elec_num~ | in | Number of electrons | | ~int64_t~ | ~elec_num~ | in | Number of electrons |
| ~int64_t~ | ~nucl_num~ | in | Number of nuclei | | ~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_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~ | ~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_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 | | ~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~ | ~elec_coord[3][elec_num]~ | in | Electron coordinates |
| ~double~ | ~nucl_coord[3][elec_num]~ | in | Nuclear coordinates | | ~double~ | ~nucl_coord[3][elec_num]~ | in | Nuclear coordinates |
| ~double~ | ~expo[prim_num]~ | in | Exponents of the primitives | | ~double~ | ~expo[prim_num]~ | in | Exponents of the primitives |
| ~double~ | ~coef_normalized[prim_num]~ | in | Coefficients 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 | | ~double~ | ~shell_vgl[5][elec_num][shell_num]~ | out | Value, gradients and Laplacian of the shells |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes #+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_ao_basis_shell_gaussian_vgl_f(context, & integer function qmckl_compute_ao_basis_shell_gaussian_vgl_f(context, &
prim_num, shell_num, elec_num, nucl_num, walk_num, & prim_num, shell_num, elec_num, nucl_num, &
nucleus_shell_num, nucleus_index, shell_prim_index, shell_prim_num, & nucleus_shell_num, nucleus_index, shell_prim_index, shell_prim_num, &
elec_coord, nucl_coord, expo, coef_normalized, shell_vgl) & elec_coord, nucl_coord, expo, coef_normalized, shell_vgl) &
result(info) result(info)
@ -2339,18 +2324,17 @@ integer function qmckl_compute_ao_basis_shell_gaussian_vgl_f(context, &
integer*8 , intent(in) :: shell_num integer*8 , intent(in) :: shell_num
integer*8 , intent(in) :: nucl_num integer*8 , intent(in) :: nucl_num
integer*8 , intent(in) :: elec_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_shell_num(nucl_num)
integer*8 , intent(in) :: nucleus_index(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_index(shell_num)
integer*8 , intent(in) :: shell_prim_num(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) :: elec_coord(elec_num,3)
double precision , intent(in) :: nucl_coord(nucl_num,3) double precision , intent(in) :: nucl_coord(nucl_num,3)
double precision , intent(in) :: expo(prim_num) double precision , intent(in) :: expo(prim_num)
double precision , intent(in) :: coef_normalized(prim_num) double precision , intent(in) :: coef_normalized(prim_num)
double precision , intent(out) :: shell_vgl(shell_num,elec_num,walk_num,5) double precision , intent(out) :: shell_vgl(shell_num,elec_num,5)
integer*8 :: inucl, iprim, iwalk, ielec, ishell integer*8 :: inucl, iprim, ielec, ishell
double precision :: x, y, z, two_a, ar2, r2, v, cutoff double precision :: x, y, z, two_a, ar2, r2, v, cutoff
info = QMCKL_SUCCESS info = QMCKL_SUCCESS
@ -2361,23 +2345,22 @@ integer function qmckl_compute_ao_basis_shell_gaussian_vgl_f(context, &
do inucl=1,nucl_num do inucl=1,nucl_num
do iwalk = 1, walk_num
do ielec = 1, elec_num do ielec = 1, elec_num
x = elec_coord(ielec,1,iwalk) - nucl_coord(inucl,1) x = elec_coord(ielec,1) - nucl_coord(inucl,1)
y = elec_coord(ielec,2,iwalk) - nucl_coord(inucl,2) y = elec_coord(ielec,2) - nucl_coord(inucl,2)
z = elec_coord(ielec,3,iwalk) - nucl_coord(inucl,3) z = elec_coord(ielec,3) - nucl_coord(inucl,3)
r2 = x*x + y*y + z*z r2 = x*x + y*y + z*z
do ishell=nucleus_index(inucl)+1, nucleus_index(inucl)+nucleus_shell_num(inucl) do ishell=nucleus_index(inucl)+1, nucleus_index(inucl)+nucleus_shell_num(inucl)
! C is zero-based, so shift bounds by one ! C is zero-based, so shift bounds by one
shell_vgl(ishell, ielec, iwalk, 1) = 0.d0 shell_vgl(ishell, ielec, 1) = 0.d0
shell_vgl(ishell, ielec, iwalk, 2) = 0.d0 shell_vgl(ishell, ielec, 2) = 0.d0
shell_vgl(ishell, ielec, iwalk, 3) = 0.d0 shell_vgl(ishell, ielec, 3) = 0.d0
shell_vgl(ishell, ielec, iwalk, 4) = 0.d0 shell_vgl(ishell, ielec, 4) = 0.d0
shell_vgl(ishell, ielec, iwalk, 5) = 0.d0 shell_vgl(ishell, ielec, 5) = 0.d0
do iprim = shell_prim_index(ishell)+1, shell_prim_index(ishell)+shell_prim_num(ishell) do iprim = shell_prim_index(ishell)+1, shell_prim_index(ishell)+shell_prim_num(ishell)
@ -2389,20 +2372,20 @@ integer function qmckl_compute_ao_basis_shell_gaussian_vgl_f(context, &
v = coef_normalized(iprim) * dexp(-ar2) v = coef_normalized(iprim) * dexp(-ar2)
two_a = -2.d0 * expo(iprim) * v two_a = -2.d0 * expo(iprim) * v
shell_vgl(ishell, ielec, iwalk, 1) = & shell_vgl(ishell, ielec, 1) = &
shell_vgl(ishell, ielec, iwalk, 1) + v shell_vgl(ishell, ielec, 1) + v
shell_vgl(ishell, ielec, iwalk, 2) = & shell_vgl(ishell, ielec, 2) = &
shell_vgl(ishell, ielec, iwalk, 2) + two_a * x shell_vgl(ishell, ielec, 2) + two_a * x
shell_vgl(ishell, ielec, iwalk, 3) = & shell_vgl(ishell, ielec, 3) = &
shell_vgl(ishell, ielec, iwalk, 3) + two_a * y shell_vgl(ishell, ielec, 3) + two_a * y
shell_vgl(ishell, ielec, iwalk, 4) = & shell_vgl(ishell, ielec, 4) = &
shell_vgl(ishell, ielec, iwalk, 4) + two_a * z shell_vgl(ishell, ielec, 4) + two_a * z
shell_vgl(ishell, ielec, iwalk, 5) = & shell_vgl(ishell, ielec, 5) = &
shell_vgl(ishell, ielec, iwalk, 5) + two_a * (3.d0 - 2.d0*ar2) shell_vgl(ishell, ielec, 5) + two_a * (3.d0 - 2.d0*ar2)
end do end do
@ -2410,7 +2393,6 @@ integer function qmckl_compute_ao_basis_shell_gaussian_vgl_f(context, &
end do end do
end do end do
end do
end function qmckl_compute_ao_basis_shell_gaussian_vgl_f end function qmckl_compute_ao_basis_shell_gaussian_vgl_f
#+end_src #+end_src
@ -2425,7 +2407,6 @@ end function qmckl_compute_ao_basis_shell_gaussian_vgl_f
const int64_t shell_num, const int64_t shell_num,
const int64_t elec_num, const int64_t elec_num,
const int64_t nucl_num, const int64_t nucl_num,
const int64_t walk_num,
const int64_t* nucleus_shell_num, const int64_t* nucleus_shell_num,
const int64_t* nucleus_index, const int64_t* nucleus_index,
const int64_t* shell_prim_index, const int64_t* shell_prim_index,
@ -2447,7 +2428,6 @@ end function qmckl_compute_ao_basis_shell_gaussian_vgl_f
shell_num, & shell_num, &
elec_num, & elec_num, &
nucl_num, & nucl_num, &
walk_num, &
nucleus_shell_num, & nucleus_shell_num, &
nucleus_index, & nucleus_index, &
shell_prim_index, & shell_prim_index, &
@ -2467,16 +2447,15 @@ end function qmckl_compute_ao_basis_shell_gaussian_vgl_f
integer (c_int64_t) , intent(in) , value :: shell_num integer (c_int64_t) , intent(in) , value :: shell_num
integer (c_int64_t) , intent(in) , value :: elec_num integer (c_int64_t) , intent(in) , value :: elec_num
integer (c_int64_t) , intent(in) , value :: nucl_num integer (c_int64_t) , intent(in) , value :: nucl_num
integer (c_int64_t) , intent(in) , value :: walk_num
integer (c_int64_t) , intent(in) :: nucleus_shell_num(nucl_num) integer (c_int64_t) , intent(in) :: nucleus_shell_num(nucl_num)
integer (c_int64_t) , intent(in) :: nucleus_index(nucl_num) integer (c_int64_t) , intent(in) :: nucleus_index(nucl_num)
integer (c_int64_t) , intent(in) :: shell_prim_index(shell_num) integer (c_int64_t) , intent(in) :: shell_prim_index(shell_num)
integer (c_int64_t) , intent(in) :: shell_prim_num(shell_num) integer (c_int64_t) , intent(in) :: shell_prim_num(shell_num)
real (c_double ) , intent(in) :: elec_coord(elec_num,3,walk_num) real (c_double ) , intent(in) :: elec_coord(elec_num,3)
real (c_double ) , intent(in) :: nucl_coord(elec_num,3) real (c_double ) , intent(in) :: nucl_coord(elec_num,3)
real (c_double ) , intent(in) :: expo(prim_num) real (c_double ) , intent(in) :: expo(prim_num)
real (c_double ) , intent(in) :: coef_normalized(prim_num) real (c_double ) , intent(in) :: coef_normalized(prim_num)
real (c_double ) , intent(out) :: shell_vgl(shell_num,elec_num,walk_num,5) real (c_double ) , intent(out) :: shell_vgl(shell_num,elec_num,5)
integer(c_int32_t), external :: qmckl_compute_ao_basis_shell_gaussian_vgl_f integer(c_int32_t), external :: qmckl_compute_ao_basis_shell_gaussian_vgl_f
info = qmckl_compute_ao_basis_shell_gaussian_vgl_f & info = qmckl_compute_ao_basis_shell_gaussian_vgl_f &
@ -2485,7 +2464,6 @@ end function qmckl_compute_ao_basis_shell_gaussian_vgl_f
shell_num, & shell_num, &
elec_num, & elec_num, &
nucl_num, & nucl_num, &
walk_num, &
nucleus_shell_num, & nucleus_shell_num, &
nucleus_index, & nucleus_index, &
shell_prim_index, & shell_prim_index, &
@ -2527,7 +2505,7 @@ elec_15_w2 = np.array( [ -2.20180344582,-1.9113150239, 2.2193744778600002 ] )
nucl_1 = np.array( [ 1.096243353458458e+00, 8.907054016973815e-01, 7.777092280258892e-01 ] ) nucl_1 = np.array( [ 1.096243353458458e+00, 8.907054016973815e-01, 7.777092280258892e-01 ] )
nucl_2 = np.array( [ 1.168459237342663e+00, 1.125660720053393e+00, 2.833370314829343e+00 ] ) nucl_2 = np.array( [ 1.168459237342663e+00, 1.125660720053393e+00, 2.833370314829343e+00 ] )
#double prim_vgl[prim_num][5][walk_num][elec_num]; #double prim_vgl[prim_num][5][elec_num];
x = elec_26_w1 ; y = nucl_1 x = elec_26_w1 ; y = nucl_1
a = [( 8.236000E+03, -1.130000E-04 * 6.1616545431994848e+02 ), a = [( 8.236000E+03, -1.130000E-04 * 6.1616545431994848e+02 ),
( 1.235000E+03, -8.780000E-04 * 1.4847738511079908e+02 ), ( 1.235000E+03, -8.780000E-04 * 1.4847738511079908e+02 ),
@ -2540,40 +2518,20 @@ a = [( 8.236000E+03, -1.130000E-04 * 6.1616545431994848e+02 ),
( 3.643000E-01, 5.986840E-01 * 3.3419848027174592e-01 ), ( 3.643000E-01, 5.986840E-01 * 3.3419848027174592e-01 ),
( 1.285000E-01, 3.953890E-01 * 1.5296336817449557e-01 )] ( 1.285000E-01, 3.953890E-01 * 1.5296336817449557e-01 )]
print ( "[1][0][0][26] : %25.15e"% f(a,x,y)) print ( "[1][0][26] : %25.15e"% f(a,x,y))
print ( "[1][1][0][26] : %25.15e"% df(a,x,y,1)) print ( "[1][1][26] : %25.15e"% df(a,x,y,1))
print ( "[1][2][0][26] : %25.15e"% df(a,x,y,2)) print ( "[1][2][26] : %25.15e"% df(a,x,y,2))
print ( "[1][3][0][26] : %25.15e"% df(a,x,y,3)) print ( "[1][3][26] : %25.15e"% df(a,x,y,3))
print ( "[1][4][0][26] : %25.15e"% lf(a,x,y)) print ( "[1][4][26] : %25.15e"% lf(a,x,y))
x = elec_15_w2 ; y = nucl_2
a = [(3.387000E+01, 6.068000E-03 *1.0006253235944540e+01),
(5.095000E+00, 4.530800E-02 *2.4169531573445120e+00),
(1.159000E+00, 2.028220E-01 *7.9610924849766440e-01),
(3.258000E-01, 5.039030E-01 *3.0734305383061117e-01),
(1.027000E-01, 3.834210E-01 *1.2929684417481876e-01)]
print ( "[0][1][15][14] : %25.15e"% f(a,x,y))
print ( "[1][1][15][14] : %25.15e"% df(a,x,y,1))
print ( "[2][1][15][14] : %25.15e"% df(a,x,y,2))
print ( "[3][1][15][14] : %25.15e"% df(a,x,y,3))
print ( "[4][1][15][14] : %25.15e"% lf(a,x,y))
#+end_src #+end_src
#+RESULTS: #+RESULTS:
#+begin_example : [1][0][26] : 3.564393437193867e-02
[1][0][0][26] : 3.564393437193867e-02 : [1][1][26] : -6.030177988891605e-03
[1][1][0][26] : -6.030177988891605e-03 : [1][2][26] : -3.074832579871845e-02
[1][2][0][26] : -3.074832579871845e-02 : [1][3][26] : 2.809546963133958e-02
[1][3][0][26] : 2.809546963133958e-02 : [1][4][26] : 1.903338597841753e-02
[1][4][0][26] : 1.903338597841753e-02
[0][1][15][14] : 5.928089771361000e-03
[1][1][15][14] : 4.355862298893037e-03
[2][1][15][14] : 3.925108924950765e-03
[3][1][15][14] : 7.935527764416084e-04
[4][1][15][14] : 2.697495005143935e-03
#+end_example
*** Test *** Test
@ -2599,34 +2557,23 @@ rc = qmckl_set_electron_coord (context, 'N', elec_coord);
assert(rc == QMCKL_SUCCESS); assert(rc == QMCKL_SUCCESS);
double shell_vgl[5][walk_num][elec_num][shell_num]; double shell_vgl[5][elec_num][shell_num];
rc = qmckl_get_ao_basis_shell_vgl(context, &(shell_vgl[0][0][0][0])); rc = qmckl_get_ao_basis_shell_vgl(context, &(shell_vgl[0][0][0]));
assert (rc == QMCKL_SUCCESS); assert (rc == QMCKL_SUCCESS);
printf(" shell_vgl[1][0][0][26] %25.15e\n", shell_vgl[0][0][26][1]); printf(" shell_vgl[1][0][26] %25.15e\n", shell_vgl[0][26][1]);
printf(" shell_vgl[1][1][0][26] %25.15e\n", shell_vgl[1][0][26][1]); printf(" shell_vgl[1][1][26] %25.15e\n", shell_vgl[1][26][1]);
printf(" shell_vgl[1][2][0][26] %25.15e\n", shell_vgl[2][0][26][1]); printf(" shell_vgl[1][2][26] %25.15e\n", shell_vgl[2][26][1]);
printf(" shell_vgl[1][3][0][26] %25.15e\n", shell_vgl[3][0][26][1]); printf(" shell_vgl[1][3][26] %25.15e\n", shell_vgl[3][26][1]);
printf(" shell_vgl[1][4][0][26] %25.15e\n", shell_vgl[4][0][26][1]); printf(" shell_vgl[1][4][26] %25.15e\n", shell_vgl[4][26][1]);
printf(" shell_vgl[14][0][1][15] %25.15e\n", shell_vgl[0][1][15][14]); assert( fabs(shell_vgl[0][26][1] - ( 3.564393437193868e-02)) < 1.e-14 );
printf(" shell_vgl[14][1][1][15] %25.15e\n", shell_vgl[1][1][15][14]); assert( fabs(shell_vgl[1][26][1] - (-6.030177987072189e-03)) < 1.e-14 );
printf(" shell_vgl[14][2][1][15] %25.15e\n", shell_vgl[2][1][15][14]); assert( fabs(shell_vgl[2][26][1] - (-3.074832579537582e-02)) < 1.e-14 );
printf(" shell_vgl[14][3][1][15] %25.15e\n", shell_vgl[3][1][15][14]); assert( fabs(shell_vgl[3][26][1] - ( 2.809546963519935e-02)) < 1.e-14 );
printf(" shell_vgl[14][4][1][15] %25.15e\n", shell_vgl[4][1][15][14]); assert( fabs(shell_vgl[4][26][1] - ( 1.896046117183968e-02)) < 1.e-14 );
assert( fabs(shell_vgl[0][0][26][1] - ( 3.564393437193868e-02)) < 1.e-14 );
assert( fabs(shell_vgl[1][0][26][1] - (-6.030177987072189e-03)) < 1.e-14 );
assert( fabs(shell_vgl[2][0][26][1] - (-3.074832579537582e-02)) < 1.e-14 );
assert( fabs(shell_vgl[3][0][26][1] - ( 2.809546963519935e-02)) < 1.e-14 );
assert( fabs(shell_vgl[4][0][26][1] - ( 1.896046117183968e-02)) < 1.e-14 );
assert( fabs(shell_vgl[0][1][15][14] - ( 5.928089771361000e-03)) < 1.e-14 );
assert( fabs(shell_vgl[1][1][15][14] - ( 4.355862296021654e-03)) < 1.e-14 );
assert( fabs(shell_vgl[2][1][15][14] - ( 3.925108924923650e-03)) < 1.e-14 );
assert( fabs(shell_vgl[3][1][15][14] - ( 7.935527784022099e-04)) < 1.e-14 );
assert( fabs(shell_vgl[4][1][15][14] - ( 2.708246573703548e-03)) < 1.e-14 );
} }
#+end_src #+end_src
@ -3242,7 +3189,7 @@ qmckl_exit_code qmckl_get_ao_vgl(qmckl_context context, double* const ao_vgl) {
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL); assert (ctx != NULL);
size_t sze = ctx->ao_basis.ao_num * 5 * ctx->electron.num * ctx->electron.walk_num; size_t sze = ctx->ao_basis.ao_num * 5 * ctx->electron.num;
memcpy(ao_vgl, ctx->ao_basis.ao_vgl, sze * sizeof(double)); memcpy(ao_vgl, ctx->ao_basis.ao_vgl, sze * sizeof(double));
return QMCKL_SUCCESS; return QMCKL_SUCCESS;
@ -3309,8 +3256,7 @@ qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context)
if (ctx->ao_basis.ao_vgl == NULL) { if (ctx->ao_basis.ao_vgl == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->ao_basis.ao_num * 5 * ctx->electron.num * mem_info.size = ctx->ao_basis.ao_num * 5 * ctx->electron.num * sizeof(double);
ctx->electron.walk_num * sizeof(double);
double* ao_vgl = (double*) qmckl_malloc(context, mem_info); double* ao_vgl = (double*) qmckl_malloc(context, mem_info);
if (ao_vgl == NULL) { if (ao_vgl == NULL) {
@ -3327,7 +3273,6 @@ qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context)
ctx->ao_basis.shell_num, ctx->ao_basis.shell_num,
ctx->electron.num, ctx->electron.num,
ctx->nucleus.num, ctx->nucleus.num,
ctx->electron.walk_num,
ctx->electron.coord_new, ctx->electron.coord_new,
ctx->nucleus.coord, ctx->nucleus.coord,
ctx->ao_basis.nucleus_index, ctx->ao_basis.nucleus_index,
@ -3362,8 +3307,7 @@ qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context)
| ~int64_t~ | ~shell_num~ | in | Number of shells | | ~int64_t~ | ~shell_num~ | in | Number of shells |
| ~int64_t~ | ~elec_num~ | in | Number of electrons | | ~int64_t~ | ~elec_num~ | in | Number of electrons |
| ~int64_t~ | ~nucl_num~ | in | Number of nuclei | | ~int64_t~ | ~nucl_num~ | in | Number of nuclei |
| ~int64_t~ | ~walk_num~ | in | Number of walkers | | ~double~ | ~elec_coord[3][elec_num]~ | in | Electron coordinates |
| ~double~ | ~elec_coord[walk_num][3][elec_num]~ | in | Electron coordinates |
| ~double~ | ~nucl_coord[3][nucl_num]~ | in | Nuclear 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_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 | | ~int64_t~ | ~nucleus_shell_num[nucl_num]~ | in | Number of shells per nucleus |
@ -3371,12 +3315,12 @@ qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context)
| ~int32_t~ | ~nucleus_max_ang_mom[nucl_num]~ | in | Maximum angular momentum per nucleus | | ~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 | | ~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~ | ~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~ | ~shell_vgl[5][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 | | ~double~ | ~ao_vgl[5][elec_num][ao_num]~ | out | Value, gradients and Laplacian of the AOs |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes #+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_ao_vgl_f(context, & integer function qmckl_compute_ao_vgl_f(context, &
ao_num, shell_num, elec_num, nucl_num, walk_num, & ao_num, shell_num, elec_num, nucl_num, &
elec_coord, nucl_coord, nucleus_index, nucleus_shell_num, & elec_coord, nucl_coord, nucleus_index, nucleus_shell_num, &
nucleus_range, nucleus_max_ang_mom, shell_ang_mom, & nucleus_range, nucleus_max_ang_mom, shell_ang_mom, &
ao_factor, shell_vgl, ao_vgl) & ao_factor, shell_vgl, ao_vgl) &
@ -3388,8 +3332,7 @@ integer function qmckl_compute_ao_vgl_f(context, &
integer*8 , intent(in) :: shell_num integer*8 , intent(in) :: shell_num
integer*8 , intent(in) :: elec_num integer*8 , intent(in) :: elec_num
integer*8 , intent(in) :: nucl_num integer*8 , intent(in) :: nucl_num
integer*8 , intent(in) :: walk_num double precision , intent(in) :: elec_coord(elec_num,3)
double precision , intent(in) :: elec_coord(elec_num,3,walk_num)
double precision , intent(in) :: nucl_coord(nucl_num,3) double precision , intent(in) :: nucl_coord(nucl_num,3)
integer*8 , intent(in) :: nucleus_index(nucl_num) integer*8 , intent(in) :: nucleus_index(nucl_num)
integer*8 , intent(in) :: nucleus_shell_num(nucl_num) integer*8 , intent(in) :: nucleus_shell_num(nucl_num)
@ -3397,13 +3340,13 @@ integer function qmckl_compute_ao_vgl_f(context, &
integer , intent(in) :: nucleus_max_ang_mom(nucl_num) integer , intent(in) :: nucleus_max_ang_mom(nucl_num)
integer , intent(in) :: shell_ang_mom(shell_num) integer , intent(in) :: shell_ang_mom(shell_num)
double precision , intent(in) :: ao_factor(ao_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(in) :: shell_vgl(shell_num,elec_num,5)
double precision , intent(out) :: ao_vgl(ao_num,elec_num,walk_num,5) double precision , intent(out) :: ao_vgl(ao_num,elec_num,5)
double precision :: e_coord(3), n_coord(3) double precision :: e_coord(3), n_coord(3)
integer*8 :: n_poly integer*8 :: n_poly
integer :: l, il, k integer :: l, il, k
integer*8 :: ielec, inucl, ishell, iwalk integer*8 :: ielec, inucl, ishell
integer :: lstart(0:20) integer :: lstart(0:20)
double precision :: x, y, z, r2 double precision :: x, y, z, r2
double precision :: cutoff double precision :: cutoff
@ -3425,11 +3368,10 @@ integer function qmckl_compute_ao_vgl_f(context, &
! TODO : Use numerical precision here ! TODO : Use numerical precision here
cutoff = -dlog(1.d-15) cutoff = -dlog(1.d-15)
do iwalk = 1,walk_num
do ielec = 1, elec_num do ielec = 1, elec_num
e_coord(1) = elec_coord(ielec,1,iwalk) e_coord(1) = elec_coord(ielec,1)
e_coord(2) = elec_coord(ielec,2,iwalk) e_coord(2) = elec_coord(ielec,2)
e_coord(3) = elec_coord(ielec,3,iwalk) e_coord(3) = elec_coord(ielec,3)
k=1 k=1
do inucl=1,nucl_num do inucl=1,nucl_num
n_coord(1) = nucl_coord(inucl,1) n_coord(1) = nucl_coord(inucl,1)
@ -3457,35 +3399,35 @@ integer function qmckl_compute_ao_vgl_f(context, &
l = shell_ang_mom(ishell) l = shell_ang_mom(ishell)
do il = lstart(l), lstart(l+1)-1 do il = lstart(l), lstart(l+1)-1
! Value ! Value
ao_vgl(k,ielec,iwalk,1) = & ao_vgl(k,ielec,1) = &
poly_vgl(1,il) * shell_vgl(ishell,ielec,iwalk,1) * ao_factor(k) poly_vgl(1,il) * shell_vgl(ishell,ielec,1) * ao_factor(k)
! Grad_x ! Grad_x
ao_vgl(k,ielec,iwalk,2) = ( & ao_vgl(k,ielec,2) = ( &
poly_vgl(2,il) * shell_vgl(ishell,ielec,iwalk,1) + & poly_vgl(2,il) * shell_vgl(ishell,ielec,1) + &
poly_vgl(1,il) * shell_vgl(ishell,ielec,iwalk,2) & poly_vgl(1,il) * shell_vgl(ishell,ielec,2) &
) * ao_factor(k) ) * ao_factor(k)
! Grad_y ! Grad_y
ao_vgl(k,ielec,iwalk,3) = ( & ao_vgl(k,ielec,3) = ( &
poly_vgl(3,il) * shell_vgl(ishell,ielec,iwalk,1) + & poly_vgl(3,il) * shell_vgl(ishell,ielec,1) + &
poly_vgl(1,il) * shell_vgl(ishell,ielec,iwalk,3) & poly_vgl(1,il) * shell_vgl(ishell,ielec,3) &
) * ao_factor(k) ) * ao_factor(k)
! Grad_z ! Grad_z
ao_vgl(k,ielec,iwalk,4) = ( & ao_vgl(k,ielec,4) = ( &
poly_vgl(4,il) * shell_vgl(ishell,ielec,iwalk,1) + & poly_vgl(4,il) * shell_vgl(ishell,ielec,1) + &
poly_vgl(1,il) * shell_vgl(ishell,ielec,iwalk,4) & poly_vgl(1,il) * shell_vgl(ishell,ielec,4) &
) * ao_factor(k) ) * ao_factor(k)
! Lapl_z ! Lapl_z
ao_vgl(k,ielec,iwalk,5) = ( & ao_vgl(k,ielec,5) = ( &
poly_vgl(5,il) * shell_vgl(ishell,ielec,iwalk,1) + & poly_vgl(5,il) * shell_vgl(ishell,ielec,1) + &
poly_vgl(1,il) * shell_vgl(ishell,ielec,iwalk,5) + & poly_vgl(1,il) * shell_vgl(ishell,ielec,5) + &
2.d0 * ( & 2.d0 * ( &
poly_vgl(2,il) * shell_vgl(ishell,ielec,iwalk,2) + & poly_vgl(2,il) * shell_vgl(ishell,ielec,2) + &
poly_vgl(3,il) * shell_vgl(ishell,ielec,iwalk,3) + & poly_vgl(3,il) * shell_vgl(ishell,ielec,3) + &
poly_vgl(4,il) * shell_vgl(ishell,ielec,iwalk,4) ) & poly_vgl(4,il) * shell_vgl(ishell,ielec,4) ) &
) * ao_factor(k) ) * ao_factor(k)
k = k+1 k = k+1
@ -3493,7 +3435,6 @@ integer function qmckl_compute_ao_vgl_f(context, &
end do end do
end do end do
end do end do
end do
deallocate(poly_vgl, powers) deallocate(poly_vgl, powers)
end function qmckl_compute_ao_vgl_f end function qmckl_compute_ao_vgl_f
@ -3509,7 +3450,6 @@ end function qmckl_compute_ao_vgl_f
const int64_t shell_num, const int64_t shell_num,
const int64_t elec_num, const int64_t elec_num,
const int64_t nucl_num, const int64_t nucl_num,
const int64_t walk_num,
const double* elec_coord, const double* elec_coord,
const double* nucl_coord, const double* nucl_coord,
const int64_t* nucleus_index, const int64_t* nucleus_index,
@ -3532,7 +3472,6 @@ end function qmckl_compute_ao_vgl_f
shell_num, & shell_num, &
elec_num, & elec_num, &
nucl_num, & nucl_num, &
walk_num, &
elec_coord, & elec_coord, &
nucl_coord, & nucl_coord, &
nucleus_index, & nucleus_index, &
@ -3553,8 +3492,7 @@ end function qmckl_compute_ao_vgl_f
integer (c_int64_t) , intent(in) , value :: shell_num integer (c_int64_t) , intent(in) , value :: shell_num
integer (c_int64_t) , intent(in) , value :: elec_num integer (c_int64_t) , intent(in) , value :: elec_num
integer (c_int64_t) , intent(in) , value :: nucl_num integer (c_int64_t) , intent(in) , value :: nucl_num
integer (c_int64_t) , intent(in) , value :: walk_num real (c_double ) , intent(in) :: elec_coord(elec_num,3)
real (c_double ) , intent(in) :: elec_coord(elec_num,3,walk_num)
real (c_double ) , intent(in) :: nucl_coord(nucl_num,3) real (c_double ) , intent(in) :: nucl_coord(nucl_num,3)
integer (c_int64_t) , intent(in) :: nucleus_index(nucl_num) integer (c_int64_t) , intent(in) :: nucleus_index(nucl_num)
integer (c_int64_t) , intent(in) :: nucleus_shell_num(nucl_num) integer (c_int64_t) , intent(in) :: nucleus_shell_num(nucl_num)
@ -3562,8 +3500,8 @@ end function qmckl_compute_ao_vgl_f
integer (c_int32_t) , intent(in) :: nucleus_max_ang_mom(nucl_num) integer (c_int32_t) , intent(in) :: nucleus_max_ang_mom(nucl_num)
integer (c_int32_t) , intent(in) :: shell_ang_mom(shell_num) integer (c_int32_t) , intent(in) :: shell_ang_mom(shell_num)
real (c_double ) , intent(in) :: ao_factor(ao_num) real (c_double ) , intent(in) :: ao_factor(ao_num)
real (c_double ) , intent(in) :: shell_vgl(shell_num,elec_num,walk_num,5) real (c_double ) , intent(in) :: shell_vgl(shell_num,elec_num,5)
real (c_double ) , intent(out) :: ao_vgl(ao_num,elec_num,walk_num,5) real (c_double ) , intent(out) :: ao_vgl(ao_num,elec_num,5)
integer(c_int32_t), external :: qmckl_compute_ao_vgl_f integer(c_int32_t), external :: qmckl_compute_ao_vgl_f
info = qmckl_compute_ao_vgl_f & info = qmckl_compute_ao_vgl_f &
@ -3572,7 +3510,6 @@ end function qmckl_compute_ao_vgl_f
shell_num, & shell_num, &
elec_num, & elec_num, &
nucl_num, & nucl_num, &
walk_num, &
elec_coord, & elec_coord, &
nucl_coord, & nucl_coord, &
nucleus_index, & nucleus_index, &
@ -3589,6 +3526,7 @@ end function qmckl_compute_ao_vgl_f
#+begin_src python :results output :exports none #+begin_src python :results output :exports none
import numpy as np import numpy as np
from math import sqrt
def f(a,x,y): def f(a,x,y):
return np.sum( [c * np.exp( -b*(np.linalg.norm(x-y))**2) for b,c in a] ) return np.sum( [c * np.exp( -b*(np.linalg.norm(x-y))**2) for b,c in a] )
@ -3614,7 +3552,7 @@ elec_26_w1 = np.array( [ 1.49050402641, 2.90106987953, -1.05920815468 ] )
elec_15_w2 = np.array( [ -2.20180344582,-1.9113150239, 2.2193744778600002 ] ) elec_15_w2 = np.array( [ -2.20180344582,-1.9113150239, 2.2193744778600002 ] )
nucl_1 = np.array( [ -2.302574592081335e+00, -3.542027060505035e-01, -5.334129934317614e-02] ) nucl_1 = np.array( [ -2.302574592081335e+00, -3.542027060505035e-01, -5.334129934317614e-02] )
#double prim_vgl[prim_num][5][walk_num][elec_num]; #double ao_vgl[prim_num][5][elec_num];
x = elec_26_w1 ; y = nucl_1 x = elec_26_w1 ; y = nucl_1
a = [( 403.830000, 0.001473 * 5.9876577632594533e+04), a = [( 403.830000, 0.001473 * 5.9876577632594533e+04),
( 121.170000, 0.012672 * 7.2836806319891484e+03), ( 121.170000, 0.012672 * 7.2836806319891484e+03),
@ -3625,40 +3563,41 @@ a = [( 403.830000, 0.001473 * 5.9876577632594533e+04),
( 1.763600, 0.273774 * 4.4423176930919421e+00), ( 1.763600, 0.273774 * 4.4423176930919421e+00),
( 0.706190, 0.074397 * 8.9541051939952665e-01)] ( 0.706190, 0.074397 * 8.9541051939952665e-01)]
print ( "[0][0][26][219] : %25.15e"%(f(a,x,y) * (x[0] - y[0])**2) ) norm = sqrt(3.)
print ( "[1][0][26][219] : %25.15e"%(df(a,x,y,1)* (x[0] - y[0]) * (x[1] - y[1]) + 2.*f(a,x,y) * (x[0] - y[0])) ) print ( "[0][26][219] : %25.15e"%(f(a,x,y) * (x[0] - y[0])**2) )
print ( "[1][26][219] : %25.15e"%(df(a,x,y,1)* (x[0] - y[0]) * (x[1] - y[1]) + 2.*f(a,x,y) * (x[0] - y[0])) )
print ( "[0][0][26][220] : %25.15e"%(f(a,x,y) * (x[0] - y[0]) * (x[1] - y[1])) ) print ( "[0][26][220] : %25.15e"%(norm*f(a,x,y) * (x[0] - y[0]) * (x[1] - y[1]) ))
print ( "[1][0][26][220] : %25.15e"%(df(a,x,y,1)* (x[0] - y[0]) * (x[1] - y[1]) + f(a,x,y) * (x[1] - y[1])) ) print ( "[1][26][220] : %25.15e"%(norm*df(a,x,y,1)* (x[0] - y[0]) * (x[1] - y[1]) + norm*f(a,x,y) * (x[1] - y[1])) )
print ( "[0][0][26][221] : %25.15e"%(f(a,x,y) * (x[0] - y[0]) * (x[2] - y[2])) ) print ( "[0][26][221] : %25.15e"%(norm*f(a,x,y) * (x[0] - y[0]) * (x[2] - y[2])) )
print ( "[1][0][26][221] : %25.15e"%(df(a,x,y,1)* (x[0] - y[0]) * (x[2] - y[2]) + f(a,x,y) * (x[2] - y[2])) ) print ( "[1][26][221] : %25.15e"%(norm*df(a,x,y,1)* (x[0] - y[0]) * (x[2] - y[2]) + norm*f(a,x,y) * (x[2] - y[2])) )
print ( "[0][0][26][222] : %25.15e"%(f(a,x,y) * (x[1] - y[1]) * (x[1] - y[1])) ) print ( "[0][26][222] : %25.15e"%(f(a,x,y) * (x[1] - y[1]) * (x[1] - y[1])) )
print ( "[1][0][26][222] : %25.15e"%(df(a,x,y,1)* (x[1] - y[1]) * (x[1] - y[1])) ) print ( "[1][26][222] : %25.15e"%(df(a,x,y,1)* (x[1] - y[1]) * (x[1] - y[1])) )
print ( "[0][0][26][223] : %25.15e"%(f(a,x,y) * (x[1] - y[1]) * (x[2] - y[2])) ) print ( "[0][26][223] : %25.15e"%(norm*f(a,x,y) * (x[1] - y[1]) * (x[2] - y[2])) )
print ( "[1][0][26][223] : %25.15e"%(df(a,x,y,1)* (x[1] - y[1]) * (x[2] - y[2])) ) print ( "[1][26][223] : %25.15e"%(norm*df(a,x,y,1)* (x[1] - y[1]) * (x[2] - y[2])) )
print ( "[0][0][26][224] : %25.15e"%(f(a,x,y) * (x[2] - y[2]) * (x[2] - y[2])) ) print ( "[0][26][224] : %25.15e"%(f(a,x,y) * (x[2] - y[2]) * (x[2] - y[2])) )
print ( "[1][0][26][224] : %25.15e"%(df(a,x,y,1)* (x[2] - y[2]) * (x[2] - y[2])) ) print ( "[1][26][224] : %25.15e"%(df(a,x,y,1)* (x[2] - y[2]) * (x[2] - y[2])) )
#+end_src #+end_src
#+RESULTS: #+RESULTS:
#+begin_example #+begin_example
[0][0][26][219] : 1.020302912653649e-08 [0][26][219] : 1.020302912653649e-08
[1][0][26][219] : -4.153046808203204e-08 [1][26][219] : -4.153046808203204e-08
[0][0][26][220] : 8.756380857379661e-09 [0][26][220] : 1.516649653540510e-08
[1][0][26][220] : -4.460176677299534e-08 [1][26][220] : -7.725252615816528e-08
[0][0][26][221] : -2.705688401075445e-09 [0][26][221] : -4.686389780112468e-09
[1][0][26][221] : 1.378177639720419e-08 [1][26][221] : 2.387073693851122e-08
[0][0][26][222] : 7.514847283937212e-09 [0][26][222] : 7.514847283937212e-09
[1][0][26][222] : -4.025905373647693e-08 [1][26][222] : -4.025905373647693e-08
[0][0][26][223] : -2.322059246071533e-09 [0][26][223] : -4.021924592380977e-09
[1][0][26][223] : 1.243989457599443e-08 [1][26][223] : 2.154652944642284e-08
[0][0][26][224] : 7.175074806631758e-10 [0][26][224] : 7.175074806631758e-10
[1][0][26][224] : -3.843880138733679e-09 [1][26][224] : -3.843880138733679e-09
#+end_example #+end_example
*** Test *** Test
@ -3686,38 +3625,38 @@ rc = qmckl_set_electron_coord (context, 'N', elec_coord);
assert(rc == QMCKL_SUCCESS); assert(rc == QMCKL_SUCCESS);
double ao_vgl[5][walk_num][elec_num][ao_num]; double ao_vgl[5][elec_num][ao_num];
rc = qmckl_get_ao_vgl(context, &(ao_vgl[0][0][0][0])); rc = qmckl_get_ao_vgl(context, &(ao_vgl[0][0][0]));
assert (rc == QMCKL_SUCCESS); assert (rc == QMCKL_SUCCESS);
printf("\n"); printf("\n");
printf(" ao_vgl ao_vgl[0][0][26][219] %25.15e\n", ao_vgl[0][0][26][219]); printf(" ao_vgl ao_vgl[0][26][219] %25.15e\n", ao_vgl[0][26][219]);
printf(" ao_vgl ao_vgl[1][0][26][219] %25.15e\n", ao_vgl[1][0][26][219]); printf(" ao_vgl ao_vgl[1][26][219] %25.15e\n", ao_vgl[1][26][219]);
printf(" ao_vgl ao_vgl[0][0][26][220] %25.15e\n", ao_vgl[0][0][26][220]); printf(" ao_vgl ao_vgl[0][26][220] %25.15e\n", ao_vgl[0][26][220]);
printf(" ao_vgl ao_vgl[1][0][26][220] %25.15e\n", ao_vgl[1][0][26][220]); printf(" ao_vgl ao_vgl[1][26][220] %25.15e\n", ao_vgl[1][26][220]);
printf(" ao_vgl ao_vgl[0][0][26][221] %25.15e\n", ao_vgl[0][0][26][221]); printf(" ao_vgl ao_vgl[0][26][221] %25.15e\n", ao_vgl[0][26][221]);
printf(" ao_vgl ao_vgl[1][0][26][221] %25.15e\n", ao_vgl[1][0][26][221]); printf(" ao_vgl ao_vgl[1][26][221] %25.15e\n", ao_vgl[1][26][221]);
printf(" ao_vgl ao_vgl[0][0][26][222] %25.15e\n", ao_vgl[0][0][26][222]); printf(" ao_vgl ao_vgl[0][26][222] %25.15e\n", ao_vgl[0][26][222]);
printf(" ao_vgl ao_vgl[1][0][26][222] %25.15e\n", ao_vgl[1][0][26][222]); printf(" ao_vgl ao_vgl[1][26][222] %25.15e\n", ao_vgl[1][26][222]);
printf(" ao_vgl ao_vgl[0][0][26][223] %25.15e\n", ao_vgl[0][0][26][223]); printf(" ao_vgl ao_vgl[0][26][223] %25.15e\n", ao_vgl[0][26][223]);
printf(" ao_vgl ao_vgl[1][0][26][223] %25.15e\n", ao_vgl[1][0][26][223]); printf(" ao_vgl ao_vgl[1][26][223] %25.15e\n", ao_vgl[1][26][223]);
printf(" ao_vgl ao_vgl[0][0][26][224] %25.15e\n", ao_vgl[0][0][26][224]); printf(" ao_vgl ao_vgl[0][26][224] %25.15e\n", ao_vgl[0][26][224]);
printf(" ao_vgl ao_vgl[1][0][26][224] %25.15e\n", ao_vgl[1][0][26][224]); printf(" ao_vgl ao_vgl[1][26][224] %25.15e\n", ao_vgl[1][26][224]);
printf("\n"); printf("\n");
assert( fabs(ao_vgl[0][0][26][219] - ( 1.020298798341620e-08)) < 1.e-14 ); assert( fabs(ao_vgl[0][26][219] - ( 1.020298798341620e-08)) < 1.e-14 );
assert( fabs(ao_vgl[1][0][26][219] - ( -4.928035238010602e-08)) < 1.e-14 ); assert( fabs(ao_vgl[1][26][219] - (-4.928035238010602e-08)) < 1.e-14 );
assert( fabs(ao_vgl[0][0][26][220] - ( 8.756345547784206e-09)) < 1.e-14 ); assert( fabs(ao_vgl[0][26][220] - ( 1.516643537739178e-08)) < 1.e-14 );
assert( fabs(ao_vgl[1][0][26][220] - ( -4.460158690983819e-08)) < 1.e-14 ); assert( fabs(ao_vgl[1][26][220] - (-7.725221462603871e-08)) < 1.e-14 );
assert( fabs(ao_vgl[0][0][26][221] - ( -2.705677490544664e-09)) < 1.e-14 ); assert( fabs(ao_vgl[0][26][221] - (-4.686370882518819e-09)) < 1.e-14 );
assert( fabs(ao_vgl[1][0][26][221] - ( 1.378172082017231e-08)) < 1.e-14 ); assert( fabs(ao_vgl[1][26][221] - ( 2.387064067626827e-08)) < 1.e-14 );
assert( fabs(ao_vgl[0][0][26][222] - ( 7.514816980753531e-09)) < 1.e-14 ); assert( fabs(ao_vgl[0][26][222] - ( 7.514816980753531e-09)) < 1.e-14 );
assert( fabs(ao_vgl[1][0][26][222] - ( -4.025889138635182e-08)) < 1.e-14 ); assert( fabs(ao_vgl[1][26][222] - (-4.025889138635182e-08)) < 1.e-14 );
assert( fabs(ao_vgl[0][0][26][223] - ( -2.322049882502961e-09)) < 1.e-14 ); assert( fabs(ao_vgl[0][26][223] - (-4.021908374204471e-09)) < 1.e-14 );
assert( fabs(ao_vgl[1][0][26][223] - ( 1.243984441042288e-08)) < 1.e-14 ); assert( fabs(ao_vgl[1][26][223] - ( 2.154644255710413e-08)) < 1.e-14 );
assert( fabs(ao_vgl[0][0][26][224] - ( 7.175045873560788e-10)) < 1.e-14 ); assert( fabs(ao_vgl[0][26][224] - ( 7.175045873560788e-10)) < 1.e-14 );
assert( fabs(ao_vgl[1][0][26][224] - ( -3.843864637762753e-09)) < 1.e-14 ); assert( fabs(ao_vgl[1][26][224] - (-3.843864637762753e-09)) < 1.e-14 );
} }
#+end_src #+end_src

File diff suppressed because it is too large Load Diff