mirror of
https://github.com/TREX-CoE/qmckl.git
synced 2024-08-24 22:21:46 +02:00
Merge branch 'master' into master
This commit is contained in:
commit
6d19d416fa
477
org/qmckl_ao.org
477
org/qmckl_ao.org
@ -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")
|
||||||
@ -116,26 +119,29 @@ 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 |
|
||||||
| ~prim_factor_sorted~ | ~[prim_num]~ | Normalization factors of the sorted primtives |
|
| ~prim_factor_sorted~ | ~[prim_num]~ | Normalization factors of the sorted primtives |
|
||||||
|
|
||||||
|
|
||||||
For H_2 with the following basis set,
|
For H_2 with the following basis set,
|
||||||
|
|
||||||
|
#+NAME: basis
|
||||||
#+BEGIN_EXAMPLE
|
#+BEGIN_EXAMPLE
|
||||||
HYDROGEN
|
HYDROGEN
|
||||||
S 5
|
S 5
|
||||||
@ -158,6 +164,7 @@ D 1
|
|||||||
|
|
||||||
we have:
|
we have:
|
||||||
|
|
||||||
|
#+NAME: params
|
||||||
#+BEGIN_EXAMPLE
|
#+BEGIN_EXAMPLE
|
||||||
type = 'G'
|
type = 'G'
|
||||||
shell_num = 12
|
shell_num = 12
|
||||||
@ -1865,7 +1872,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;
|
||||||
@ -1904,7 +1911,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) {
|
||||||
@ -1922,7 +1929,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,
|
||||||
@ -1953,20 +1959,20 @@ qmckl_exit_code qmckl_provide_ao_basis_primitive_vgl(qmckl_context context)
|
|||||||
:END:
|
:END:
|
||||||
|
|
||||||
#+NAME: qmckl_ao_basis_primitive_gaussian_vgl_args
|
#+NAME: qmckl_ao_basis_primitive_gaussian_vgl_args
|
||||||
|
|
||||||
| qmckl_context | context | in | Global state |
|
| qmckl_context | context | in | Global state |
|
||||||
| 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
|
||||||
@ -1975,14 +1981,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
|
||||||
@ -1993,11 +1998,11 @@ 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
|
||||||
@ -2006,16 +2011,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
|
||||||
@ -2026,7 +2030,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,
|
||||||
@ -2039,7 +2042,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
|
||||||
@ -2049,16 +2060,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
|
||||||
@ -2087,40 +2105,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
|
||||||
|
|
||||||
@ -2147,22 +2150,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 );
|
|
||||||
|
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -2172,13 +2169,12 @@ 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];
|
||||||
@ -2195,7 +2191,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
|
||||||
@ -2225,7 +2220,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;
|
||||||
@ -2284,8 +2279,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) {
|
||||||
@ -2304,7 +2298,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,
|
||||||
@ -2339,25 +2332,25 @@ qmckl_exit_code qmckl_provide_ao_basis_shell_vgl(qmckl_context context)
|
|||||||
:END:
|
:END:
|
||||||
|
|
||||||
#+NAME: qmckl_ao_basis_shell_gaussian_vgl_args
|
#+NAME: qmckl_ao_basis_shell_gaussian_vgl_args
|
||||||
|
|
||||||
| ~qmckl_context~ | ~context~ | in | Global state |
|
| ~qmckl_context~ | ~context~ | in | Global state |
|
||||||
| ~int64_t~ | ~prim_num~ | in | Number of primitives |
|
| ~int64_t~ | ~prim_num~ | in | Number of primitives |
|
||||||
| ~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)
|
||||||
@ -2368,18 +2361,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
|
||||||
@ -2390,23 +2382,23 @@ 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)
|
||||||
|
|
||||||
@ -2418,27 +2410,28 @@ 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, iwalk, 1) + v
|
|
||||||
|
|
||||||
shell_vgl(ishell, ielec, iwalk, 2) = &
|
shell_vgl(ishell, ielec, 1) = &
|
||||||
shell_vgl(ishell, ielec, iwalk, 2) + two_a * x
|
shell_vgl(ishell, ielec, 1) + v
|
||||||
|
|
||||||
shell_vgl(ishell, ielec, iwalk, 3) = &
|
shell_vgl(ishell, ielec, 2) = &
|
||||||
shell_vgl(ishell, ielec, iwalk, 3) + two_a * y
|
shell_vgl(ishell, ielec, 2) + two_a * x
|
||||||
|
|
||||||
shell_vgl(ishell, ielec, iwalk, 4) = &
|
shell_vgl(ishell, ielec, 3) = &
|
||||||
shell_vgl(ishell, ielec, iwalk, 4) + two_a * z
|
shell_vgl(ishell, ielec, 3) + two_a * y
|
||||||
|
|
||||||
shell_vgl(ishell, ielec, iwalk, 5) = &
|
shell_vgl(ishell, ielec, 4) = &
|
||||||
shell_vgl(ishell, ielec, iwalk, 5) + two_a * (3.d0 - 2.d0*ar2)
|
shell_vgl(ishell, ielec, 4) + two_a * z
|
||||||
|
|
||||||
|
shell_vgl(ishell, ielec, 5) = &
|
||||||
|
shell_vgl(ishell, ielec, 5) + two_a * (3.d0 - 2.d0*ar2)
|
||||||
|
|
||||||
end do
|
end do
|
||||||
|
|
||||||
end do
|
end do
|
||||||
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
|
||||||
@ -2454,7 +2447,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,
|
||||||
@ -2476,7 +2468,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, &
|
||||||
@ -2496,16 +2487,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 &
|
||||||
@ -2514,7 +2504,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, &
|
||||||
@ -2556,7 +2545,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 ),
|
||||||
@ -2569,40 +2558,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
|
||||||
|
|
||||||
@ -2628,34 +2597,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
|
||||||
@ -3309,7 +3267,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;
|
||||||
@ -3376,8 +3334,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) {
|
||||||
@ -3394,7 +3351,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,
|
||||||
@ -3424,13 +3380,13 @@ qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context)
|
|||||||
:END:
|
:END:
|
||||||
|
|
||||||
#+NAME: qmckl_ao_vgl_args
|
#+NAME: qmckl_ao_vgl_args
|
||||||
|
|
||||||
| ~qmckl_context~ | ~context~ | in | Global state |
|
| ~qmckl_context~ | ~context~ | in | Global state |
|
||||||
| ~int64_t~ | ~ao_num~ | in | Number of AOs |
|
| ~int64_t~ | ~ao_num~ | in | Number of AOs |
|
||||||
| ~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 |
|
||||||
@ -3438,12 +3394,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) &
|
||||||
@ -3455,8 +3411,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)
|
||||||
@ -3464,13 +3419,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
|
||||||
@ -3492,11 +3447,11 @@ 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)
|
||||||
@ -3524,35 +3479,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
|
||||||
@ -3560,7 +3515,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
|
||||||
@ -3576,7 +3530,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,
|
||||||
@ -3599,7 +3552,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, &
|
||||||
@ -3620,8 +3572,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)
|
||||||
@ -3629,8 +3580,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 &
|
||||||
@ -3639,7 +3590,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, &
|
||||||
@ -3656,6 +3606,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] )
|
||||||
@ -3681,7 +3632,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),
|
||||||
@ -3692,40 +3643,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
|
||||||
@ -3753,38 +3705,39 @@ 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[1][0][26][219] - ( -4.928035238010602e-08)) < 1.e-14 );
|
assert( fabs(ao_vgl[0][26][219] - ( 1.020298798341620e-08)) < 1.e-14 );
|
||||||
assert( fabs(ao_vgl[0][0][26][220] - ( 8.756345547784206e-09)) < 1.e-14 );
|
assert( fabs(ao_vgl[1][26][219] - (-4.928035238010602e-08)) < 1.e-14 );
|
||||||
assert( fabs(ao_vgl[1][0][26][220] - ( -4.460158690983819e-08)) < 1.e-14 );
|
assert( fabs(ao_vgl[0][26][220] - ( 1.516643537739178e-08)) < 1.e-14 );
|
||||||
assert( fabs(ao_vgl[0][0][26][221] - ( -2.705677490544664e-09)) < 1.e-14 );
|
assert( fabs(ao_vgl[1][26][220] - (-7.725221462603871e-08)) < 1.e-14 );
|
||||||
assert( fabs(ao_vgl[1][0][26][221] - ( 1.378172082017231e-08)) < 1.e-14 );
|
assert( fabs(ao_vgl[0][26][221] - (-4.686370882518819e-09)) < 1.e-14 );
|
||||||
assert( fabs(ao_vgl[0][0][26][222] - ( 7.514816980753531e-09)) < 1.e-14 );
|
assert( fabs(ao_vgl[1][26][221] - ( 2.387064067626827e-08)) < 1.e-14 );
|
||||||
assert( fabs(ao_vgl[1][0][26][222] - ( -4.025889138635182e-08)) < 1.e-14 );
|
assert( fabs(ao_vgl[0][26][222] - ( 7.514816980753531e-09)) < 1.e-14 );
|
||||||
assert( fabs(ao_vgl[0][0][26][223] - ( -2.322049882502961e-09)) < 1.e-14 );
|
assert( fabs(ao_vgl[1][26][222] - (-4.025889138635182e-08)) < 1.e-14 );
|
||||||
assert( fabs(ao_vgl[1][0][26][223] - ( 1.243984441042288e-08)) < 1.e-14 );
|
assert( fabs(ao_vgl[0][26][223] - (-4.021908374204471e-09)) < 1.e-14 );
|
||||||
assert( fabs(ao_vgl[0][0][26][224] - ( 7.175045873560788e-10)) < 1.e-14 );
|
assert( fabs(ao_vgl[1][26][223] - ( 2.154644255710413e-08)) < 1.e-14 );
|
||||||
assert( fabs(ao_vgl[1][0][26][224] - ( -3.843864637762753e-09)) < 1.e-14 );
|
assert( fabs(ao_vgl[0][26][224] - ( 7.175045873560788e-10)) < 1.e-14 );
|
||||||
|
assert( fabs(ao_vgl[1][26][224] - (-3.843864637762753e-09)) < 1.e-14 );
|
||||||
}
|
}
|
||||||
|
|
||||||
#+end_src
|
#+end_src
|
||||||
|
322
org/qmckl_blas.org
Normal file
322
org/qmckl_blas.org
Normal file
@ -0,0 +1,322 @@
|
|||||||
|
#+TITLE: BLAS functions
|
||||||
|
#+SETUPFILE: ../tools/theme.setup
|
||||||
|
#+INCLUDE: ../tools/lib.org
|
||||||
|
|
||||||
|
* Headers :noexport:
|
||||||
|
#+begin_src elisp :noexport :results none
|
||||||
|
(org-babel-lob-ingest "../tools/lib.org")
|
||||||
|
#+end_src
|
||||||
|
|
||||||
|
#+begin_src c :comments link :tangle (eval c_test) :noweb yes
|
||||||
|
#include "qmckl.h"
|
||||||
|
#include "assert.h"
|
||||||
|
#ifdef HAVE_CONFIG_H
|
||||||
|
#include "config.h"
|
||||||
|
#endif
|
||||||
|
int main() {
|
||||||
|
qmckl_context context;
|
||||||
|
context = qmckl_context_create();
|
||||||
|
|
||||||
|
#+end_src
|
||||||
|
|
||||||
|
* Matrix operations
|
||||||
|
|
||||||
|
** ~qmckl_dgemm~
|
||||||
|
|
||||||
|
Matrix multiply l$C_{ij} = \beta C_{ij} + \alpha \sum_{k} A_{ik} \cdot B_{kj}$ using Fortran ~matmul~ function.
|
||||||
|
TODO: Add description about the external library dependence.
|
||||||
|
|
||||||
|
#+NAME: qmckl_dgemm_args
|
||||||
|
| qmckl_context | context | in | Global state |
|
||||||
|
| bool | TransA | in | Number of rows of the input matrix |
|
||||||
|
| bool | TransB | in | Number of rows of the input matrix |
|
||||||
|
| int64_t | m | in | Number of rows of the input matrix |
|
||||||
|
| int64_t | n | in | Number of columns of the input matrix |
|
||||||
|
| int64_t | k | in | Number of columns of the input matrix |
|
||||||
|
| double | alpha | in | Number of columns of the input matrix |
|
||||||
|
| double | A[][lda] | in | Array containing the $m \times n$ matrix $A$ |
|
||||||
|
| int64_t | lda | in | Leading dimension of array ~A~ |
|
||||||
|
| double | B[][ldb] | in | Array containing the $n \times m$ matrix $B$ |
|
||||||
|
| int64_t | ldb | in | Leading dimension of array ~B~ |
|
||||||
|
| double | beta | in | Array containing the $n \times m$ matrix $B$ |
|
||||||
|
| double | C[][ldc] | out | Array containing the $n \times m$ matrix $B$ |
|
||||||
|
| int64_t | ldc | in | Leading dimension of array ~B~ |
|
||||||
|
|
||||||
|
*** Requirements
|
||||||
|
|
||||||
|
- ~context~ is not ~QMCKL_NULL_CONTEXT~
|
||||||
|
- ~m > 0~
|
||||||
|
- ~n > 0~
|
||||||
|
- ~k > 0~
|
||||||
|
- ~lda >= m~
|
||||||
|
- ~ldb >= n~
|
||||||
|
- ~ldc >= n~
|
||||||
|
- ~A~ is allocated with at least $m \times k \times 8$ bytes
|
||||||
|
- ~B~ is allocated with at least $k \times n \times 8$ bytes
|
||||||
|
- ~C~ is allocated with at least $m \times n \times 8$ bytes
|
||||||
|
|
||||||
|
*** C header
|
||||||
|
|
||||||
|
#+CALL: generate_c_header(table=qmckl_dgemm_args,rettyp="qmckl_exit_code",fname="qmckl_dgemm")
|
||||||
|
|
||||||
|
#+RESULTS:
|
||||||
|
#+BEGIN_src c :tangle (eval h_func) :comments org
|
||||||
|
qmckl_exit_code qmckl_dgemm (
|
||||||
|
const qmckl_context context,
|
||||||
|
const bool TransA,
|
||||||
|
const bool TransB,
|
||||||
|
const int64_t m,
|
||||||
|
const int64_t n,
|
||||||
|
const int64_t k,
|
||||||
|
const double alpha,
|
||||||
|
const double* A,
|
||||||
|
const int64_t lda,
|
||||||
|
const double* B,
|
||||||
|
const int64_t ldb,
|
||||||
|
const double beta,
|
||||||
|
double* const C,
|
||||||
|
const int64_t ldc );
|
||||||
|
#+END_src
|
||||||
|
|
||||||
|
|
||||||
|
*** Source
|
||||||
|
#+begin_src f90 :tangle (eval f)
|
||||||
|
integer function qmckl_dgemm_f(context, TransA, TransB, m, n, k, alpha, A, LDA, B, LDB, beta, C, LDC) &
|
||||||
|
result(info)
|
||||||
|
use qmckl
|
||||||
|
implicit none
|
||||||
|
integer(qmckl_context) , intent(in) :: context
|
||||||
|
logical*8 , intent(in) :: TransA, TransB
|
||||||
|
integer*8 , intent(in) :: m, n, k
|
||||||
|
real*8 , intent(in) :: alpha, beta
|
||||||
|
integer*8 , intent(in) :: lda
|
||||||
|
real*8 , intent(in) :: A(m,k)
|
||||||
|
integer*8 , intent(in) :: ldb
|
||||||
|
real*8 , intent(in) :: B(k,n)
|
||||||
|
integer*8 , intent(in) :: ldc
|
||||||
|
real*8 , intent(out) :: C(m,n)
|
||||||
|
real*8, allocatable :: AT(:,:), BT(:,:), CT(:,:)
|
||||||
|
|
||||||
|
integer*8 :: i,j,l, LDA_2, LDB_2
|
||||||
|
|
||||||
|
info = QMCKL_SUCCESS
|
||||||
|
|
||||||
|
if (TransA) then
|
||||||
|
allocate(AT(k,m))
|
||||||
|
do i = 1, m
|
||||||
|
do j = 1, k
|
||||||
|
AT(j,i) = A(i,j)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
LDA_2 = M
|
||||||
|
else
|
||||||
|
LDA_2 = LDA
|
||||||
|
endif
|
||||||
|
|
||||||
|
if (TransB) then
|
||||||
|
allocate(BT(n,k))
|
||||||
|
do i = 1, k
|
||||||
|
do j = 1, n
|
||||||
|
BT(j,i) = B(i,j)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
LDB_2 = K
|
||||||
|
else
|
||||||
|
LDB_2 = LDB
|
||||||
|
endif
|
||||||
|
|
||||||
|
if (context == QMCKL_NULL_CONTEXT) then
|
||||||
|
info = QMCKL_INVALID_CONTEXT
|
||||||
|
return
|
||||||
|
endif
|
||||||
|
|
||||||
|
if (m <= 0_8) then
|
||||||
|
info = QMCKL_INVALID_ARG_4
|
||||||
|
return
|
||||||
|
endif
|
||||||
|
|
||||||
|
if (n <= 0_8) then
|
||||||
|
info = QMCKL_INVALID_ARG_5
|
||||||
|
return
|
||||||
|
endif
|
||||||
|
|
||||||
|
if (k <= 0_8) then
|
||||||
|
info = QMCKL_INVALID_ARG_6
|
||||||
|
return
|
||||||
|
endif
|
||||||
|
|
||||||
|
if (LDA_2 .ne. m) then
|
||||||
|
info = QMCKL_INVALID_ARG_9
|
||||||
|
return
|
||||||
|
endif
|
||||||
|
|
||||||
|
if (LDB_2 .ne. k) then
|
||||||
|
info = QMCKL_INVALID_ARG_10
|
||||||
|
return
|
||||||
|
endif
|
||||||
|
|
||||||
|
if (LDC .ne. m) then
|
||||||
|
info = QMCKL_INVALID_ARG_13
|
||||||
|
return
|
||||||
|
endif
|
||||||
|
|
||||||
|
if (TransA) then
|
||||||
|
C = matmul(AT,B)
|
||||||
|
else if (TransB) then
|
||||||
|
C = matmul(A,BT)
|
||||||
|
else
|
||||||
|
C = matmul(A,B)
|
||||||
|
endif
|
||||||
|
end function qmckl_dgemm_f
|
||||||
|
#+end_src
|
||||||
|
|
||||||
|
*** C interface :noexport:
|
||||||
|
|
||||||
|
#+CALL: generate_c_interface(table=qmckl_dgemm_args,rettyp="qmckl_exit_code",fname="qmckl_dgemm")
|
||||||
|
|
||||||
|
#+RESULTS:
|
||||||
|
#+BEGIN_src f90 :tangle (eval f) :comments org :exports none
|
||||||
|
integer(c_int32_t) function qmckl_dgemm &
|
||||||
|
(context, TransA, TransB, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc) &
|
||||||
|
bind(C) result(info)
|
||||||
|
|
||||||
|
use, intrinsic :: iso_c_binding
|
||||||
|
implicit none
|
||||||
|
|
||||||
|
integer (c_int64_t) , intent(in) , value :: context
|
||||||
|
logical*8 , intent(in) , value :: TransA
|
||||||
|
logical*8 , intent(in) , value :: TransB
|
||||||
|
integer (c_int64_t) , intent(in) , value :: m
|
||||||
|
integer (c_int64_t) , intent(in) , value :: n
|
||||||
|
integer (c_int64_t) , intent(in) , value :: k
|
||||||
|
real (c_double ) , intent(in) , value :: alpha
|
||||||
|
integer (c_int64_t) , intent(in) , value :: lda
|
||||||
|
real (c_double ) , intent(in) :: A(lda,*)
|
||||||
|
integer (c_int64_t) , intent(in) , value :: ldb
|
||||||
|
real (c_double ) , intent(in) :: B(ldb,*)
|
||||||
|
real (c_double ) , intent(in) , value :: beta
|
||||||
|
integer (c_int64_t) , intent(in) , value :: ldc
|
||||||
|
real (c_double ) , intent(out) :: C(ldc,*)
|
||||||
|
|
||||||
|
integer(c_int32_t), external :: qmckl_dgemm_f
|
||||||
|
info = qmckl_dgemm_f &
|
||||||
|
(context, TransA, TransB, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc)
|
||||||
|
|
||||||
|
end function qmckl_dgemm
|
||||||
|
#+END_src
|
||||||
|
|
||||||
|
|
||||||
|
#+CALL: generate_f_interface(table=qmckl_dgemm_args,rettyp="qmckl_exit_code",fname="qmckl_dgemm")
|
||||||
|
|
||||||
|
#+RESULTS:
|
||||||
|
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
|
||||||
|
interface
|
||||||
|
integer(c_int32_t) function qmckl_dgemm &
|
||||||
|
(context, TransA, TransB, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc) &
|
||||||
|
bind(C)
|
||||||
|
use, intrinsic :: iso_c_binding
|
||||||
|
import
|
||||||
|
implicit none
|
||||||
|
|
||||||
|
integer (c_int64_t) , intent(in) , value :: context
|
||||||
|
logical*8 , intent(in) , value :: TransA
|
||||||
|
logical*8 , intent(in) , value :: TransB
|
||||||
|
integer (c_int64_t) , intent(in) , value :: m
|
||||||
|
integer (c_int64_t) , intent(in) , value :: n
|
||||||
|
integer (c_int64_t) , intent(in) , value :: k
|
||||||
|
real (c_double ) , intent(in) , value :: alpha
|
||||||
|
integer (c_int64_t) , intent(in) , value :: lda
|
||||||
|
real (c_double ) , intent(in) :: A(lda,*)
|
||||||
|
integer (c_int64_t) , intent(in) , value :: ldb
|
||||||
|
real (c_double ) , intent(in) :: B(ldb,*)
|
||||||
|
real (c_double ) , intent(in) , value :: beta
|
||||||
|
integer (c_int64_t) , intent(in) , value :: ldc
|
||||||
|
real (c_double ) , intent(out) :: C(ldc,*)
|
||||||
|
|
||||||
|
end function qmckl_dgemm
|
||||||
|
end interface
|
||||||
|
#+end_src
|
||||||
|
|
||||||
|
*** Test :noexport:
|
||||||
|
#+begin_src f90 :tangle (eval f_test)
|
||||||
|
integer(qmckl_exit_code) function test_qmckl_dgemm(context) bind(C)
|
||||||
|
use qmckl
|
||||||
|
implicit none
|
||||||
|
integer(qmckl_context), intent(in), value :: context
|
||||||
|
|
||||||
|
double precision, allocatable :: A(:,:), B(:,:), C(:,:), D(:,:)
|
||||||
|
integer*8 :: m, n, k, LDA, LDB, LDC
|
||||||
|
integer*8 :: i,j,l
|
||||||
|
logical*8 :: TransA, TransB
|
||||||
|
double precision :: x, alpha, beta
|
||||||
|
|
||||||
|
TransA = .False.
|
||||||
|
TransB = .False.
|
||||||
|
m = 1_8
|
||||||
|
k = 4_8
|
||||||
|
n = 6_8
|
||||||
|
LDA = m
|
||||||
|
LDB = k
|
||||||
|
LDC = m
|
||||||
|
|
||||||
|
allocate( A(LDA,k), B(LDB,n) , C(LDC,n), D(LDC,n))
|
||||||
|
|
||||||
|
A = 0.d0
|
||||||
|
B = 0.d0
|
||||||
|
C = 0.d0
|
||||||
|
D = 0.d0
|
||||||
|
alpha = 1.0d0
|
||||||
|
beta = 0.0d0
|
||||||
|
do j=1,k
|
||||||
|
do i=1,m
|
||||||
|
A(i,j) = -10.d0 + dble(i+j)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
|
||||||
|
do j=1,n
|
||||||
|
do i=1,k
|
||||||
|
B(i,j) = -10.d0 + dble(i+j)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
|
||||||
|
test_qmckl_dgemm = qmckl_dgemm(context, TransA, TransB, m, n, k, alpha, A, LDA, B, LDB, beta, C, LDC)
|
||||||
|
|
||||||
|
if (test_qmckl_dgemm /= QMCKL_SUCCESS) return
|
||||||
|
|
||||||
|
test_qmckl_dgemm = QMCKL_FAILURE
|
||||||
|
|
||||||
|
x = 0.d0
|
||||||
|
do j=1,n
|
||||||
|
do i=1,m
|
||||||
|
do l=1,k
|
||||||
|
D(i,j) = D(i,j) + A(i,l)*B(l,j)
|
||||||
|
end do
|
||||||
|
x = x + (D(i,j) - C(i,j))**2
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
|
||||||
|
if (dabs(x) <= 1.d-15) then
|
||||||
|
test_qmckl_dgemm = QMCKL_SUCCESS
|
||||||
|
endif
|
||||||
|
|
||||||
|
deallocate(A,B,C,D)
|
||||||
|
end function test_qmckl_dgemm
|
||||||
|
#+end_src
|
||||||
|
|
||||||
|
#+begin_src c :comments link :tangle (eval c_test)
|
||||||
|
qmckl_exit_code test_qmckl_dgemm(qmckl_context context);
|
||||||
|
assert(QMCKL_SUCCESS == test_qmckl_dgemm(context));
|
||||||
|
#+end_src
|
||||||
|
|
||||||
|
* End of files :noexport:
|
||||||
|
|
||||||
|
#+begin_src c :comments link :tangle (eval c_test)
|
||||||
|
assert (qmckl_context_destroy(context) == QMCKL_SUCCESS);
|
||||||
|
return 0;
|
||||||
|
}
|
||||||
|
|
||||||
|
#+end_src
|
||||||
|
|
||||||
|
|
||||||
|
# -*- mode: org -*-
|
||||||
|
# vim: syntax=c
|
@ -31,10 +31,12 @@ int main() {
|
|||||||
#include "qmckl_nucleus_private_type.h"
|
#include "qmckl_nucleus_private_type.h"
|
||||||
#include "qmckl_electron_private_type.h"
|
#include "qmckl_electron_private_type.h"
|
||||||
#include "qmckl_ao_private_type.h"
|
#include "qmckl_ao_private_type.h"
|
||||||
|
#include "qmckl_mo_private_type.h"
|
||||||
#include "qmckl_jastrow_private_type.h"
|
#include "qmckl_jastrow_private_type.h"
|
||||||
#include "qmckl_nucleus_private_func.h"
|
#include "qmckl_nucleus_private_func.h"
|
||||||
#include "qmckl_electron_private_func.h"
|
#include "qmckl_electron_private_func.h"
|
||||||
#include "qmckl_ao_private_func.h"
|
#include "qmckl_ao_private_func.h"
|
||||||
|
#include "qmckl_mo_private_func.h"
|
||||||
#+end_src
|
#+end_src
|
||||||
|
|
||||||
#+begin_src c :tangle (eval c)
|
#+begin_src c :tangle (eval c)
|
||||||
@ -119,6 +121,7 @@ typedef struct qmckl_context_struct {
|
|||||||
qmckl_nucleus_struct nucleus;
|
qmckl_nucleus_struct nucleus;
|
||||||
qmckl_electron_struct electron;
|
qmckl_electron_struct electron;
|
||||||
qmckl_ao_basis_struct ao_basis;
|
qmckl_ao_basis_struct ao_basis;
|
||||||
|
qmckl_mo_basis_struct mo_basis;
|
||||||
qmckl_jastrow_struct jastrow;
|
qmckl_jastrow_struct jastrow;
|
||||||
|
|
||||||
/* To be implemented:
|
/* To be implemented:
|
||||||
|
@ -532,11 +532,11 @@ qmckl_set_electron_num(qmckl_context context,
|
|||||||
"up_num <= 0");
|
"up_num <= 0");
|
||||||
}
|
}
|
||||||
|
|
||||||
if (down_num <= 0) {
|
if (down_num < 0) {
|
||||||
return qmckl_failwith( context,
|
return qmckl_failwith( context,
|
||||||
QMCKL_INVALID_ARG_3,
|
QMCKL_INVALID_ARG_3,
|
||||||
"qmckl_set_electron_num",
|
"qmckl_set_electron_num",
|
||||||
"down_num <= 0");
|
"down_num < 0");
|
||||||
}
|
}
|
||||||
|
|
||||||
int32_t mask = 1 << 0;
|
int32_t mask = 1 << 0;
|
||||||
|
@ -202,7 +202,7 @@ for i in range(elec_num):
|
|||||||
type_nucl_num = 1
|
type_nucl_num = 1
|
||||||
aord_num = 5
|
aord_num = 5
|
||||||
bord_num = 5
|
bord_num = 5
|
||||||
cord_num = 23
|
cord_num = 5
|
||||||
dim_cord_vect= 23
|
dim_cord_vect= 23
|
||||||
type_nucl_vector = [ 1, 1]
|
type_nucl_vector = [ 1, 1]
|
||||||
aord_vector = [
|
aord_vector = [
|
||||||
@ -609,7 +609,7 @@ qmckl_exit_code qmckl_get_jastrow_cord_vector (const qmckl_context context, doub
|
|||||||
}
|
}
|
||||||
|
|
||||||
assert (ctx->jastrow.cord_vector != NULL);
|
assert (ctx->jastrow.cord_vector != NULL);
|
||||||
memcpy(cord_vector, ctx->jastrow.cord_vector, ctx->jastrow.cord_num*sizeof(double));
|
memcpy(cord_vector, ctx->jastrow.cord_vector, ctx->jastrow.dim_cord_vect*sizeof(double));
|
||||||
return QMCKL_SUCCESS;
|
return QMCKL_SUCCESS;
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -860,19 +860,22 @@ qmckl_exit_code qmckl_set_jastrow_cord_vector(qmckl_context context, double cons
|
|||||||
|
|
||||||
int32_t mask = 1 << 5;
|
int32_t mask = 1 << 5;
|
||||||
|
|
||||||
int64_t cord_num;
|
qmckl_exit_code rc = qmckl_provide_dim_cord_vect(context);
|
||||||
qmckl_exit_code rc = qmckl_get_jastrow_cord_num(context, &cord_num);
|
if (rc != QMCKL_SUCCESS) return rc;
|
||||||
|
|
||||||
|
int64_t dim_cord_vect;
|
||||||
|
rc = qmckl_get_jastrow_dim_cord_vect(context, &dim_cord_vect);
|
||||||
if (rc != QMCKL_SUCCESS) return rc;
|
if (rc != QMCKL_SUCCESS) return rc;
|
||||||
|
|
||||||
int64_t type_nucl_num;
|
int64_t type_nucl_num;
|
||||||
rc = qmckl_get_jastrow_type_nucl_num(context, &type_nucl_num);
|
rc = qmckl_get_jastrow_type_nucl_num(context, &type_nucl_num);
|
||||||
if (rc != QMCKL_SUCCESS) return rc;
|
if (rc != QMCKL_SUCCESS) return rc;
|
||||||
|
|
||||||
if (cord_num == 0) {
|
if (dim_cord_vect == 0) {
|
||||||
return qmckl_failwith( context,
|
return qmckl_failwith( context,
|
||||||
QMCKL_FAILURE,
|
QMCKL_FAILURE,
|
||||||
"qmckl_set_jastrow_coefficient",
|
"qmckl_set_jastrow_coefficient",
|
||||||
"cord_num is not set");
|
"dim_cord_vect is not set");
|
||||||
}
|
}
|
||||||
|
|
||||||
if (cord_vector == NULL) {
|
if (cord_vector == NULL) {
|
||||||
@ -892,7 +895,7 @@ qmckl_exit_code qmckl_set_jastrow_cord_vector(qmckl_context context, double cons
|
|||||||
}
|
}
|
||||||
|
|
||||||
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 = cord_num * type_nucl_num * sizeof(double);
|
mem_info.size = dim_cord_vect * type_nucl_num * sizeof(double);
|
||||||
double* new_array = (double*) qmckl_malloc(context, mem_info);
|
double* new_array = (double*) qmckl_malloc(context, mem_info);
|
||||||
|
|
||||||
if(new_array == NULL) {
|
if(new_array == NULL) {
|
||||||
@ -1324,20 +1327,20 @@ end function qmckl_compute_asymp_jasb_f
|
|||||||
#+CALL: generate_c_header(table=qmckl_asymp_jasb_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
#+CALL: generate_c_header(table=qmckl_asymp_jasb_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
||||||
|
|
||||||
#+RESULTS:
|
#+RESULTS:
|
||||||
#+begin_src c :tangle (eval h_private_func) :comments org
|
#+BEGIN_src c :tangle (eval h_func) :comments org
|
||||||
qmckl_exit_code qmckl_compute_asymp_jasb (
|
qmckl_exit_code qmckl_compute_asymp_jasb (
|
||||||
const qmckl_context context,
|
const qmckl_context context,
|
||||||
const int64_t bord_num,
|
const int64_t bord_num,
|
||||||
const double* bord_vector,
|
const double* bord_vector,
|
||||||
const double rescale_factor_kappa_ee,
|
const double rescale_factor_kappa_ee,
|
||||||
double* const asymp_jasb );
|
double* const asymp_jasb );
|
||||||
#+end_src
|
#+END_src
|
||||||
|
|
||||||
|
|
||||||
#+CALL: generate_c_interface(table=qmckl_asymp_jasb_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
#+CALL: generate_c_interface(table=qmckl_asymp_jasb_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
||||||
|
|
||||||
#+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_asymp_jasb &
|
integer(c_int32_t) function qmckl_compute_asymp_jasb &
|
||||||
(context, bord_num, bord_vector, rescale_factor_kappa_ee, asymp_jasb) &
|
(context, bord_num, bord_vector, rescale_factor_kappa_ee, asymp_jasb) &
|
||||||
bind(C) result(info)
|
bind(C) result(info)
|
||||||
@ -1356,7 +1359,7 @@ end function qmckl_compute_asymp_jasb_f
|
|||||||
(context, bord_num, bord_vector, rescale_factor_kappa_ee, asymp_jasb)
|
(context, bord_num, bord_vector, rescale_factor_kappa_ee, asymp_jasb)
|
||||||
|
|
||||||
end function qmckl_compute_asymp_jasb
|
end function qmckl_compute_asymp_jasb
|
||||||
#+end_src
|
#+END_src
|
||||||
|
|
||||||
*** Test
|
*** Test
|
||||||
#+name: asymp_jasb
|
#+name: asymp_jasb
|
||||||
@ -1380,11 +1383,6 @@ print("asymp_jasb[1] : ", asymp_jasb[1])
|
|||||||
#+end_src
|
#+end_src
|
||||||
|
|
||||||
#+RESULTS: asymp_jasb
|
#+RESULTS: asymp_jasb
|
||||||
: asym_one : 0.6634291325000664
|
|
||||||
: asymp_jasb[0] : 1.043287918508297
|
|
||||||
: asymp_jasb[1] : 0.7115733522582638
|
|
||||||
|
|
||||||
#+RESULTS:
|
|
||||||
: asym_one : 0.43340325572525706
|
: asym_one : 0.43340325572525706
|
||||||
: asymp_jasb[0] : 0.5323750557252571
|
: asymp_jasb[0] : 0.5323750557252571
|
||||||
: asymp_jasb[1] : 0.31567342786262853
|
: asymp_jasb[1] : 0.31567342786262853
|
||||||
@ -1416,6 +1414,8 @@ rc = qmckl_set_jastrow_aord_vector(context, aord_vector);
|
|||||||
assert(rc == QMCKL_SUCCESS);
|
assert(rc == QMCKL_SUCCESS);
|
||||||
rc = qmckl_set_jastrow_bord_vector(context, bord_vector);
|
rc = qmckl_set_jastrow_bord_vector(context, bord_vector);
|
||||||
assert(rc == QMCKL_SUCCESS);
|
assert(rc == QMCKL_SUCCESS);
|
||||||
|
rc = qmckl_set_jastrow_bord_vector(context, bord_vector);
|
||||||
|
assert(rc == QMCKL_SUCCESS);
|
||||||
rc = qmckl_set_jastrow_cord_vector(context, cord_vector);
|
rc = qmckl_set_jastrow_cord_vector(context, cord_vector);
|
||||||
assert(rc == QMCKL_SUCCESS);
|
assert(rc == QMCKL_SUCCESS);
|
||||||
rc = qmckl_set_jastrow_dependencies(context);
|
rc = qmckl_set_jastrow_dependencies(context);
|
||||||
@ -2114,28 +2114,16 @@ print("factor_ee_deriv_e[0][0]:",factor_ee_deriv_e[0][0])
|
|||||||
print("factor_ee_deriv_e[1][0]:",factor_ee_deriv_e[1][0])
|
print("factor_ee_deriv_e[1][0]:",factor_ee_deriv_e[1][0])
|
||||||
print("factor_ee_deriv_e[2][0]:",factor_ee_deriv_e[2][0])
|
print("factor_ee_deriv_e[2][0]:",factor_ee_deriv_e[2][0])
|
||||||
print("factor_ee_deriv_e[3][0]:",factor_ee_deriv_e[3][0])
|
print("factor_ee_deriv_e[3][0]:",factor_ee_deriv_e[3][0])
|
||||||
print(factor_ee_deriv_e)
|
|
||||||
|
|
||||||
#+end_src
|
#+end_src
|
||||||
|
|
||||||
#+RESULTS:
|
#+RESULTS:
|
||||||
#+begin_example
|
: asym_one : 0.43340325572525706
|
||||||
asym_one : 0.43340325572525706
|
: asymp_jasb[0] : 0.5323750557252571
|
||||||
asymp_jasb[0] : 0.5323750557252571
|
: asymp_jasb[1] : 0.31567342786262853
|
||||||
asymp_jasb[1] : 0.31567342786262853
|
: factor_ee_deriv_e[0][0]: 0.16364894652107934
|
||||||
factor_ee_deriv_e[0][0]: 0.16364894652107934
|
: factor_ee_deriv_e[1][0]: -0.6927548119830084
|
||||||
factor_ee_deriv_e[1][0]: -0.6927548119830084
|
: factor_ee_deriv_e[2][0]: 0.073267755223968
|
||||||
factor_ee_deriv_e[2][0]: 0.073267755223968
|
: factor_ee_deriv_e[3][0]: 1.5111672803213185
|
||||||
factor_ee_deriv_e[3][0]: 1.5111672803213185
|
|
||||||
[[ 0.16364895 0.60354957 -0.19825547 0.02359797 -0.13123153 -0.18789233
|
|
||||||
0.07762515 -0.42459184 0.27920265 -0.2056531 ]
|
|
||||||
[-0.69275481 0.15690393 0.09831069 0.18490587 0.04361723 0.3250686
|
|
||||||
0.12657961 -0.01736522 -0.40149005 0.17622416]
|
|
||||||
[ 0.07326776 -0.27532276 0.22396943 0.18771633 -0.34506246 0.07298062
|
|
||||||
0.63302352 -0.00910198 -0.30238713 -0.25908332]
|
|
||||||
[ 1.51116728 1.5033247 0.00325003 2.89377255 0.1338393 2.15893795
|
|
||||||
1.74732003 0.23561147 2.67455607 0.82810434]]
|
|
||||||
#+end_example
|
|
||||||
|
|
||||||
|
|
||||||
#+begin_src c :tangle (eval c_test)
|
#+begin_src c :tangle (eval c_test)
|
||||||
@ -3039,6 +3027,13 @@ integer function qmckl_compute_een_rescaled_e_f(context, walk_num, elec_num, cor
|
|||||||
end do
|
end do
|
||||||
end do
|
end do
|
||||||
end do
|
end do
|
||||||
|
|
||||||
|
do l = 0, cord_num
|
||||||
|
do j = 1, elec_num
|
||||||
|
een_rescaled_e(l, j, j, nw) = 0.0d0
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
|
||||||
end do
|
end do
|
||||||
|
|
||||||
end function qmckl_compute_een_rescaled_e_f
|
end function qmckl_compute_een_rescaled_e_f
|
||||||
@ -3124,6 +3119,10 @@ for l in range(1,cord_num+1):
|
|||||||
een_rescaled_e[j, i, l] = x
|
een_rescaled_e[j, i, l] = x
|
||||||
k = k + 1
|
k = k + 1
|
||||||
|
|
||||||
|
for l in range(0,cord_num+1):
|
||||||
|
for j in range(0, elec_num):
|
||||||
|
een_rescaled_e[j,j,l] = 0.0
|
||||||
|
|
||||||
print(" een_rescaled_e[0, 2, 1] = ",een_rescaled_e[0, 2, 1])
|
print(" een_rescaled_e[0, 2, 1] = ",een_rescaled_e[0, 2, 1])
|
||||||
print(" een_rescaled_e[0, 3, 1] = ",een_rescaled_e[0, 3, 1])
|
print(" een_rescaled_e[0, 3, 1] = ",een_rescaled_e[0, 3, 1])
|
||||||
print(" een_rescaled_e[0, 4, 1] = ",een_rescaled_e[0, 4, 1])
|
print(" een_rescaled_e[0, 4, 1] = ",een_rescaled_e[0, 4, 1])
|
||||||
@ -3135,7 +3134,7 @@ print(" een_rescaled_e[1, 5, 2] = ",een_rescaled_e[1, 5, 2])
|
|||||||
#+RESULTS:
|
#+RESULTS:
|
||||||
: een_rescaled_e[0, 2, 1] = 0.08084493981483197
|
: een_rescaled_e[0, 2, 1] = 0.08084493981483197
|
||||||
: een_rescaled_e[0, 3, 1] = 0.1066745707571846
|
: een_rescaled_e[0, 3, 1] = 0.1066745707571846
|
||||||
: een_rescaled_e[0, 4, 1] = 0.01754273169464735
|
: een_rescaled_e[0, 4, 1] = 0.017542731694647366
|
||||||
: een_rescaled_e[1, 3, 2] = 0.02214680362033448
|
: een_rescaled_e[1, 3, 2] = 0.02214680362033448
|
||||||
: een_rescaled_e[1, 4, 2] = 0.0005700154999202759
|
: een_rescaled_e[1, 4, 2] = 0.0005700154999202759
|
||||||
: een_rescaled_e[1, 5, 2] = 0.3424402276009091
|
: een_rescaled_e[1, 5, 2] = 0.3424402276009091
|
||||||
@ -3159,9 +3158,10 @@ assert(fabs(een_rescaled_e[0][1][5][2]-0.3424402276009091) < 1.e-12);
|
|||||||
|
|
||||||
** Electron-electron rescaled distances for each order and derivatives
|
** Electron-electron rescaled distances for each order and derivatives
|
||||||
|
|
||||||
~een_rescaled_e~ stores the table of the rescaled distances between all
|
~een_rescaled_e_deriv_e~ stores the table of the derivatives of the
|
||||||
pairs of electrons and raised to the power \(p\) defined by ~cord_num~.
|
rescaled distances between all pairs of electrons and raised to the
|
||||||
Here we take its derivatives required for the een jastrow.
|
power \(p\) defined by ~cord_num~. Here we take its derivatives
|
||||||
|
required for the een jastrow.
|
||||||
|
|
||||||
TODO: write formulae
|
TODO: write formulae
|
||||||
|
|
||||||
@ -3419,7 +3419,7 @@ end function qmckl_compute_factor_een_rescaled_e_deriv_e_f
|
|||||||
#+end_src
|
#+end_src
|
||||||
|
|
||||||
*** Test
|
*** Test
|
||||||
|
#+name: een_e_deriv_e
|
||||||
#+begin_src python :results output :exports none :noweb yes
|
#+begin_src python :results output :exports none :noweb yes
|
||||||
import numpy as np
|
import numpy as np
|
||||||
|
|
||||||
@ -3431,6 +3431,16 @@ for i in range(elec_num):
|
|||||||
for j in range(elec_num):
|
for j in range(elec_num):
|
||||||
elec_dist[i, j] = np.linalg.norm(elec_coord[i] - elec_coord[j])
|
elec_dist[i, j] = np.linalg.norm(elec_coord[i] - elec_coord[j])
|
||||||
|
|
||||||
|
elec_dist_deriv_e = np.zeros(shape=(4,elec_num, elec_num),dtype=float)
|
||||||
|
for j in range(elec_num):
|
||||||
|
for i in range(elec_num):
|
||||||
|
rij_inv = 1.0 / elec_dist[i, j]
|
||||||
|
for ii in range(3):
|
||||||
|
elec_dist_deriv_e[ii, i, j] = -(elec_coord[j][ii] - elec_coord[i][ii]) * rij_inv
|
||||||
|
elec_dist_deriv_e[3, i, j] = 2.0 * rij_inv
|
||||||
|
elec_dist_deriv_e[:, j, j] = 0.0
|
||||||
|
|
||||||
|
|
||||||
kappa = 1.0
|
kappa = 1.0
|
||||||
|
|
||||||
een_rescaled_e_ij = np.zeros(shape=(elec_num * (elec_num - 1)//2, cord_num+1), dtype=float)
|
een_rescaled_e_ij = np.zeros(shape=(elec_num * (elec_num - 1)//2, cord_num+1), dtype=float)
|
||||||
@ -3458,25 +3468,49 @@ for l in range(1,cord_num+1):
|
|||||||
een_rescaled_e[j, i, l] = x
|
een_rescaled_e[j, i, l] = x
|
||||||
k = k + 1
|
k = k + 1
|
||||||
|
|
||||||
print(" een_rescaled_e[0, 2, 1] = ",een_rescaled_e[0, 2, 1])
|
een_rescaled_e_deriv_e = np.zeros(shape=(elec_num,4,elec_num,cord_num+1),dtype=float)
|
||||||
print(" een_rescaled_e[0, 3, 1] = ",een_rescaled_e[0, 3, 1])
|
for l in range(0,cord_num+1):
|
||||||
print(" een_rescaled_e[0, 4, 1] = ",een_rescaled_e[0, 4, 1])
|
kappa_l = -1.0 * kappa * l
|
||||||
print(" een_rescaled_e[1, 3, 2] = ",een_rescaled_e[1, 3, 2])
|
for j in range(0,elec_num):
|
||||||
print(" een_rescaled_e[1, 4, 2] = ",een_rescaled_e[1, 4, 2])
|
for i in range(0,elec_num):
|
||||||
print(" een_rescaled_e[1, 5, 2] = ",een_rescaled_e[1, 5, 2])
|
for ii in range(0,4):
|
||||||
|
een_rescaled_e_deriv_e[i,ii,j,l] = kappa_l * elec_dist_deriv_e[ii,i,j]
|
||||||
|
een_rescaled_e_deriv_e[i,3,j,l] = een_rescaled_e_deriv_e[i,3,j,l] + \
|
||||||
|
een_rescaled_e_deriv_e[i,0,j,l] * een_rescaled_e_deriv_e[i,0,j,l] + \
|
||||||
|
een_rescaled_e_deriv_e[i,1,j,l] * een_rescaled_e_deriv_e[i,1,j,l] + \
|
||||||
|
een_rescaled_e_deriv_e[i,2,j,l] * een_rescaled_e_deriv_e[i,2,j,l]
|
||||||
|
|
||||||
|
for ii in range(0,4):
|
||||||
|
een_rescaled_e_deriv_e[i,ii,j,l] = een_rescaled_e_deriv_e[i,ii,j,l] * een_rescaled_e[i,j,l]
|
||||||
|
|
||||||
|
#print(" een_rescaled_e_deriv_e[1, 1, 3, 1] = ",een_rescaled_e_deriv_e[0, 0, 2, 1])
|
||||||
|
#print(" een_rescaled_e_deriv_e[1, 1, 4, 1] = ",een_rescaled_e_deriv_e[0, 0, 3, 1])
|
||||||
|
#print(" een_rescaled_e_deriv_e[1, 1, 5, 1] = ",een_rescaled_e_deriv_e[0, 0, 4, 1])
|
||||||
|
#print(" een_rescaled_e_deriv_e[2, 1, 4, 2] = ",een_rescaled_e_deriv_e[1, 0, 3, 2])
|
||||||
|
#print(" een_rescaled_e_deriv_e[2, 1, 5, 2] = ",een_rescaled_e_deriv_e[1, 0, 4, 2])
|
||||||
|
#print(" een_rescaled_e_deriv_e[2, 1, 6, 2] = ",een_rescaled_e_deriv_e[1, 0, 5, 2])
|
||||||
#+end_src
|
#+end_src
|
||||||
|
|
||||||
#+RESULTS:
|
#+RESULTS: een_e_deriv_e
|
||||||
: een_rescaled_e[0, 2, 1] = 0.08084493981483197
|
: een_rescaled_e_deriv_e[1, 1, 3, 1] = 0.05991352796887283
|
||||||
: een_rescaled_e[0, 3, 1] = 0.1066745707571846
|
: een_rescaled_e_deriv_e[1, 1, 4, 1] = 0.011714035071545248
|
||||||
: een_rescaled_e[0, 4, 1] = 0.01754273169464735
|
: een_rescaled_e_deriv_e[1, 1, 5, 1] = 0.00441398875758468
|
||||||
: een_rescaled_e[1, 3, 2] = 0.02214680362033448
|
: een_rescaled_e_deriv_e[2, 1, 4, 2] = 0.013553180060167595
|
||||||
: een_rescaled_e[1, 4, 2] = 0.0005700154999202759
|
: een_rescaled_e_deriv_e[2, 1, 5, 2] = 0.00041342909359870457
|
||||||
: een_rescaled_e[1, 5, 2] = 0.3424402276009091
|
: een_rescaled_e_deriv_e[2, 1, 6, 2] = 0.5880599146214673
|
||||||
|
|
||||||
#+begin_src c :tangle (eval c_test)
|
#+begin_src c :tangle (eval c_test)
|
||||||
//assert(qmckl_electron_provided(context));
|
//assert(qmckl_electron_provided(context));
|
||||||
|
double een_rescaled_e_deriv_e[walk_num][elec_num][4][elec_num][(cord_num + 1)];
|
||||||
|
rc = qmckl_get_jastrow_een_rescaled_e_deriv_e(context, &(een_rescaled_e_deriv_e[0][0][0][0][0]));
|
||||||
|
|
||||||
|
// value of (0,0,0,2,1)
|
||||||
|
assert(fabs(een_rescaled_e_deriv_e[0][0][0][2][1] + 0.05991352796887283 ) < 1.e-12);
|
||||||
|
assert(fabs(een_rescaled_e_deriv_e[0][0][0][3][1] + 0.011714035071545248 ) < 1.e-12);
|
||||||
|
assert(fabs(een_rescaled_e_deriv_e[0][0][0][4][1] + 0.00441398875758468 ) < 1.e-12);
|
||||||
|
assert(fabs(een_rescaled_e_deriv_e[0][1][0][3][2] + 0.013553180060167595 ) < 1.e-12);
|
||||||
|
assert(fabs(een_rescaled_e_deriv_e[0][1][0][4][2] + 0.00041342909359870457) < 1.e-12);
|
||||||
|
assert(fabs(een_rescaled_e_deriv_e[0][1][0][5][2] + 0.5880599146214673 ) < 1.e-12);
|
||||||
#+end_src
|
#+end_src
|
||||||
|
|
||||||
** Electron-nucleus rescaled distances for each order
|
** Electron-nucleus rescaled distances for each order
|
||||||
@ -4076,6 +4110,14 @@ for i in range(elec_num):
|
|||||||
for a in range(nucl_num):
|
for a in range(nucl_num):
|
||||||
elnuc_dist[i, a] = np.linalg.norm(elec_coord[i] - nucl_coord[:,a])
|
elnuc_dist[i, a] = np.linalg.norm(elec_coord[i] - nucl_coord[:,a])
|
||||||
|
|
||||||
|
elnuc_dist_deriv_e = np.zeros(shape=(4, elec_num, nucl_num),dtype=float)
|
||||||
|
for a in range(nucl_num):
|
||||||
|
for i in range(elec_num):
|
||||||
|
rij_inv = 1.0 / elnuc_dist[i, a]
|
||||||
|
for ii in range(3):
|
||||||
|
elnuc_dist_deriv_e[ii, i, a] = (elec_coord[i][ii] - nucl_coord[ii][a]) * rij_inv
|
||||||
|
elnuc_dist_deriv_e[3, i, a] = 2.0 * rij_inv
|
||||||
|
|
||||||
kappa = 1.0
|
kappa = 1.0
|
||||||
|
|
||||||
een_rescaled_n = np.zeros(shape=(nucl_num, elec_num, cord_num + 1), dtype=float)
|
een_rescaled_n = np.zeros(shape=(nucl_num, elec_num, cord_num + 1), dtype=float)
|
||||||
@ -4090,24 +4132,50 @@ for l in range(2,cord_num+1):
|
|||||||
for i in range(elec_num):
|
for i in range(elec_num):
|
||||||
een_rescaled_n[a, i, l] = een_rescaled_n[a, i, l - 1] * een_rescaled_n[a, i, 1]
|
een_rescaled_n[a, i, l] = een_rescaled_n[a, i, l - 1] * een_rescaled_n[a, i, 1]
|
||||||
|
|
||||||
print(" een_rescaled_n[0, 2, 1] = ",een_rescaled_n[0, 2, 1])
|
een_rescaled_n_deriv_e = np.zeros(shape=(elec_num,4,nucl_num,cord_num+1),dtype=float)
|
||||||
print(" een_rescaled_n[0, 3, 1] = ",een_rescaled_n[0, 3, 1])
|
for l in range(0,cord_num+1):
|
||||||
print(" een_rescaled_n[0, 4, 1] = ",een_rescaled_n[0, 4, 1])
|
kappa_l = -1.0 * kappa * l
|
||||||
print(" een_rescaled_n[1, 3, 2] = ",een_rescaled_n[1, 3, 2])
|
for j in range(0,elec_num):
|
||||||
print(" een_rescaled_n[1, 4, 2] = ",een_rescaled_n[1, 4, 2])
|
for a in range(0,nucl_num):
|
||||||
print(" een_rescaled_n[1, 5, 2] = ",een_rescaled_n[1, 5, 2])
|
for ii in range(0,4):
|
||||||
|
een_rescaled_n_deriv_e[j,ii,a,l] = kappa_l * elnuc_dist_deriv_e[ii,j,a]
|
||||||
|
een_rescaled_n_deriv_e[j,3,a,l] = een_rescaled_n_deriv_e[j,3,a,l] + \
|
||||||
|
een_rescaled_n_deriv_e[j,0,a,l] * een_rescaled_n_deriv_e[j,0,a,l] + \
|
||||||
|
een_rescaled_n_deriv_e[j,1,a,l] * een_rescaled_n_deriv_e[j,1,a,l] + \
|
||||||
|
een_rescaled_n_deriv_e[j,2,a,l] * een_rescaled_n_deriv_e[j,2,a,l]
|
||||||
|
|
||||||
|
for ii in range(0,4):
|
||||||
|
een_rescaled_n_deriv_e[j,ii,a,l] = een_rescaled_n_deriv_e[j,ii,a,l] * een_rescaled_n[a,j,l]
|
||||||
|
|
||||||
|
print(" een_rescaled_n_deriv_e[1, 1, 3, 1] = ",een_rescaled_n_deriv_e[2, 0, 0, 1])
|
||||||
|
print(" een_rescaled_n_deriv_e[1, 1, 4, 1] = ",een_rescaled_n_deriv_e[3, 0, 0, 1])
|
||||||
|
print(" een_rescaled_n_deriv_e[1, 1, 5, 1] = ",een_rescaled_n_deriv_e[4, 0, 0, 1])
|
||||||
|
print(" een_rescaled_n_deriv_e[2, 1, 4, 2] = ",een_rescaled_n_deriv_e[3, 0, 1, 2])
|
||||||
|
print(" een_rescaled_n_deriv_e[2, 1, 5, 2] = ",een_rescaled_n_deriv_e[4, 0, 1, 2])
|
||||||
|
print(" een_rescaled_n_deriv_e[2, 1, 6, 2] = ",een_rescaled_n_deriv_e[5, 0, 1, 2])
|
||||||
#+end_src
|
#+end_src
|
||||||
|
|
||||||
#+RESULTS:
|
#+RESULTS:
|
||||||
: een_rescaled_n[0, 2, 1] = 0.10612983920006765
|
: een_rescaled_n_deriv_e[1, 1, 3, 1] = -0.07633444246999128
|
||||||
: een_rescaled_n[0, 3, 1] = 0.135652809635553
|
: een_rescaled_n_deriv_e[1, 1, 4, 1] = 0.00033282346259738276
|
||||||
: een_rescaled_n[0, 4, 1] = 0.023391817607642338
|
: een_rescaled_n_deriv_e[1, 1, 5, 1] = -0.004775370547333061
|
||||||
: een_rescaled_n[1, 3, 2] = 0.880957224822116
|
: een_rescaled_n_deriv_e[2, 1, 4, 2] = 0.1362654644223866
|
||||||
: een_rescaled_n[1, 4, 2] = 0.027185942659395074
|
: een_rescaled_n_deriv_e[2, 1, 5, 2] = -0.0231253431662794
|
||||||
: een_rescaled_n[1, 5, 2] = 0.01343938025140174
|
: een_rescaled_n_deriv_e[2, 1, 6, 2] = 0.001593334817691633
|
||||||
|
|
||||||
#+begin_src c :tangle (eval c_test)
|
#+begin_src c :tangle (eval c_test)
|
||||||
//assert(qmckl_electron_provided(context));
|
assert(qmckl_electron_provided(context));
|
||||||
|
|
||||||
|
double een_rescaled_n_deriv_e[walk_num][elec_num][4][nucl_num][(cord_num + 1)];
|
||||||
|
rc = qmckl_get_jastrow_een_rescaled_n_deriv_e(context, &(een_rescaled_n_deriv_e[0][0][0][0][0]));
|
||||||
|
|
||||||
|
// value of (0,2,1)
|
||||||
|
assert(fabs(een_rescaled_n_deriv_e[0][2][0][0][1]+0.07633444246999128 ) < 1.e-12);
|
||||||
|
assert(fabs(een_rescaled_n_deriv_e[0][3][0][0][1]-0.00033282346259738276) < 1.e-12);
|
||||||
|
assert(fabs(een_rescaled_n_deriv_e[0][4][0][0][1]+0.004775370547333061 ) < 1.e-12);
|
||||||
|
assert(fabs(een_rescaled_n_deriv_e[0][3][0][1][2]-0.1362654644223866 ) < 1.e-12);
|
||||||
|
assert(fabs(een_rescaled_n_deriv_e[0][4][0][1][2]+0.0231253431662794 ) < 1.e-12);
|
||||||
|
assert(fabs(een_rescaled_n_deriv_e[0][5][0][1][2]-0.001593334817691633 ) < 1.e-12);
|
||||||
|
|
||||||
#+end_src
|
#+end_src
|
||||||
|
|
||||||
@ -4264,7 +4332,6 @@ qmckl_exit_code qmckl_provide_cord_vect_full(qmckl_context context)
|
|||||||
qmckl_exit_code rc =
|
qmckl_exit_code rc =
|
||||||
qmckl_compute_cord_vect_full(context,
|
qmckl_compute_cord_vect_full(context,
|
||||||
ctx->nucleus.num,
|
ctx->nucleus.num,
|
||||||
ctx->jastrow.cord_num,
|
|
||||||
ctx->jastrow.dim_cord_vect,
|
ctx->jastrow.dim_cord_vect,
|
||||||
ctx->jastrow.type_nucl_num,
|
ctx->jastrow.type_nucl_num,
|
||||||
ctx->jastrow.type_nucl_vector,
|
ctx->jastrow.type_nucl_vector,
|
||||||
@ -4426,26 +4493,24 @@ end function qmckl_compute_dim_cord_vect_f
|
|||||||
#+NAME: qmckl_factor_cord_vect_full_args
|
#+NAME: qmckl_factor_cord_vect_full_args
|
||||||
| qmckl_context | context | in | Global state |
|
| qmckl_context | context | in | Global state |
|
||||||
| int64_t | nucl_num | in | Number of atoms |
|
| int64_t | nucl_num | in | Number of atoms |
|
||||||
| int64_t | cord_num | in | Order of polynomials |
|
|
||||||
| int64_t | dim_cord_vect | in | dimension of cord full table |
|
| int64_t | dim_cord_vect | in | dimension of cord full table |
|
||||||
| int64_t | type_nucl_num | in | dimension of cord full table |
|
| int64_t | type_nucl_num | in | dimension of cord full table |
|
||||||
| int64_t | type_nucl_vector[nucl_num] | in | dimension of cord full table |
|
| int64_t | type_nucl_vector[nucl_num] | in | dimension of cord full table |
|
||||||
| double | cord_vector[cord_num][type_nucl_num] | in | dimension of cord full table |
|
| double | cord_vector[dim_cord_vect][type_nucl_num] | in | dimension of cord full table |
|
||||||
| double | cord_vect_full[dim_cord_vect][nucl_num] | out | Full list of coefficients |
|
| double | cord_vect_full[dim_cord_vect][nucl_num] | out | Full list of coefficients |
|
||||||
|
|
||||||
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
|
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
|
||||||
integer function qmckl_compute_cord_vect_full_f(context, nucl_num, cord_num, dim_cord_vect, type_nucl_num, &
|
integer function qmckl_compute_cord_vect_full_f(context, nucl_num, dim_cord_vect, type_nucl_num, &
|
||||||
type_nucl_vector, cord_vector, cord_vect_full) &
|
type_nucl_vector, cord_vector, cord_vect_full) &
|
||||||
result(info)
|
result(info)
|
||||||
use qmckl
|
use qmckl
|
||||||
implicit none
|
implicit none
|
||||||
integer(qmckl_context), intent(in) :: context
|
integer(qmckl_context), intent(in) :: context
|
||||||
integer*8 , intent(in) :: nucl_num
|
integer*8 , intent(in) :: nucl_num
|
||||||
integer*8 , intent(in) :: cord_num
|
|
||||||
integer*8 , intent(in) :: dim_cord_vect
|
integer*8 , intent(in) :: dim_cord_vect
|
||||||
integer*8 , intent(in) :: type_nucl_num
|
integer*8 , intent(in) :: type_nucl_num
|
||||||
integer*8 , intent(in) :: type_nucl_vector(nucl_num)
|
integer*8 , intent(in) :: type_nucl_vector(nucl_num)
|
||||||
double precision , intent(in) :: cord_vector(cord_num, type_nucl_num)
|
double precision , intent(in) :: cord_vector(type_nucl_num, dim_cord_vect)
|
||||||
double precision , intent(out) :: cord_vect_full(nucl_num,dim_cord_vect)
|
double precision , intent(out) :: cord_vect_full(nucl_num,dim_cord_vect)
|
||||||
double precision :: x
|
double precision :: x
|
||||||
integer*8 :: i, a, k, l, nw
|
integer*8 :: i, a, k, l, nw
|
||||||
@ -4462,11 +4527,6 @@ integer function qmckl_compute_cord_vect_full_f(context, nucl_num, cord_num, dim
|
|||||||
return
|
return
|
||||||
endif
|
endif
|
||||||
|
|
||||||
if (cord_num <= 0) then
|
|
||||||
info = QMCKL_INVALID_ARG_3
|
|
||||||
return
|
|
||||||
endif
|
|
||||||
|
|
||||||
if (type_nucl_num <= 0) then
|
if (type_nucl_num <= 0) then
|
||||||
info = QMCKL_INVALID_ARG_4
|
info = QMCKL_INVALID_ARG_4
|
||||||
return
|
return
|
||||||
@ -4479,7 +4539,7 @@ integer function qmckl_compute_cord_vect_full_f(context, nucl_num, cord_num, dim
|
|||||||
|
|
||||||
|
|
||||||
do a = 1, nucl_num
|
do a = 1, nucl_num
|
||||||
cord_vect_full(1:dim_cord_vect,a) = cord_vector(1:dim_cord_vect,type_nucl_vector(a))
|
cord_vect_full(a,1:dim_cord_vect) = cord_vector(type_nucl_vector(a),1:dim_cord_vect)
|
||||||
end do
|
end do
|
||||||
|
|
||||||
end function qmckl_compute_cord_vect_full_f
|
end function qmckl_compute_cord_vect_full_f
|
||||||
@ -4492,7 +4552,6 @@ end function qmckl_compute_cord_vect_full_f
|
|||||||
qmckl_exit_code qmckl_compute_cord_vect_full (
|
qmckl_exit_code qmckl_compute_cord_vect_full (
|
||||||
const qmckl_context context,
|
const qmckl_context context,
|
||||||
const int64_t nucl_num,
|
const int64_t nucl_num,
|
||||||
const int64_t cord_num,
|
|
||||||
const int64_t dim_cord_vect,
|
const int64_t dim_cord_vect,
|
||||||
const int64_t type_nucl_num,
|
const int64_t type_nucl_num,
|
||||||
const int64_t* type_nucl_vector,
|
const int64_t* type_nucl_vector,
|
||||||
@ -4506,14 +4565,7 @@ end function qmckl_compute_cord_vect_full_f
|
|||||||
#+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_cord_vect_full &
|
integer(c_int32_t) function qmckl_compute_cord_vect_full &
|
||||||
(context, &
|
(context, nucl_num, dim_cord_vect, type_nucl_num, type_nucl_vector, cord_vector, cord_vect_full) &
|
||||||
nucl_num, &
|
|
||||||
cord_num, &
|
|
||||||
dim_cord_vect, &
|
|
||||||
type_nucl_num, &
|
|
||||||
type_nucl_vector, &
|
|
||||||
cord_vector, &
|
|
||||||
cord_vect_full) &
|
|
||||||
bind(C) result(info)
|
bind(C) result(info)
|
||||||
|
|
||||||
use, intrinsic :: iso_c_binding
|
use, intrinsic :: iso_c_binding
|
||||||
@ -4521,23 +4573,15 @@ end function qmckl_compute_cord_vect_full_f
|
|||||||
|
|
||||||
integer (c_int64_t) , intent(in) , value :: context
|
integer (c_int64_t) , intent(in) , value :: context
|
||||||
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 :: cord_num
|
|
||||||
integer (c_int64_t) , intent(in) , value :: dim_cord_vect
|
integer (c_int64_t) , intent(in) , value :: dim_cord_vect
|
||||||
integer (c_int64_t) , intent(in) , value :: type_nucl_num
|
integer (c_int64_t) , intent(in) , value :: type_nucl_num
|
||||||
integer (c_int64_t) , intent(in) :: type_nucl_vector(nucl_num)
|
integer (c_int64_t) , intent(in) :: type_nucl_vector(nucl_num)
|
||||||
real (c_double ) , intent(in) :: cord_vector(type_nucl_num,cord_num)
|
real (c_double ) , intent(in) :: cord_vector(type_nucl_num,dim_cord_vect)
|
||||||
real (c_double ) , intent(out) :: cord_vect_full(nucl_num,dim_cord_vect)
|
real (c_double ) , intent(out) :: cord_vect_full(nucl_num,dim_cord_vect)
|
||||||
|
|
||||||
integer(c_int32_t), external :: qmckl_compute_cord_vect_full_f
|
integer(c_int32_t), external :: qmckl_compute_cord_vect_full_f
|
||||||
info = qmckl_compute_cord_vect_full_f &
|
info = qmckl_compute_cord_vect_full_f &
|
||||||
(context, &
|
(context, nucl_num, dim_cord_vect, type_nucl_num, type_nucl_vector, cord_vector, cord_vect_full)
|
||||||
nucl_num, &
|
|
||||||
cord_num, &
|
|
||||||
dim_cord_vect, &
|
|
||||||
type_nucl_num, &
|
|
||||||
type_nucl_vector, &
|
|
||||||
cord_vector, &
|
|
||||||
cord_vect_full)
|
|
||||||
|
|
||||||
end function qmckl_compute_cord_vect_full
|
end function qmckl_compute_cord_vect_full
|
||||||
#+end_src
|
#+end_src
|
||||||
@ -4644,9 +4688,9 @@ end function qmckl_compute_lkpm_combined_index_f
|
|||||||
end function qmckl_compute_lkpm_combined_index
|
end function qmckl_compute_lkpm_combined_index
|
||||||
#+end_src
|
#+end_src
|
||||||
|
|
||||||
|
|
||||||
*** Test
|
*** Test
|
||||||
|
|
||||||
|
#+name: helper_funcs
|
||||||
#+begin_src python :results output :exports none :noweb yes
|
#+begin_src python :results output :exports none :noweb yes
|
||||||
import numpy as np
|
import numpy as np
|
||||||
|
|
||||||
@ -4673,21 +4717,46 @@ for l in range(2,cord_num+1):
|
|||||||
for i in range(elec_num):
|
for i in range(elec_num):
|
||||||
een_rescaled_n[a, i, l] = een_rescaled_n[a, i, l - 1] * een_rescaled_n[a, i, 1]
|
een_rescaled_n[a, i, l] = een_rescaled_n[a, i, l - 1] * een_rescaled_n[a, i, 1]
|
||||||
|
|
||||||
print(" een_rescaled_n[0, 2, 1] = ",een_rescaled_n[0, 2, 1])
|
elec_dist = np.zeros(shape=(elec_num, elec_num),dtype=float)
|
||||||
print(" een_rescaled_n[0, 3, 1] = ",een_rescaled_n[0, 3, 1])
|
for i in range(elec_num):
|
||||||
print(" een_rescaled_n[0, 4, 1] = ",een_rescaled_n[0, 4, 1])
|
for j in range(elec_num):
|
||||||
print(" een_rescaled_n[1, 3, 2] = ",een_rescaled_n[1, 3, 2])
|
elec_dist[i, j] = np.linalg.norm(elec_coord[i] - elec_coord[j])
|
||||||
print(" een_rescaled_n[1, 4, 2] = ",een_rescaled_n[1, 4, 2])
|
|
||||||
print(" een_rescaled_n[1, 5, 2] = ",een_rescaled_n[1, 5, 2])
|
kappa = 1.0
|
||||||
|
|
||||||
|
een_rescaled_e_ij = np.zeros(shape=(elec_num * (elec_num - 1)//2, cord_num+1), dtype=float)
|
||||||
|
een_rescaled_e_ij[:,0] = 1.0
|
||||||
|
|
||||||
|
k = 0
|
||||||
|
for j in range(elec_num):
|
||||||
|
for i in range(j):
|
||||||
|
een_rescaled_e_ij[k, 1] = np.exp(-kappa * elec_dist[i, j])
|
||||||
|
k = k + 1
|
||||||
|
|
||||||
|
for l in range(2, cord_num + 1):
|
||||||
|
for k in range(elec_num * (elec_num - 1)//2):
|
||||||
|
een_rescaled_e_ij[k, l] = een_rescaled_e_ij[k, l - 1] * een_rescaled_e_ij[k, 1]
|
||||||
|
|
||||||
|
een_rescaled_e = np.zeros(shape=(elec_num, elec_num, cord_num + 1), dtype=float)
|
||||||
|
een_rescaled_e[:,:,0] = 1.0
|
||||||
|
|
||||||
|
for l in range(1,cord_num+1):
|
||||||
|
k = 0
|
||||||
|
for j in range(elec_num):
|
||||||
|
for i in range(j):
|
||||||
|
x = een_rescaled_e_ij[k, l]
|
||||||
|
een_rescaled_e[i, j, l] = x
|
||||||
|
een_rescaled_e[j, i, l] = x
|
||||||
|
k = k + 1
|
||||||
|
|
||||||
|
for l in range(0,cord_num+1):
|
||||||
|
for j in range(0, elec_num):
|
||||||
|
een_rescaled_e[j,j,l] = 0.0
|
||||||
|
|
||||||
|
lkpm_of_cindex = np.array(lkpm_combined_index).T
|
||||||
#+end_src
|
#+end_src
|
||||||
|
|
||||||
#+RESULTS:
|
#+RESULTS: helper_funcs
|
||||||
: een_rescaled_n[0, 2, 1] = 0.10612983920006765
|
|
||||||
: een_rescaled_n[0, 3, 1] = 0.135652809635553
|
|
||||||
: een_rescaled_n[0, 4, 1] = 0.023391817607642338
|
|
||||||
: een_rescaled_n[1, 3, 2] = 0.880957224822116
|
|
||||||
: een_rescaled_n[1, 4, 2] = 0.027185942659395074
|
|
||||||
: een_rescaled_n[1, 5, 2] = 0.01343938025140174
|
|
||||||
|
|
||||||
#+begin_src c :tangle (eval c_test)
|
#+begin_src c :tangle (eval c_test)
|
||||||
//assert(qmckl_electron_provided(context));
|
//assert(qmckl_electron_provided(context));
|
||||||
@ -4834,10 +4903,10 @@ integer function qmckl_compute_factor_een_f(context, walk_num, elec_num, nucl_nu
|
|||||||
implicit none
|
implicit none
|
||||||
integer(qmckl_context), intent(in) :: context
|
integer(qmckl_context), intent(in) :: context
|
||||||
integer*8 , intent(in) :: walk_num, elec_num, cord_num, nucl_num, dim_cord_vect
|
integer*8 , intent(in) :: walk_num, elec_num, cord_num, nucl_num, dim_cord_vect
|
||||||
integer*8 , intent(in) :: lkpm_combined_index(4,dim_cord_vect)
|
integer*8 , intent(in) :: lkpm_combined_index(dim_cord_vect,4)
|
||||||
double precision , intent(in) :: cord_vect_full(dim_cord_vect, nucl_num)
|
double precision , intent(in) :: cord_vect_full(nucl_num, dim_cord_vect)
|
||||||
double precision , intent(in) :: een_rescaled_e(walk_num, elec_num, elec_num, 0:cord_num)
|
double precision , intent(in) :: een_rescaled_e(0:cord_num, elec_num, elec_num, walk_num)
|
||||||
double precision , intent(in) :: een_rescaled_n(walk_num, elec_num, nucl_num, 0:cord_num)
|
double precision , intent(in) :: een_rescaled_n(0:cord_num, nucl_num, elec_num, walk_num)
|
||||||
double precision , intent(out) :: factor_een(walk_num)
|
double precision , intent(out) :: factor_een(walk_num)
|
||||||
|
|
||||||
integer*8 :: i, a, j, l, k, p, m, n, nw
|
integer*8 :: i, a, j, l, k, p, m, n, nw
|
||||||
@ -4874,23 +4943,27 @@ integer function qmckl_compute_factor_een_f(context, walk_num, elec_num, nucl_nu
|
|||||||
|
|
||||||
do nw =1, walk_num
|
do nw =1, walk_num
|
||||||
do n = 1, dim_cord_vect
|
do n = 1, dim_cord_vect
|
||||||
l = lkpm_combined_index(1, n)
|
l = lkpm_combined_index(n, 1)
|
||||||
k = lkpm_combined_index(2, n)
|
k = lkpm_combined_index(n, 2)
|
||||||
p = lkpm_combined_index(3, n)
|
p = lkpm_combined_index(n, 3)
|
||||||
m = lkpm_combined_index(4, n)
|
m = lkpm_combined_index(n, 4)
|
||||||
|
|
||||||
do a = 1, nucl_num
|
do a = 1, nucl_num
|
||||||
accu2 = 0.0d0
|
accu2 = 0.0d0
|
||||||
cn = cord_vect_full(n, a)
|
cn = cord_vect_full(a, n)
|
||||||
do j = 1, elec_num
|
do j = 1, elec_num
|
||||||
accu = 0.0d0
|
accu = 0.0d0
|
||||||
do i = 1, elec_num
|
do i = 1, elec_num
|
||||||
accu = accu + een_rescaled_e(nw, i, j, k) * &
|
accu = accu + een_rescaled_e(k,i,j,nw) * &
|
||||||
een_rescaled_n(nw, i, a, m)
|
een_rescaled_n(m,a,i,nw)
|
||||||
|
!if(nw .eq. 1) then
|
||||||
|
! print *,l,k,p,m,j,i,een_rescaled_e(k,i,j,nw), een_rescaled_n(m,a,i,nw), accu
|
||||||
|
!endif
|
||||||
end do
|
end do
|
||||||
accu2 = accu2 + accu * een_rescaled_n(nw, j, a, m + l)
|
accu2 = accu2 + accu * een_rescaled_n(m + l,a,j,nw)
|
||||||
|
!print *, l,m,nw,accu, accu2, een_rescaled_n(m + l, a, j, nw), cn, factor_een(nw)
|
||||||
end do
|
end do
|
||||||
factor_een(nw) = factor_een(nw) + accu2 + cn
|
factor_een(nw) = factor_een(nw) + accu2 * cn
|
||||||
end do
|
end do
|
||||||
end do
|
end do
|
||||||
end do
|
end do
|
||||||
@ -4916,7 +4989,6 @@ end function qmckl_compute_factor_een_f
|
|||||||
double* const factor_een );
|
double* const factor_een );
|
||||||
#+end_src
|
#+end_src
|
||||||
|
|
||||||
|
|
||||||
#+CALL: generate_c_interface(table=qmckl_factor_een_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
#+CALL: generate_c_interface(table=qmckl_factor_een_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
||||||
|
|
||||||
#+RESULTS:
|
#+RESULTS:
|
||||||
@ -4973,100 +5045,44 @@ import numpy as np
|
|||||||
|
|
||||||
<<jastrow_data>>
|
<<jastrow_data>>
|
||||||
|
|
||||||
|
<<helper_funcs>>
|
||||||
|
|
||||||
kappa = 1.0
|
kappa = 1.0
|
||||||
|
|
||||||
elec_coord = np.array(elec_coord)[0]
|
factor_een = 0.0
|
||||||
nucl_coord = np.array(nucl_coord)
|
|
||||||
elnuc_dist = np.zeros(shape=(elec_num, nucl_num),dtype=float)
|
|
||||||
for i in range(elec_num):
|
|
||||||
for j in range(nucl_num):
|
|
||||||
elnuc_dist[i, j] = np.linalg.norm(elec_coord[i] - nucl_coord[:,j])
|
|
||||||
|
|
||||||
elnuc_dist_deriv_e = np.zeros(shape=(4, elec_num, nucl_num),dtype=float)
|
for n in range(0, dim_cord_vect):
|
||||||
for a in range(nucl_num):
|
l = lkpm_of_cindex[0,n]
|
||||||
for i in range(elec_num):
|
k = lkpm_of_cindex[1,n]
|
||||||
rij_inv = 1.0 / elnuc_dist[i, a]
|
p = lkpm_of_cindex[2,n]
|
||||||
for ii in range(3):
|
m = lkpm_of_cindex[3,n]
|
||||||
elnuc_dist_deriv_e[ii, i, a] = (elec_coord[i][ii] - nucl_coord[ii][a]) * rij_inv
|
|
||||||
elnuc_dist_deriv_e[3, i, a] = 2.0 * rij_inv
|
|
||||||
|
|
||||||
en_distance_rescaled_deriv_e = np.zeros(shape=(4,elec_num,nucl_num),dtype=float)
|
|
||||||
for a in range(nucl_num):
|
|
||||||
for i in range(elec_num):
|
|
||||||
f = 1.0 - kappa * en_distance_rescaled[i][a]
|
|
||||||
for ii in range(4):
|
|
||||||
en_distance_rescaled_deriv_e[ii][i][a] = elnuc_dist_deriv_e[ii][i][a]
|
|
||||||
en_distance_rescaled_deriv_e[3][i][a] = en_distance_rescaled_deriv_e[3][i][a] + \
|
|
||||||
(-kappa * en_distance_rescaled_deriv_e[0][i][a] * en_distance_rescaled_deriv_e[0][i][a]) + \
|
|
||||||
(-kappa * en_distance_rescaled_deriv_e[1][i][a] * en_distance_rescaled_deriv_e[1][i][a]) + \
|
|
||||||
(-kappa * en_distance_rescaled_deriv_e[2][i][a] * en_distance_rescaled_deriv_e[2][i][a])
|
|
||||||
for ii in range(4):
|
|
||||||
en_distance_rescaled_deriv_e[ii][i][a] = en_distance_rescaled_deriv_e[ii][i][a] * f
|
|
||||||
|
|
||||||
third = 1.0 / 3.0
|
|
||||||
factor_en_deriv_e = np.zeros(shape=(4,elec_num),dtype=float)
|
|
||||||
dx = np.zeros(shape=(4),dtype=float)
|
|
||||||
pow_ser_g = np.zeros(shape=(3),dtype=float)
|
|
||||||
for a in range(nucl_num):
|
|
||||||
for i in range(elec_num):
|
|
||||||
x = en_distance_rescaled[i][a]
|
|
||||||
if abs(x) < 1e-18:
|
|
||||||
continue
|
|
||||||
pow_ser_g = np.zeros(shape=(3),dtype=float)
|
|
||||||
den = 1.0 + aord_vector[1][type_nucl_vector[a]-1] * x
|
|
||||||
invden = 1.0 / den
|
|
||||||
invden2 = invden * invden
|
|
||||||
invden3 = invden2 * invden
|
|
||||||
xinv = 1.0 / (x + 1.0E-18)
|
|
||||||
|
|
||||||
for ii in range(4):
|
|
||||||
dx[ii] = en_distance_rescaled_deriv_e[ii][i][a]
|
|
||||||
|
|
||||||
lap1 = 0.0
|
|
||||||
lap2 = 0.0
|
|
||||||
lap3 = 0.0
|
|
||||||
for ii in range(3):
|
|
||||||
x = en_distance_rescaled[i][a]
|
|
||||||
if x < 1e-18:
|
|
||||||
continue
|
|
||||||
for p in range(2,aord_num+1):
|
|
||||||
y = p * aord_vector[(p-1) + 1][type_nucl_vector[a]-1] * x
|
|
||||||
pow_ser_g[ii] = pow_ser_g[ii] + y * dx[ii]
|
|
||||||
lap1 = lap1 + (p - 1) * y * xinv * dx[ii] * dx[ii]
|
|
||||||
lap2 = lap2 + y
|
|
||||||
x = x * en_distance_rescaled[i][a]
|
|
||||||
|
|
||||||
lap3 = lap3 - 2.0 * aord_vector[1][type_nucl_vector[a]-1] * dx[ii] * dx[ii]
|
|
||||||
|
|
||||||
factor_en_deriv_e[ii][i] = factor_en_deriv_e[ii][i] + aord_vector[0][type_nucl_vector[a]-1] * \
|
|
||||||
dx[ii] * invden2 + pow_ser_g[ii]
|
|
||||||
|
|
||||||
ii = 3
|
|
||||||
lap2 = lap2 * dx[ii] * third
|
|
||||||
lap3 = lap3 + den * dx[ii]
|
|
||||||
lap3 = lap3 * (aord_vector[0][type_nucl_vector[a]-1] * invden3)
|
|
||||||
factor_en_deriv_e[ii][i] = factor_en_deriv_e[ii][i] + lap1 + lap2 + lap3
|
|
||||||
|
|
||||||
print("factor_en_deriv_e[0][0]:",factor_en_deriv_e[0][0])
|
|
||||||
print("factor_en_deriv_e[1][0]:",factor_en_deriv_e[1][0])
|
|
||||||
print("factor_en_deriv_e[2][0]:",factor_en_deriv_e[2][0])
|
|
||||||
print("factor_en_deriv_e[3][0]:",factor_en_deriv_e[3][0])
|
|
||||||
|
|
||||||
|
for a in range(0, nucl_num):
|
||||||
|
accu2 = 0.0
|
||||||
|
cn = cord_vector_full[a][n]
|
||||||
|
for j in range(0, elec_num):
|
||||||
|
accu = 0.0
|
||||||
|
for i in range(0, elec_num):
|
||||||
|
accu = accu + een_rescaled_e[i,j,k] * \
|
||||||
|
een_rescaled_n[a,i,m]
|
||||||
|
accu2 = accu2 + accu * een_rescaled_n[a,j,m+l]
|
||||||
|
factor_een = factor_een + accu2 * cn
|
||||||
|
|
||||||
|
print("factor_een:",factor_een)
|
||||||
#+end_src
|
#+end_src
|
||||||
|
|
||||||
#+RESULTS:
|
#+RESULTS:
|
||||||
: factor_en_deriv_e[0][0]: 0.11609919541763383
|
: factor_een: -0.37407972141304213
|
||||||
: factor_en_deriv_e[1][0]: -0.23301394780804574
|
|
||||||
: factor_en_deriv_e[2][0]: 0.17548337641865783
|
|
||||||
: factor_en_deriv_e[3][0]: -0.9667363412285741
|
|
||||||
|
|
||||||
|
|
||||||
#+begin_src c :tangle (eval c_test)
|
#+begin_src c :tangle (eval c_test)
|
||||||
/* Check if Jastrow is properly initialized */
|
/* Check if Jastrow is properly initialized */
|
||||||
//assert(qmckl_jastrow_provided(context));
|
assert(qmckl_jastrow_provided(context));
|
||||||
//
|
|
||||||
|
|
||||||
|
double factor_een[walk_num];
|
||||||
|
rc = qmckl_get_jastrow_factor_een(context, &(factor_een[0]));
|
||||||
|
|
||||||
|
assert(fabs(factor_een[0] + 0.37407972141304213) < 1e-12);
|
||||||
#+end_src
|
#+end_src
|
||||||
|
|
||||||
** Electron-electron-nucleus Jastrow \(f_{een}\) derivative
|
** Electron-electron-nucleus Jastrow \(f_{een}\) derivative
|
||||||
@ -5221,12 +5237,12 @@ integer function qmckl_compute_factor_een_deriv_e_f(context, walk_num, elec_num,
|
|||||||
implicit none
|
implicit none
|
||||||
integer(qmckl_context), intent(in) :: context
|
integer(qmckl_context), intent(in) :: context
|
||||||
integer*8 , intent(in) :: walk_num, elec_num, cord_num, nucl_num, dim_cord_vect
|
integer*8 , intent(in) :: walk_num, elec_num, cord_num, nucl_num, dim_cord_vect
|
||||||
integer*8 , intent(in) :: lkpm_combined_index(4,dim_cord_vect)
|
integer*8 , intent(in) :: lkpm_combined_index(dim_cord_vect, 4)
|
||||||
double precision , intent(in) :: cord_vect_full(dim_cord_vect, nucl_num)
|
double precision , intent(in) :: cord_vect_full(nucl_num, dim_cord_vect)
|
||||||
double precision , intent(in) :: een_rescaled_e(walk_num, elec_num, elec_num, 0:cord_num)
|
double precision , intent(in) :: een_rescaled_e(0:cord_num, elec_num, elec_num, walk_num)
|
||||||
double precision , intent(in) :: een_rescaled_n(walk_num, elec_num, nucl_num, 0:cord_num)
|
double precision , intent(in) :: een_rescaled_n(0:cord_num, nucl_num, elec_num, walk_num)
|
||||||
double precision , intent(in) :: een_rescaled_e_deriv_e(walk_num, elec_num, 4, elec_num, 0:cord_num)
|
double precision , intent(in) :: een_rescaled_e_deriv_e(0:cord_num, elec_num, 4, elec_num, walk_num)
|
||||||
double precision , intent(in) :: een_rescaled_n_deriv_e(walk_num, elec_num, 4, nucl_num, 0:cord_num)
|
double precision , intent(in) :: een_rescaled_n_deriv_e(0:cord_num, nucl_num, 4, elec_num, walk_num)
|
||||||
double precision , intent(out) :: factor_een_deriv_e(elec_num, 4, walk_num)
|
double precision , intent(out) :: factor_een_deriv_e(elec_num, 4, walk_num)
|
||||||
|
|
||||||
integer*8 :: i, a, j, l, k, p, m, n, nw
|
integer*8 :: i, a, j, l, k, p, m, n, nw
|
||||||
@ -5264,41 +5280,41 @@ integer function qmckl_compute_factor_een_deriv_e_f(context, walk_num, elec_num,
|
|||||||
|
|
||||||
do nw =1, walk_num
|
do nw =1, walk_num
|
||||||
do n = 1, dim_cord_vect
|
do n = 1, dim_cord_vect
|
||||||
l = lkpm_combined_index(1, n)
|
l = lkpm_combined_index(n, 1)
|
||||||
k = lkpm_combined_index(2, n)
|
k = lkpm_combined_index(n, 2)
|
||||||
p = lkpm_combined_index(3, n)
|
p = lkpm_combined_index(n, 3)
|
||||||
m = lkpm_combined_index(4, n)
|
m = lkpm_combined_index(n, 4)
|
||||||
|
|
||||||
do a = 1, nucl_num
|
do a = 1, nucl_num
|
||||||
cn = cord_vect_full(n, a)
|
cn = cord_vect_full(a, n)
|
||||||
do j = 1, elec_num
|
do j = 1, elec_num
|
||||||
accu = 0.0d0
|
accu = 0.0d0
|
||||||
accu2 = 0.0d0
|
accu2 = 0.0d0
|
||||||
daccu = 0.0d0
|
daccu = 0.0d0
|
||||||
daccu2 = 0.0d0
|
daccu2 = 0.0d0
|
||||||
do i = 1, elec_num
|
do i = 1, elec_num
|
||||||
accu = accu + een_rescaled_e(nw, i, j, k) * &
|
accu = accu + een_rescaled_e(k, i, j, nw) * &
|
||||||
een_rescaled_n(nw, i, a, m)
|
een_rescaled_n(m, a, i, nw)
|
||||||
accu2 = accu2 + een_rescaled_e(nw, i, j, k) * &
|
accu2 = accu2 + een_rescaled_e(k, i, j, nw) * &
|
||||||
een_rescaled_n(nw, i, a, m + l)
|
een_rescaled_n(m + l, a, i, nw)
|
||||||
daccu(1:4) = daccu(1:4) + een_rescaled_e_deriv_e(nw, j, 1:4, i, k) * &
|
daccu(1:4) = daccu(1:4) + een_rescaled_e_deriv_e(k, j, 1:4, i, nw) * &
|
||||||
een_rescaled_n(nw, i, a, m)
|
een_rescaled_n(m, a, i, nw)
|
||||||
daccu2(1:4) = daccu2(1:4) + een_rescaled_e_deriv_e(nw, j, 1:4, i, k) * &
|
daccu2(1:4) = daccu2(1:4) + een_rescaled_e_deriv_e(k, j, 1:4, i, nw) * &
|
||||||
een_rescaled_n(nw, i, a, m + l)
|
een_rescaled_n(m + l, a, i, nw)
|
||||||
end do
|
end do
|
||||||
factor_een_deriv_e(j, 1:4, nw) = factor_een_deriv_e(j, 1:4, nw) + &
|
factor_een_deriv_e(j, 1:4, nw) = factor_een_deriv_e(j, 1:4, nw) + &
|
||||||
(accu * een_rescaled_n_deriv_e(nw, j, 1:4, a, m + l) &
|
(accu * een_rescaled_n_deriv_e(m + l, a, 1:4, j, nw) &
|
||||||
+ daccu(1:4) * een_rescaled_n(nw, j, a, m + l) &
|
+ daccu(1:4) * een_rescaled_n(m + l, a, j, nw) &
|
||||||
+ daccu2(1:4) * een_rescaled_n(nw, j, a, m) &
|
+ daccu2(1:4) * een_rescaled_n(m, a, j, nw) &
|
||||||
+ accu2 * een_rescaled_n_deriv_e(nw, j, 1:4, a, m)) * cn
|
+ accu2 * een_rescaled_n_deriv_e(m, a, 1:4, j, nw)) * cn
|
||||||
|
|
||||||
factor_een_deriv_e(j, 4, nw) = factor_een_deriv_e(j, 4, nw) + 2.0d0 * ( &
|
factor_een_deriv_e(j, 4, nw) = factor_een_deriv_e(j, 4, nw) + 2.0d0 * ( &
|
||||||
daccu (1) * een_rescaled_n_deriv_e(nw, j, 1, a, m + l) + &
|
daccu (1) * een_rescaled_n_deriv_e(m + l, a, 1, j, nw) + &
|
||||||
daccu (2) * een_rescaled_n_deriv_e(nw, j, 2, a, m + l) + &
|
daccu (2) * een_rescaled_n_deriv_e(m + l, a, 2, j, nw) + &
|
||||||
daccu (3) * een_rescaled_n_deriv_e(nw, j, 3, a, m + l) + &
|
daccu (3) * een_rescaled_n_deriv_e(m + l, a, 3, j, nw) + &
|
||||||
daccu2(1) * een_rescaled_n_deriv_e(nw, j, 1, a, m ) + &
|
daccu2(1) * een_rescaled_n_deriv_e(m, a, 1, j, nw ) + &
|
||||||
daccu2(2) * een_rescaled_n_deriv_e(nw, j, 2, a, m ) + &
|
daccu2(2) * een_rescaled_n_deriv_e(m, a, 2, j, nw ) + &
|
||||||
daccu2(3) * een_rescaled_n_deriv_e(nw, j, 3, a, m ) ) * cn
|
daccu2(3) * een_rescaled_n_deriv_e(m, a, 3, j, nw ) ) * cn
|
||||||
|
|
||||||
end do
|
end do
|
||||||
end do
|
end do
|
||||||
@ -5391,98 +5407,60 @@ import numpy as np
|
|||||||
|
|
||||||
<<jastrow_data>>
|
<<jastrow_data>>
|
||||||
|
|
||||||
|
<<een_e_deriv_e>>
|
||||||
|
|
||||||
|
<<helper_funcs>>
|
||||||
|
|
||||||
kappa = 1.0
|
kappa = 1.0
|
||||||
|
|
||||||
elec_coord = np.array(elec_coord)[0]
|
factor_een = 0.0
|
||||||
nucl_coord = np.array(nucl_coord)
|
|
||||||
elnuc_dist = np.zeros(shape=(elec_num, nucl_num),dtype=float)
|
|
||||||
for i in range(elec_num):
|
|
||||||
for j in range(nucl_num):
|
|
||||||
elnuc_dist[i, j] = np.linalg.norm(elec_coord[i] - nucl_coord[:,j])
|
|
||||||
|
|
||||||
elnuc_dist_deriv_e = np.zeros(shape=(4, elec_num, nucl_num),dtype=float)
|
daccu = np.zeros(4, dtype=float)
|
||||||
for a in range(nucl_num):
|
daccu2 = np.zeros(4, dtype=float)
|
||||||
for i in range(elec_num):
|
een_rescaled_e_deriv_e_t = een_rescaled_e_deriv_e.T
|
||||||
rij_inv = 1.0 / elnuc_dist[i, a]
|
print(een_rescaled_e_deriv_e_t.shape)
|
||||||
for ii in range(3):
|
for n in range(0, dim_cord_vect):
|
||||||
elnuc_dist_deriv_e[ii, i, a] = (elec_coord[i][ii] - nucl_coord[ii][a]) * rij_inv
|
l = lkpm_of_cindex[0,n]
|
||||||
elnuc_dist_deriv_e[3, i, a] = 2.0 * rij_inv
|
k = lkpm_of_cindex[1,n]
|
||||||
|
p = lkpm_of_cindex[2,n]
|
||||||
|
m = lkpm_of_cindex[3,n]
|
||||||
|
|
||||||
en_distance_rescaled_deriv_e = np.zeros(shape=(4,elec_num,nucl_num),dtype=float)
|
for a in range(0, nucl_num):
|
||||||
for a in range(nucl_num):
|
cn = cord_vector_full[a][n]
|
||||||
for i in range(elec_num):
|
for j in range(0, elec_num):
|
||||||
f = 1.0 - kappa * en_distance_rescaled[i][a]
|
accu = 0.0
|
||||||
for ii in range(4):
|
accu2 = 0.0
|
||||||
en_distance_rescaled_deriv_e[ii][i][a] = elnuc_dist_deriv_e[ii][i][a]
|
daccu = 0.0
|
||||||
en_distance_rescaled_deriv_e[3][i][a] = en_distance_rescaled_deriv_e[3][i][a] + \
|
daccu2 = 0.0
|
||||||
(-kappa * en_distance_rescaled_deriv_e[0][i][a] * en_distance_rescaled_deriv_e[0][i][a]) + \
|
for i in range(0, elec_num):
|
||||||
(-kappa * en_distance_rescaled_deriv_e[1][i][a] * en_distance_rescaled_deriv_e[1][i][a]) + \
|
accu = accu + een_rescaled_e[i,j,k] * \
|
||||||
(-kappa * en_distance_rescaled_deriv_e[2][i][a] * en_distance_rescaled_deriv_e[2][i][a])
|
een_rescaled_n[a,i,m]
|
||||||
for ii in range(4):
|
accu2 = accu2 + een_rescaled_e[i,j,k] * \
|
||||||
en_distance_rescaled_deriv_e[ii][i][a] = en_distance_rescaled_deriv_e[ii][i][a] * f
|
een_rescaled_n[a,i,m+l]
|
||||||
|
# daccu[0:4] = daccu[0:4] + een_rescaled_e_deriv_e_t[k,j,0:4,i,k] * \
|
||||||
third = 1.0 / 3.0
|
# een_rescaled_n[a,i,m]
|
||||||
factor_en_deriv_e = np.zeros(shape=(4,elec_num),dtype=float)
|
# daccu[0:4] = daccu[0:4] + een_rescaled_e_deriv_e_t[k,j,0:4,i,k] * \
|
||||||
dx = np.zeros(shape=(4),dtype=float)
|
# een_rescaled_n[a,i,m]
|
||||||
pow_ser_g = np.zeros(shape=(3),dtype=float)
|
accu2 = accu2 + accu * een_rescaled_n[a,j,m+l]
|
||||||
for a in range(nucl_num):
|
# factor_een = factor_een + accu2 * cn
|
||||||
for i in range(elec_num):
|
|
||||||
x = en_distance_rescaled[i][a]
|
|
||||||
if abs(x) < 1e-18:
|
|
||||||
continue
|
|
||||||
pow_ser_g = np.zeros(shape=(3),dtype=float)
|
|
||||||
den = 1.0 + aord_vector[1][type_nucl_vector[a]-1] * x
|
|
||||||
invden = 1.0 / den
|
|
||||||
invden2 = invden * invden
|
|
||||||
invden3 = invden2 * invden
|
|
||||||
xinv = 1.0 / (x + 1.0E-18)
|
|
||||||
|
|
||||||
for ii in range(4):
|
|
||||||
dx[ii] = en_distance_rescaled_deriv_e[ii][i][a]
|
|
||||||
|
|
||||||
lap1 = 0.0
|
|
||||||
lap2 = 0.0
|
|
||||||
lap3 = 0.0
|
|
||||||
for ii in range(3):
|
|
||||||
x = en_distance_rescaled[i][a]
|
|
||||||
if x < 1e-18:
|
|
||||||
continue
|
|
||||||
for p in range(2,aord_num+1):
|
|
||||||
y = p * aord_vector[(p-1) + 1][type_nucl_vector[a]-1] * x
|
|
||||||
pow_ser_g[ii] = pow_ser_g[ii] + y * dx[ii]
|
|
||||||
lap1 = lap1 + (p - 1) * y * xinv * dx[ii] * dx[ii]
|
|
||||||
lap2 = lap2 + y
|
|
||||||
x = x * en_distance_rescaled[i][a]
|
|
||||||
|
|
||||||
lap3 = lap3 - 2.0 * aord_vector[1][type_nucl_vector[a]-1] * dx[ii] * dx[ii]
|
|
||||||
|
|
||||||
factor_en_deriv_e[ii][i] = factor_en_deriv_e[ii][i] + aord_vector[0][type_nucl_vector[a]-1] * \
|
|
||||||
dx[ii] * invden2 + pow_ser_g[ii]
|
|
||||||
|
|
||||||
ii = 3
|
|
||||||
lap2 = lap2 * dx[ii] * third
|
|
||||||
lap3 = lap3 + den * dx[ii]
|
|
||||||
lap3 = lap3 * (aord_vector[0][type_nucl_vector[a]-1] * invden3)
|
|
||||||
factor_en_deriv_e[ii][i] = factor_en_deriv_e[ii][i] + lap1 + lap2 + lap3
|
|
||||||
|
|
||||||
print("factor_en_deriv_e[0][0]:",factor_en_deriv_e[0][0])
|
|
||||||
print("factor_en_deriv_e[1][0]:",factor_en_deriv_e[1][0])
|
|
||||||
print("factor_en_deriv_e[2][0]:",factor_en_deriv_e[2][0])
|
|
||||||
print("factor_en_deriv_e[3][0]:",factor_en_deriv_e[3][0])
|
|
||||||
|
|
||||||
|
print("factor_een:",factor_een)
|
||||||
|
|
||||||
#+end_src
|
#+end_src
|
||||||
|
|
||||||
#+RESULTS:
|
#+RESULTS:
|
||||||
: factor_en_deriv_e[0][0]: 0.11609919541763383
|
: (6, 10, 4, 10)
|
||||||
: factor_en_deriv_e[1][0]: -0.23301394780804574
|
: factor_een: 0.0
|
||||||
: factor_en_deriv_e[2][0]: 0.17548337641865783
|
|
||||||
: factor_en_deriv_e[3][0]: -0.9667363412285741
|
|
||||||
|
|
||||||
|
|
||||||
#+begin_src c :tangle (eval c_test)
|
#+begin_src c :tangle (eval c_test)
|
||||||
///* Check if Jastrow is properly initialized */
|
/* Check if Jastrow is properly initialized */
|
||||||
|
assert(qmckl_jastrow_provided(context));
|
||||||
|
|
||||||
|
double factor_een_deriv_e[walk_num][elec_num];
|
||||||
|
rc = qmckl_get_jastrow_factor_een_deriv_e(context, &(factor_een_deriv_e[0][0]));
|
||||||
|
|
||||||
|
assert(fabs(factor_een_deriv_e[0][0] + 0.0005481671107226865) < 1e-12);
|
||||||
#+end_src
|
#+end_src
|
||||||
|
|
||||||
* End of files :noexport:
|
* End of files :noexport:
|
||||||
|
880
org/qmckl_mo.org
Normal file
880
org/qmckl_mo.org
Normal file
@ -0,0 +1,880 @@
|
|||||||
|
#+TITLE: Molecular Orbitals
|
||||||
|
#+SETUPFILE: ../tools/theme.setup
|
||||||
|
#+INCLUDE: ../tools/lib.org
|
||||||
|
|
||||||
|
The molecular orbitals (MOs) are defined in the basis of AOs along with a AO to MO
|
||||||
|
coefficient matrix \[C\]. Using these coefficients (e.g. from Hartree Fock SCF method)
|
||||||
|
the MOs are defined as follows:
|
||||||
|
|
||||||
|
\[
|
||||||
|
\phi_i(\mathbf{r}) = C_i * \chi_i (\mathbf{r})
|
||||||
|
\]
|
||||||
|
|
||||||
|
|
||||||
|
In this section we demonstrate how to use the QMCkl specific DGEMM
|
||||||
|
function to calculate the MOs.
|
||||||
|
|
||||||
|
|
||||||
|
* Headers :noexport:
|
||||||
|
#+begin_src elisp :noexport :results none
|
||||||
|
(org-babel-lob-ingest "../tools/lib.org")
|
||||||
|
#+end_src
|
||||||
|
|
||||||
|
|
||||||
|
#+begin_src c :tangle (eval h_private_type)
|
||||||
|
#ifndef QMCKL_MO_HPT
|
||||||
|
#define QMCKL_MO_HPT
|
||||||
|
|
||||||
|
#include <stdbool.h>
|
||||||
|
#+end_src
|
||||||
|
|
||||||
|
#+begin_src c :tangle (eval c_test) :noweb yes
|
||||||
|
#include "qmckl.h"
|
||||||
|
#include "assert.h"
|
||||||
|
#ifdef HAVE_CONFIG_H
|
||||||
|
#include "config.h"
|
||||||
|
#endif
|
||||||
|
|
||||||
|
#include <stdio.h>
|
||||||
|
#include <math.h>
|
||||||
|
#include "chbrclf.h"
|
||||||
|
#include "qmckl_ao_private_func.h"
|
||||||
|
#include "qmckl_mo_private_func.h"
|
||||||
|
|
||||||
|
int main() {
|
||||||
|
qmckl_context context;
|
||||||
|
context = qmckl_context_create();
|
||||||
|
|
||||||
|
qmckl_exit_code rc;
|
||||||
|
#+end_src
|
||||||
|
|
||||||
|
#+begin_src c :tangle (eval c)
|
||||||
|
#ifdef HAVE_CONFIG_H
|
||||||
|
#include "config.h"
|
||||||
|
#endif
|
||||||
|
|
||||||
|
#ifdef HAVE_STDINT_H
|
||||||
|
#include <stdint.h>
|
||||||
|
#elif HAVE_INTTYPES_H
|
||||||
|
#include <inttypes.h>
|
||||||
|
#endif
|
||||||
|
|
||||||
|
#include <stdlib.h>
|
||||||
|
#include <string.h>
|
||||||
|
#include <stdbool.h>
|
||||||
|
#include <assert.h>
|
||||||
|
|
||||||
|
#include "qmckl.h"
|
||||||
|
#include "qmckl_context_private_type.h"
|
||||||
|
#include "qmckl_memory_private_type.h"
|
||||||
|
#include "qmckl_memory_private_func.h"
|
||||||
|
#include "qmckl_ao_private_type.h"
|
||||||
|
#include "qmckl_ao_private_func.h"
|
||||||
|
#include "qmckl_mo_private_type.h"
|
||||||
|
#include "qmckl_mo_private_func.h"
|
||||||
|
#+end_src
|
||||||
|
|
||||||
|
* Context
|
||||||
|
|
||||||
|
The following arrays are stored in the context:
|
||||||
|
|
||||||
|
|---------------------+--------------------+--------------------------------------------------------------|
|
||||||
|
| ~type~ | | Gaussian (~'G'~) or Slater (~'S'~) |
|
||||||
|
| ~mo_num~ | | Number of MOs |
|
||||||
|
| ~coefficient~ | ~[mo_num, ao_num]~ | Orbital coefficients |
|
||||||
|
|
||||||
|
Computed data:
|
||||||
|
|
||||||
|
|---------------+-----------------------------------+----------------------------------------------------------------------------------------|
|
||||||
|
|---------------+-----------------------------------+----------------------------------------------------------------------------------------|
|
||||||
|
| ~mo_vgl~ | ~[5][walk_num][elec_num][mo_num]~ | Value, gradients, Laplacian of the MOs at electron positions |
|
||||||
|
| ~mo_vgl_date~ | ~uint64_t~ | Late modification date of Value, gradients, Laplacian of the MOs at electron positions |
|
||||||
|
|---------------+-----------------------------------+----------------------------------------------------------------------------------------|
|
||||||
|
|
||||||
|
** Data structure
|
||||||
|
|
||||||
|
#+begin_src c :comments org :tangle (eval h_private_type)
|
||||||
|
typedef struct qmckl_mo_basis_struct {
|
||||||
|
int64_t mo_num;
|
||||||
|
double * coefficient;
|
||||||
|
|
||||||
|
double * mo_vgl;
|
||||||
|
int64_t mo_vgl_date;
|
||||||
|
|
||||||
|
int32_t uninitialized;
|
||||||
|
bool provided;
|
||||||
|
char type;
|
||||||
|
} qmckl_mo_basis_struct;
|
||||||
|
#+end_src
|
||||||
|
|
||||||
|
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.
|
||||||
|
|
||||||
|
** Access functions
|
||||||
|
|
||||||
|
#+begin_src c :comments org :tangle (eval h_private_func) :exports none
|
||||||
|
char qmckl_get_mo_basis_type (const qmckl_context context);
|
||||||
|
int64_t qmckl_get_mo_basis_mo_num (const qmckl_context context);
|
||||||
|
double* qmckl_get_mo_basis_coefficient (const qmckl_context context);
|
||||||
|
#+end_src
|
||||||
|
|
||||||
|
When all the data for the AOs have been provided, the following
|
||||||
|
function returns ~true~.
|
||||||
|
|
||||||
|
#+begin_src c :comments org :tangle (eval h_func)
|
||||||
|
bool qmckl_mo_basis_provided (const qmckl_context context);
|
||||||
|
#+end_src
|
||||||
|
|
||||||
|
#+NAME:post
|
||||||
|
#+begin_src c :exports none
|
||||||
|
if ( (ctx->mo_basis.uninitialized & mask) != 0) {
|
||||||
|
return NULL;
|
||||||
|
}
|
||||||
|
#+end_src
|
||||||
|
|
||||||
|
|
||||||
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
||||||
|
char qmckl_get_mo_basis_type (const qmckl_context context) {
|
||||||
|
|
||||||
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
||||||
|
return (char) 0;
|
||||||
|
}
|
||||||
|
|
||||||
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
||||||
|
assert (ctx != NULL);
|
||||||
|
|
||||||
|
int32_t mask = 1;
|
||||||
|
|
||||||
|
if ( (ctx->mo_basis.uninitialized & mask) != 0) {
|
||||||
|
return (char) 0;
|
||||||
|
}
|
||||||
|
|
||||||
|
assert (ctx->mo_basis.type != (char) 0);
|
||||||
|
return ctx->mo_basis.type;
|
||||||
|
}
|
||||||
|
|
||||||
|
int64_t qmckl_get_mo_basis_mo_num (const qmckl_context context) {
|
||||||
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
||||||
|
return (int64_t) 0;
|
||||||
|
}
|
||||||
|
|
||||||
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
||||||
|
assert (ctx != NULL);
|
||||||
|
|
||||||
|
int32_t mask = 1 << 1;
|
||||||
|
|
||||||
|
if ( (ctx->mo_basis.uninitialized & mask) != 0) {
|
||||||
|
return (int64_t) 0;
|
||||||
|
}
|
||||||
|
|
||||||
|
assert (ctx->mo_basis.mo_num > (int64_t) 0);
|
||||||
|
return ctx->mo_basis.mo_num;
|
||||||
|
}
|
||||||
|
|
||||||
|
bool qmckl_mo_basis_provided(const 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);
|
||||||
|
|
||||||
|
return ctx->mo_basis.provided;
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
#+end_src
|
||||||
|
|
||||||
|
** Initialization functions
|
||||||
|
|
||||||
|
To set the basis set, all the following functions need to be
|
||||||
|
called.
|
||||||
|
|
||||||
|
#+begin_src c :comments org :tangle (eval h_func)
|
||||||
|
qmckl_exit_code qmckl_set_mo_basis_type (qmckl_context context, const char t);
|
||||||
|
qmckl_exit_code qmckl_set_mo_basis_mo_num (qmckl_context context, const int64_t mo_num);
|
||||||
|
qmckl_exit_code qmckl_set_mo_basis_coefficient (qmckl_context context, const double * coefficient);
|
||||||
|
#+end_src
|
||||||
|
|
||||||
|
#+NAME:pre2
|
||||||
|
#+begin_src c :exports none
|
||||||
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
||||||
|
return QMCKL_NULL_CONTEXT;
|
||||||
|
}
|
||||||
|
|
||||||
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
||||||
|
#+end_src
|
||||||
|
|
||||||
|
#+NAME:post2
|
||||||
|
#+begin_src c :exports none
|
||||||
|
ctx->mo_basis.uninitialized &= ~mask;
|
||||||
|
ctx->mo_basis.provided = (ctx->mo_basis.uninitialized == 0);
|
||||||
|
if (ctx->mo_basis.provided) {
|
||||||
|
qmckl_exit_code rc_ = qmckl_finalize_basis(context);
|
||||||
|
if (rc_ != QMCKL_SUCCESS) return rc_;
|
||||||
|
}
|
||||||
|
|
||||||
|
return QMCKL_SUCCESS;
|
||||||
|
#+end_src
|
||||||
|
|
||||||
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
||||||
|
qmckl_exit_code qmckl_set_mo_basis_type(qmckl_context context, const char t) {
|
||||||
|
<<pre2>>
|
||||||
|
|
||||||
|
if (t != 'G' && t != 'S') {
|
||||||
|
return qmckl_failwith( context,
|
||||||
|
QMCKL_INVALID_ARG_2,
|
||||||
|
"qmckl_set_mo_basis_type",
|
||||||
|
NULL);
|
||||||
|
}
|
||||||
|
|
||||||
|
int32_t mask = 1;
|
||||||
|
ctx->mo_basis.type = t;
|
||||||
|
|
||||||
|
<<post2>>
|
||||||
|
}
|
||||||
|
|
||||||
|
qmckl_exit_code qmckl_set_mo_basis_mo_num(qmckl_context context, const int64_t mo_num) {
|
||||||
|
<<pre2>>
|
||||||
|
|
||||||
|
if (mo_num <= 0) {
|
||||||
|
return qmckl_failwith( context,
|
||||||
|
QMCKL_INVALID_ARG_2,
|
||||||
|
"qmckl_set_mo_basis_mo_num",
|
||||||
|
"mo_num <= 0");
|
||||||
|
}
|
||||||
|
|
||||||
|
int32_t mask = 1 << 1;
|
||||||
|
ctx->mo_basis.mo_num = mo_num;
|
||||||
|
|
||||||
|
<<post2>>
|
||||||
|
}
|
||||||
|
|
||||||
|
qmckl_exit_code qmckl_set_mo_basis_coefficient(qmckl_context context, const double* coefficient) {
|
||||||
|
<<pre2>>
|
||||||
|
|
||||||
|
int32_t mask = 1 << 2;
|
||||||
|
|
||||||
|
if (ctx->mo_basis.coefficient != NULL) {
|
||||||
|
qmckl_exit_code rc = qmckl_free(context, ctx->mo_basis.coefficient);
|
||||||
|
if (rc != QMCKL_SUCCESS) {
|
||||||
|
return qmckl_failwith( context, rc,
|
||||||
|
"qmckl_set_mo_basis_coefficient",
|
||||||
|
NULL);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
|
||||||
|
mem_info.size = ctx->ao_basis.ao_num * ctx->mo_basis.mo_num * sizeof(double);
|
||||||
|
double* new_array = (double*) qmckl_malloc(context, mem_info);
|
||||||
|
if (new_array == NULL) {
|
||||||
|
return qmckl_failwith( context,
|
||||||
|
QMCKL_ALLOCATION_FAILED,
|
||||||
|
"qmckl_set_mo_basis_coefficient",
|
||||||
|
NULL);
|
||||||
|
}
|
||||||
|
|
||||||
|
memcpy(new_array, coefficient, mem_info.size);
|
||||||
|
|
||||||
|
ctx->mo_basis.coefficient = new_array;
|
||||||
|
|
||||||
|
<<post2>>
|
||||||
|
}
|
||||||
|
|
||||||
|
#+end_src
|
||||||
|
|
||||||
|
When the basis set is completely entered, other data structures are
|
||||||
|
computed to accelerate the calculations.
|
||||||
|
|
||||||
|
* Computation
|
||||||
|
|
||||||
|
** Computation of MOs
|
||||||
|
|
||||||
|
*** Get
|
||||||
|
|
||||||
|
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
|
||||||
|
qmckl_exit_code qmckl_get_mo_basis_vgl(qmckl_context context, double* const mo_vgl);
|
||||||
|
#+end_src
|
||||||
|
|
||||||
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
||||||
|
qmckl_exit_code qmckl_get_mo_basis_vgl(qmckl_context context, double* const mo_vgl) {
|
||||||
|
|
||||||
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
||||||
|
return QMCKL_NULL_CONTEXT;
|
||||||
|
}
|
||||||
|
|
||||||
|
qmckl_exit_code rc;
|
||||||
|
|
||||||
|
rc = qmckl_provide_ao_vgl(context);
|
||||||
|
if (rc != QMCKL_SUCCESS) return rc;
|
||||||
|
|
||||||
|
rc = qmckl_provide_mo_vgl(context);
|
||||||
|
if (rc != QMCKL_SUCCESS) return rc;
|
||||||
|
|
||||||
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
||||||
|
assert (ctx != NULL);
|
||||||
|
|
||||||
|
size_t sze = 5 * ctx->electron.num * ctx->mo_basis.mo_num * ctx->electron.walk_num;
|
||||||
|
memcpy(mo_vgl, ctx->mo_basis.mo_vgl, sze * sizeof(double));
|
||||||
|
|
||||||
|
return QMCKL_SUCCESS;
|
||||||
|
}
|
||||||
|
#+end_src
|
||||||
|
|
||||||
|
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
|
||||||
|
interface
|
||||||
|
integer(c_int32_t) function qmckl_get_mo_basis_vgl (context, mo_vgl) &
|
||||||
|
bind(C)
|
||||||
|
use, intrinsic :: iso_c_binding
|
||||||
|
import
|
||||||
|
implicit none
|
||||||
|
|
||||||
|
integer (c_int64_t) , intent(in) , value :: context
|
||||||
|
double precision, intent(out) :: mo_vgl(*)
|
||||||
|
end function
|
||||||
|
end interface
|
||||||
|
#+end_src
|
||||||
|
|
||||||
|
*** Provide
|
||||||
|
|
||||||
|
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
|
||||||
|
qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context);
|
||||||
|
#+end_src
|
||||||
|
|
||||||
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
||||||
|
qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context)
|
||||||
|
{
|
||||||
|
|
||||||
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
||||||
|
return QMCKL_NULL_CONTEXT;
|
||||||
|
}
|
||||||
|
|
||||||
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
||||||
|
assert (ctx != NULL);
|
||||||
|
|
||||||
|
if (!ctx->ao_basis.provided) {
|
||||||
|
return qmckl_failwith( context,
|
||||||
|
QMCKL_NOT_PROVIDED,
|
||||||
|
"qmckl_ao_basis",
|
||||||
|
NULL);
|
||||||
|
}
|
||||||
|
|
||||||
|
if(!(ctx->electron.provided)) {
|
||||||
|
return qmckl_failwith( context,
|
||||||
|
QMCKL_NOT_PROVIDED,
|
||||||
|
"qmckl_electron",
|
||||||
|
NULL);
|
||||||
|
}
|
||||||
|
|
||||||
|
if (!ctx->mo_basis.provided) {
|
||||||
|
return qmckl_failwith( context,
|
||||||
|
QMCKL_NOT_PROVIDED,
|
||||||
|
"qmckl_mo_basis",
|
||||||
|
NULL);
|
||||||
|
}
|
||||||
|
|
||||||
|
/* Compute if necessary */
|
||||||
|
if (ctx->electron.coord_new_date > ctx->mo_basis.mo_vgl_date) {
|
||||||
|
|
||||||
|
/* Allocate array */
|
||||||
|
if (ctx->mo_basis.mo_vgl == NULL) {
|
||||||
|
|
||||||
|
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
|
||||||
|
mem_info.size = 5 * ctx->electron.num * ctx->mo_basis.mo_num *
|
||||||
|
ctx->electron.walk_num * sizeof(double);
|
||||||
|
double* mo_vgl = (double*) qmckl_malloc(context, mem_info);
|
||||||
|
|
||||||
|
if (mo_vgl == NULL) {
|
||||||
|
return qmckl_failwith( context,
|
||||||
|
QMCKL_ALLOCATION_FAILED,
|
||||||
|
"qmckl_mo_basis_vgl",
|
||||||
|
NULL);
|
||||||
|
}
|
||||||
|
ctx->mo_basis.mo_vgl = mo_vgl;
|
||||||
|
}
|
||||||
|
|
||||||
|
qmckl_exit_code rc;
|
||||||
|
if (ctx->mo_basis.type == 'G') {
|
||||||
|
rc = qmckl_compute_mo_basis_gaussian_vgl(context,
|
||||||
|
ctx->ao_basis.ao_num,
|
||||||
|
ctx->mo_basis.mo_num,
|
||||||
|
ctx->electron.num,
|
||||||
|
ctx->electron.walk_num,
|
||||||
|
ctx->mo_basis.coefficient,
|
||||||
|
ctx->ao_basis.ao_vgl,
|
||||||
|
ctx->mo_basis.mo_vgl);
|
||||||
|
} else {
|
||||||
|
return qmckl_failwith( context,
|
||||||
|
QMCKL_FAILURE,
|
||||||
|
"compute_mo_basis_vgl",
|
||||||
|
"Not yet implemented");
|
||||||
|
}
|
||||||
|
if (rc != QMCKL_SUCCESS) {
|
||||||
|
return rc;
|
||||||
|
}
|
||||||
|
|
||||||
|
ctx->mo_basis.mo_vgl_date = ctx->date;
|
||||||
|
}
|
||||||
|
|
||||||
|
return QMCKL_SUCCESS;
|
||||||
|
}
|
||||||
|
#+end_src
|
||||||
|
|
||||||
|
*** Compute
|
||||||
|
:PROPERTIES:
|
||||||
|
:Name: qmckl_compute_mo_basis_gaussian_vgl
|
||||||
|
:CRetType: qmckl_exit_code
|
||||||
|
:FRetType: qmckl_exit_code
|
||||||
|
:END:
|
||||||
|
|
||||||
|
#+NAME: qmckl_mo_basis_gaussian_vgl_args
|
||||||
|
| ~qmckl_context~ | ~context~ | in | Global state |
|
||||||
|
| ~int64_t~ | ~ao_num~ | in | Number of AOs |
|
||||||
|
| ~int64_t~ | ~mo_num~ | in | Number of MOs |
|
||||||
|
| ~int64_t~ | ~elec_num~ | in | Number of electrons |
|
||||||
|
| ~int64_t~ | ~walk_num~ | in | Number of walkers |
|
||||||
|
| ~double~ | ~coef_normalized[mo_num][ao_num]~ | in | AO to MO transformation matrix |
|
||||||
|
| ~double~ | ~ao_vgl[5][walk_num][elec_num][ao_num]~ | in | Value, gradients and Laplacian of the AOs |
|
||||||
|
| ~double~ | ~mo_vgl[5][walk_num][elec_num][mo_num]~ | out | Value, gradients and Laplacian of the MOs |
|
||||||
|
|
||||||
|
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
|
||||||
|
integer function qmckl_compute_mo_basis_gaussian_vgl_f(context, &
|
||||||
|
ao_num, mo_num, elec_num, walk_num, &
|
||||||
|
coef_normalized, ao_vgl, mo_vgl) &
|
||||||
|
result(info)
|
||||||
|
use qmckl
|
||||||
|
implicit none
|
||||||
|
integer(qmckl_context), intent(in) :: context
|
||||||
|
integer*8 , intent(in) :: ao_num, mo_num
|
||||||
|
integer*8 , intent(in) :: elec_num
|
||||||
|
integer*8 , intent(in) :: walk_num
|
||||||
|
double precision , intent(in) :: ao_vgl(ao_num,elec_num,walk_num,5)
|
||||||
|
double precision , intent(in) :: coef_normalized(ao_num,mo_num)
|
||||||
|
double precision , intent(out) :: mo_vgl(mo_num,elec_num,walk_num,5)
|
||||||
|
logical*8 :: TransA, TransB
|
||||||
|
double precision :: alpha, beta
|
||||||
|
integer :: info_qmckl_dgemm_value
|
||||||
|
integer :: info_qmckl_dgemm_Gx
|
||||||
|
integer :: info_qmckl_dgemm_Gy
|
||||||
|
integer :: info_qmckl_dgemm_Gz
|
||||||
|
integer :: info_qmckl_dgemm_lap
|
||||||
|
integer*8 :: M, N, K, LDA, LDB, LDC, i,j
|
||||||
|
|
||||||
|
integer*8 :: inucl, iprim, iwalk, ielec, ishell
|
||||||
|
double precision :: x, y, z, two_a, ar2, r2, v, cutoff
|
||||||
|
TransA = .False.
|
||||||
|
TransB = .False.
|
||||||
|
alpha = 1.0d0
|
||||||
|
beta = 0.0d0
|
||||||
|
|
||||||
|
info = QMCKL_SUCCESS
|
||||||
|
info_qmckl_dgemm_value = QMCKL_SUCCESS
|
||||||
|
info_qmckl_dgemm_Gx = QMCKL_SUCCESS
|
||||||
|
info_qmckl_dgemm_Gy = QMCKL_SUCCESS
|
||||||
|
info_qmckl_dgemm_Gz = QMCKL_SUCCESS
|
||||||
|
info_qmckl_dgemm_lap = QMCKL_SUCCESS
|
||||||
|
|
||||||
|
! Don't compute exponentials when the result will be almost zero.
|
||||||
|
! TODO : Use numerical precision here
|
||||||
|
cutoff = -dlog(1.d-15)
|
||||||
|
M = 1_8
|
||||||
|
N = mo_num * 1_8
|
||||||
|
K = ao_num * 1_8
|
||||||
|
LDA = M
|
||||||
|
LDB = K
|
||||||
|
LDC = M
|
||||||
|
|
||||||
|
do iwalk = 1, walk_num
|
||||||
|
do ielec = 1, elec_num
|
||||||
|
! Value
|
||||||
|
info_qmckl_dgemm_value = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, &
|
||||||
|
ao_vgl(:, ielec, iwalk, 1), LDA, &
|
||||||
|
coef_normalized(1:ao_num,1:mo_num),size(coef_normalized,1) * 1_8, &
|
||||||
|
beta, &
|
||||||
|
mo_vgl(:,ielec,iwalk,1),LDC)
|
||||||
|
! Grad_x
|
||||||
|
info_qmckl_dgemm_Gx = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, &
|
||||||
|
ao_vgl(:, ielec, iwalk, 2), LDA, &
|
||||||
|
coef_normalized(1:ao_num,1:mo_num),size(coef_normalized,1) * 1_8, &
|
||||||
|
beta, &
|
||||||
|
mo_vgl(:,ielec,iwalk,2),LDC)
|
||||||
|
! Grad_y
|
||||||
|
info_qmckl_dgemm_Gy = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, &
|
||||||
|
ao_vgl(:, ielec, iwalk, 3), LDA, &
|
||||||
|
coef_normalized(1:ao_num,1:mo_num),size(coef_normalized,1) * 1_8, &
|
||||||
|
beta, &
|
||||||
|
mo_vgl(:,ielec,iwalk,3),LDC)
|
||||||
|
! Grad_z
|
||||||
|
info_qmckl_dgemm_Gz = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, &
|
||||||
|
ao_vgl(:, ielec, iwalk, 4), LDA, &
|
||||||
|
coef_normalized(1:ao_num,1:mo_num),size(coef_normalized,1) * 1_8, &
|
||||||
|
beta, &
|
||||||
|
mo_vgl(:,ielec,iwalk,4),LDC)
|
||||||
|
! Lapl_z
|
||||||
|
info_qmckl_dgemm_lap = qmckl_dgemm(context, TransA, TransB, M, N, K, alpha, &
|
||||||
|
ao_vgl(:, ielec, iwalk, 5), LDA, &
|
||||||
|
coef_normalized(1:ao_num,1:mo_num),size(coef_normalized,1) * 1_8, &
|
||||||
|
beta, &
|
||||||
|
mo_vgl(:,ielec,iwalk,5),LDC)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
|
||||||
|
if(info_qmckl_dgemm_value .eq. QMCKL_SUCCESS .and. &
|
||||||
|
info_qmckl_dgemm_Gx .eq. QMCKL_SUCCESS .and. &
|
||||||
|
info_qmckl_dgemm_Gy .eq. QMCKL_SUCCESS .and. &
|
||||||
|
info_qmckl_dgemm_Gz .eq. QMCKL_SUCCESS .and. &
|
||||||
|
info_qmckl_dgemm_lap .eq. QMCKL_SUCCESS ) then
|
||||||
|
info = QMCKL_SUCCESS
|
||||||
|
else
|
||||||
|
info = QMCKL_FAILURE
|
||||||
|
end if
|
||||||
|
|
||||||
|
end function qmckl_compute_mo_basis_gaussian_vgl_f
|
||||||
|
#+end_src
|
||||||
|
|
||||||
|
#+CALL: generate_c_header(table=qmckl_mo_basis_gaussian_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_gaussian_vgl"))
|
||||||
|
|
||||||
|
#+RESULTS:
|
||||||
|
#+begin_src c :tangle (eval h_func) :comments org
|
||||||
|
qmckl_exit_code qmckl_compute_mo_basis_gaussian_vgl (
|
||||||
|
const qmckl_context context,
|
||||||
|
const int64_t ao_num,
|
||||||
|
const int64_t mo_num,
|
||||||
|
const int64_t elec_num,
|
||||||
|
const int64_t walk_num,
|
||||||
|
const double* coef_normalized,
|
||||||
|
const double* ao_vgl,
|
||||||
|
double* const mo_vgl );
|
||||||
|
#+end_src
|
||||||
|
|
||||||
|
|
||||||
|
#+CALL: generate_c_interface(table=qmckl_mo_basis_gaussian_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_gaussian_vgl"))
|
||||||
|
|
||||||
|
#+RESULTS:
|
||||||
|
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
||||||
|
integer(c_int32_t) function qmckl_compute_mo_basis_gaussian_vgl &
|
||||||
|
(context, ao_num, mo_num, elec_num, walk_num, coef_normalized, ao_vgl, mo_vgl) &
|
||||||
|
bind(C) result(info)
|
||||||
|
|
||||||
|
use, intrinsic :: iso_c_binding
|
||||||
|
implicit none
|
||||||
|
|
||||||
|
integer (c_int64_t) , intent(in) , value :: context
|
||||||
|
integer (c_int64_t) , intent(in) , value :: ao_num
|
||||||
|
integer (c_int64_t) , intent(in) , value :: mo_num
|
||||||
|
integer (c_int64_t) , intent(in) , value :: elec_num
|
||||||
|
integer (c_int64_t) , intent(in) , value :: walk_num
|
||||||
|
real (c_double ) , intent(in) :: coef_normalized(ao_num,mo_num)
|
||||||
|
real (c_double ) , intent(in) :: ao_vgl(ao_num,elec_num,walk_num,5)
|
||||||
|
real (c_double ) , intent(out) :: mo_vgl(mo_num,elec_num,walk_num,5)
|
||||||
|
|
||||||
|
integer(c_int32_t), external :: qmckl_compute_mo_basis_gaussian_vgl_f
|
||||||
|
info = qmckl_compute_mo_basis_gaussian_vgl_f &
|
||||||
|
(context, ao_num, mo_num, elec_num, walk_num, coef_normalized, ao_vgl, mo_vgl)
|
||||||
|
|
||||||
|
end function qmckl_compute_mo_basis_gaussian_vgl
|
||||||
|
#+end_src
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
*** Test
|
||||||
|
|
||||||
|
#+begin_src python :results output :exports none
|
||||||
|
import numpy as np
|
||||||
|
|
||||||
|
def f(a,x,y):
|
||||||
|
return np.sum( [c * np.exp( -b*(np.linalg.norm(x-y))**2) for b,c in a] )
|
||||||
|
|
||||||
|
def df(a,x,y,n):
|
||||||
|
h0 = 1.e-6
|
||||||
|
if n == 1: h = np.array([h0,0.,0.])
|
||||||
|
elif n == 2: h = np.array([0.,h0,0.])
|
||||||
|
elif n == 3: h = np.array([0.,0.,h0])
|
||||||
|
return ( f(a,x+h,y) - f(a,x-h,y) ) / (2.*h0)
|
||||||
|
|
||||||
|
def d2f(a,x,y,n):
|
||||||
|
h0 = 1.e-6
|
||||||
|
if n == 1: h = np.array([h0,0.,0.])
|
||||||
|
elif n == 2: h = np.array([0.,h0,0.])
|
||||||
|
elif n == 3: h = np.array([0.,0.,h0])
|
||||||
|
return ( f(a,x+h,y) - 2.*f(a,x,y) + f(a,x-h,y) ) / h0**2
|
||||||
|
|
||||||
|
def lf(a,x,y):
|
||||||
|
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_15_w2 = np.array( [ -2.20180344582,-1.9113150239, 2.2193744778600002 ] )
|
||||||
|
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 ] )
|
||||||
|
|
||||||
|
#double prim_vgl[prim_num][5][walk_num][elec_num];
|
||||||
|
x = elec_26_w1 ; y = nucl_1
|
||||||
|
a = [( 8.236000E+03, -1.130000E-04 * 6.1616545431994848e+02 ),
|
||||||
|
( 1.235000E+03, -8.780000E-04 * 1.4847738511079908e+02 ),
|
||||||
|
( 2.808000E+02, -4.540000E-03 * 4.8888635917437597e+01 ),
|
||||||
|
( 7.927000E+01, -1.813300E-02 * 1.8933972232608955e+01 ),
|
||||||
|
( 2.559000E+01, -5.576000E-02 * 8.1089160941724145e+00 ),
|
||||||
|
( 8.997000E+00, -1.268950E-01 * 3.7024003863155635e+00 ),
|
||||||
|
( 3.319000E+00, -1.703520E-01 * 1.7525302846177560e+00 ),
|
||||||
|
( 9.059000E-01, 1.403820E-01 * 6.6179013183966806e-01 ),
|
||||||
|
( 3.643000E-01, 5.986840E-01 * 3.3419848027174592e-01 ),
|
||||||
|
( 1.285000E-01, 3.953890E-01 * 1.5296336817449557e-01 )]
|
||||||
|
|
||||||
|
print ( "[1][0][0][26] : %25.15e"% f(a,x,y))
|
||||||
|
print ( "[1][1][0][26] : %25.15e"% df(a,x,y,1))
|
||||||
|
print ( "[1][2][0][26] : %25.15e"% df(a,x,y,2))
|
||||||
|
print ( "[1][3][0][26] : %25.15e"% df(a,x,y,3))
|
||||||
|
print ( "[1][4][0][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
|
||||||
|
|
||||||
|
#+begin_src c :tangle (eval c_test) :exports none
|
||||||
|
{
|
||||||
|
#define walk_num chbrclf_walk_num
|
||||||
|
#define elec_num chbrclf_elec_num
|
||||||
|
#define shell_num chbrclf_shell_num
|
||||||
|
#define ao_num chbrclf_ao_num
|
||||||
|
|
||||||
|
int64_t elec_up_num = chbrclf_elec_up_num;
|
||||||
|
int64_t elec_dn_num = chbrclf_elec_dn_num;
|
||||||
|
double* elec_coord = &(chbrclf_elec_coord[0][0][0]);
|
||||||
|
const int64_t nucl_num = chbrclf_nucl_num;
|
||||||
|
const double* nucl_charge = chbrclf_charge;
|
||||||
|
const double* nucl_coord = &(chbrclf_nucl_coord[0][0]);
|
||||||
|
|
||||||
|
rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num);
|
||||||
|
assert (rc == QMCKL_SUCCESS);
|
||||||
|
|
||||||
|
rc = qmckl_set_electron_walk_num (context, walk_num);
|
||||||
|
assert (rc == QMCKL_SUCCESS);
|
||||||
|
|
||||||
|
assert(qmckl_electron_provided(context));
|
||||||
|
|
||||||
|
rc = qmckl_set_electron_coord (context, 'N', elec_coord);
|
||||||
|
assert(rc == QMCKL_SUCCESS);
|
||||||
|
|
||||||
|
rc = qmckl_set_nucleus_num (context, nucl_num);
|
||||||
|
assert(rc == QMCKL_SUCCESS);
|
||||||
|
|
||||||
|
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]));
|
||||||
|
assert(rc == QMCKL_SUCCESS);
|
||||||
|
|
||||||
|
rc = qmckl_set_nucleus_charge(context, nucl_charge);
|
||||||
|
assert(rc == QMCKL_SUCCESS);
|
||||||
|
|
||||||
|
assert(qmckl_nucleus_provided(context));
|
||||||
|
|
||||||
|
const int64_t * nucleus_index = &(chbrclf_basis_nucleus_index[0]);
|
||||||
|
const int64_t * nucleus_shell_num = &(chbrclf_basis_nucleus_shell_num[0]);
|
||||||
|
const int32_t * shell_ang_mom = &(chbrclf_basis_shell_ang_mom[0]);
|
||||||
|
const int64_t * shell_prim_num = &(chbrclf_basis_shell_prim_num[0]);
|
||||||
|
const int64_t * shell_prim_index = &(chbrclf_basis_shell_prim_index[0]);
|
||||||
|
const double * shell_factor = &(chbrclf_basis_shell_factor[0]);
|
||||||
|
const double * exponent = &(chbrclf_basis_exponent[0]);
|
||||||
|
const double * coefficient = &(chbrclf_basis_coefficient[0]);
|
||||||
|
const double * prim_factor = &(chbrclf_basis_prim_factor[0]);
|
||||||
|
const double * ao_factor = &(chbrclf_basis_ao_factor[0]);
|
||||||
|
|
||||||
|
const char typ = 'G';
|
||||||
|
|
||||||
|
assert(!qmckl_ao_basis_provided(context));
|
||||||
|
|
||||||
|
rc = qmckl_set_ao_basis_type (context, typ);
|
||||||
|
assert(rc == QMCKL_SUCCESS);
|
||||||
|
assert(!qmckl_ao_basis_provided(context));
|
||||||
|
|
||||||
|
rc = qmckl_set_ao_basis_shell_num (context, chbrclf_shell_num);
|
||||||
|
assert(rc == QMCKL_SUCCESS);
|
||||||
|
assert(!qmckl_ao_basis_provided(context));
|
||||||
|
|
||||||
|
rc = qmckl_set_ao_basis_prim_num (context, chbrclf_prim_num);
|
||||||
|
assert(rc == QMCKL_SUCCESS);
|
||||||
|
assert(!qmckl_ao_basis_provided(context));
|
||||||
|
|
||||||
|
rc = qmckl_set_ao_basis_nucleus_index (context, nucleus_index);
|
||||||
|
assert(rc == QMCKL_SUCCESS);
|
||||||
|
assert(!qmckl_ao_basis_provided(context));
|
||||||
|
|
||||||
|
rc = qmckl_set_ao_basis_nucleus_shell_num (context, nucleus_shell_num);
|
||||||
|
assert(rc == QMCKL_SUCCESS);
|
||||||
|
assert(!qmckl_ao_basis_provided(context));
|
||||||
|
|
||||||
|
rc = qmckl_set_ao_basis_shell_ang_mom (context, shell_ang_mom);
|
||||||
|
assert(rc == QMCKL_SUCCESS);
|
||||||
|
assert(!qmckl_ao_basis_provided(context));
|
||||||
|
|
||||||
|
rc = qmckl_set_ao_basis_shell_factor (context, shell_factor);
|
||||||
|
assert(rc == QMCKL_SUCCESS);
|
||||||
|
assert(!qmckl_ao_basis_provided(context));
|
||||||
|
|
||||||
|
rc = qmckl_set_ao_basis_shell_prim_num (context, shell_prim_num);
|
||||||
|
assert(rc == QMCKL_SUCCESS);
|
||||||
|
assert(!qmckl_ao_basis_provided(context));
|
||||||
|
|
||||||
|
rc = qmckl_set_ao_basis_shell_prim_index (context, shell_prim_index);
|
||||||
|
assert(rc == QMCKL_SUCCESS);
|
||||||
|
assert(!qmckl_ao_basis_provided(context));
|
||||||
|
|
||||||
|
rc = qmckl_set_ao_basis_exponent (context, exponent);
|
||||||
|
assert(rc == QMCKL_SUCCESS);
|
||||||
|
assert(!qmckl_ao_basis_provided(context));
|
||||||
|
|
||||||
|
rc = qmckl_set_ao_basis_coefficient (context, coefficient);
|
||||||
|
assert(rc == QMCKL_SUCCESS);
|
||||||
|
assert(!qmckl_ao_basis_provided(context));
|
||||||
|
|
||||||
|
rc = qmckl_set_ao_basis_prim_factor (context, prim_factor);
|
||||||
|
assert(rc == QMCKL_SUCCESS);
|
||||||
|
|
||||||
|
rc = qmckl_set_ao_basis_ao_num(context, chbrclf_ao_num);
|
||||||
|
assert(rc == QMCKL_SUCCESS);
|
||||||
|
|
||||||
|
rc = qmckl_set_ao_basis_ao_factor (context, ao_factor);
|
||||||
|
assert(rc == QMCKL_SUCCESS);
|
||||||
|
|
||||||
|
assert(qmckl_ao_basis_provided(context));
|
||||||
|
|
||||||
|
|
||||||
|
double ao_vgl[5][walk_num][elec_num][chbrclf_ao_num];
|
||||||
|
|
||||||
|
rc = qmckl_get_ao_vgl(context, &(ao_vgl[0][0][0][0]));
|
||||||
|
assert (rc == QMCKL_SUCCESS);
|
||||||
|
|
||||||
|
/* Set up MO data */
|
||||||
|
rc = qmckl_set_mo_basis_type(context, typ);
|
||||||
|
assert (rc == QMCKL_SUCCESS);
|
||||||
|
|
||||||
|
const int64_t mo_num = chbrclf_mo_num;
|
||||||
|
rc = qmckl_set_mo_basis_mo_num(context, mo_num);
|
||||||
|
assert (rc == QMCKL_SUCCESS);
|
||||||
|
|
||||||
|
const double * mo_coefficient = &(chbrclf_mo_coef[0]);
|
||||||
|
|
||||||
|
rc = qmckl_set_mo_basis_coefficient(context, mo_coefficient);
|
||||||
|
assert (rc == QMCKL_SUCCESS);
|
||||||
|
|
||||||
|
assert(qmckl_mo_basis_provided(context));
|
||||||
|
|
||||||
|
double mo_vgl[5][walk_num][elec_num][chbrclf_mo_num];
|
||||||
|
rc = qmckl_get_mo_basis_vgl(context, &(mo_vgl[0][0][0][0]));
|
||||||
|
assert (rc == QMCKL_SUCCESS);
|
||||||
|
|
||||||
|
// Test overlap of MO
|
||||||
|
//double point_x[100];
|
||||||
|
//double point_y[100];
|
||||||
|
//double point_z[100];
|
||||||
|
//int32_t npoints=100;
|
||||||
|
//// obtain points
|
||||||
|
//double dr = 20./(npoints-1);
|
||||||
|
//double dr3 = dr*dr*dr;
|
||||||
|
//
|
||||||
|
//for (int i=0;i<npoints;++i) {
|
||||||
|
// point_x[i] = -10. + dr*i;
|
||||||
|
// point_y[i] = -10. + dr*i;
|
||||||
|
// point_z[i] = -10. + dr*i;
|
||||||
|
//}
|
||||||
|
//
|
||||||
|
//double ovlmo1 = 0.0;
|
||||||
|
//// Calculate overlap
|
||||||
|
//for (int i=0;i<npoints;++i) {
|
||||||
|
// printf(".");
|
||||||
|
// fflush(stdout);
|
||||||
|
// for (int j=0;j<npoints;++j) {
|
||||||
|
// for (int k=0;k<npoints;++k) {
|
||||||
|
// // Set point
|
||||||
|
// elec_coord[0] = point_x[i];
|
||||||
|
// elec_coord[1] = point_y[j];
|
||||||
|
// elec_coord[2] = point_z[k];
|
||||||
|
// rc = qmckl_set_electron_coord (context, 'N', elec_coord);
|
||||||
|
// assert(rc == QMCKL_SUCCESS);
|
||||||
|
//
|
||||||
|
// // Calculate value of MO (1st electron)
|
||||||
|
// double mo_vgl[5][walk_num][elec_num][chbrclf_mo_num];
|
||||||
|
// rc = qmckl_get_mo_basis_vgl(context, &(mo_vgl[0][0][0][0]));
|
||||||
|
// assert (rc == QMCKL_SUCCESS);
|
||||||
|
// ovlmo1 += mo_vgl[0][0][0][0]*mo_vgl[0][0][0][0]*dr3;
|
||||||
|
// }
|
||||||
|
// }
|
||||||
|
//}
|
||||||
|
//printf("OVL MO1 = %10.15f\n",ovlmo1);
|
||||||
|
|
||||||
|
|
||||||
|
printf("\n");
|
||||||
|
printf(" mo_vgl mo_vgl[0][0][26][219] %25.15e\n", mo_vgl[0][0][2][3]);
|
||||||
|
printf(" mo_vgl mo_vgl[1][0][26][219] %25.15e\n", mo_vgl[1][0][2][3]);
|
||||||
|
printf(" mo_vgl mo_vgl[0][0][26][220] %25.15e\n", mo_vgl[0][0][2][3]);
|
||||||
|
printf(" mo_vgl mo_vgl[1][0][26][220] %25.15e\n", mo_vgl[1][0][2][3]);
|
||||||
|
printf(" mo_vgl mo_vgl[0][0][26][221] %25.15e\n", mo_vgl[0][0][2][3]);
|
||||||
|
printf(" mo_vgl mo_vgl[1][0][26][221] %25.15e\n", mo_vgl[1][0][2][3]);
|
||||||
|
printf(" mo_vgl mo_vgl[0][0][26][222] %25.15e\n", mo_vgl[0][0][2][3]);
|
||||||
|
printf(" mo_vgl mo_vgl[1][0][26][222] %25.15e\n", mo_vgl[1][0][2][3]);
|
||||||
|
printf(" mo_vgl mo_vgl[0][0][26][223] %25.15e\n", mo_vgl[0][0][2][3]);
|
||||||
|
printf(" mo_vgl mo_vgl[1][0][26][223] %25.15e\n", mo_vgl[1][0][2][3]);
|
||||||
|
printf(" mo_vgl mo_vgl[0][0][26][224] %25.15e\n", mo_vgl[0][0][2][3]);
|
||||||
|
printf(" mo_vgl mo_vgl[1][0][26][224] %25.15e\n", mo_vgl[1][0][2][3]);
|
||||||
|
printf("\n");
|
||||||
|
}
|
||||||
|
|
||||||
|
#+end_src
|
||||||
|
|
||||||
|
|
||||||
|
* End of files :noexport:
|
||||||
|
|
||||||
|
#+begin_src c :tangle (eval h_private_type)
|
||||||
|
#endif
|
||||||
|
#+end_src
|
||||||
|
|
||||||
|
*** Test
|
||||||
|
#+begin_src c :tangle (eval c_test)
|
||||||
|
rc = qmckl_context_destroy(context);
|
||||||
|
assert (rc == QMCKL_SUCCESS);
|
||||||
|
|
||||||
|
return 0;
|
||||||
|
}
|
||||||
|
#+end_src
|
||||||
|
|
||||||
|
*** Compute file names
|
||||||
|
#+begin_src emacs-lisp
|
||||||
|
; The following is required to compute the file names
|
||||||
|
|
||||||
|
(setq pwd (file-name-directory buffer-file-name))
|
||||||
|
(setq name (file-name-nondirectory (substring buffer-file-name 0 -4)))
|
||||||
|
(setq f (concat pwd name "_f.f90"))
|
||||||
|
(setq fh (concat pwd name "_fh.f90"))
|
||||||
|
(setq c (concat pwd name ".c"))
|
||||||
|
(setq h (concat name ".h"))
|
||||||
|
(setq h_private (concat name "_private.h"))
|
||||||
|
(setq c_test (concat pwd "test_" name ".c"))
|
||||||
|
(setq f_test (concat pwd "test_" name "_f.f90"))
|
||||||
|
|
||||||
|
; Minted
|
||||||
|
(require 'ox-latex)
|
||||||
|
(setq org-latex-listings 'minted)
|
||||||
|
(add-to-list 'org-latex-packages-alist '("" "listings"))
|
||||||
|
(add-to-list 'org-latex-packages-alist '("" "color"))
|
||||||
|
|
||||||
|
#+end_src
|
||||||
|
|
||||||
|
|
||||||
|
# -*- mode: org -*-
|
||||||
|
# vim: syntax=c
|
||||||
|
|
||||||
|
|
File diff suppressed because one or more lines are too long
59032
org/qmckl_tests.org
59032
org/qmckl_tests.org
File diff suppressed because it is too large
Load Diff
@ -6,8 +6,10 @@ qmckl_numprec.org
|
|||||||
qmckl_nucleus.org
|
qmckl_nucleus.org
|
||||||
qmckl_electron.org
|
qmckl_electron.org
|
||||||
qmckl_ao.org
|
qmckl_ao.org
|
||||||
|
qmckl_mo.org
|
||||||
qmckl_jastrow.org
|
qmckl_jastrow.org
|
||||||
qmckl_sherman_morrison_woodbury.org
|
qmckl_sherman_morrison_woodbury.org
|
||||||
qmckl_distance.org
|
qmckl_distance.org
|
||||||
qmckl_utils.org
|
qmckl_utils.org
|
||||||
|
qmckl_blas.org
|
||||||
qmckl_tests.org
|
qmckl_tests.org
|
||||||
|
@ -11,11 +11,10 @@
|
|||||||
** Function to get the value of a property.
|
** Function to get the value of a property.
|
||||||
#+NAME: get_value
|
#+NAME: get_value
|
||||||
#+begin_src elisp :var key="Type"
|
#+begin_src elisp :var key="Type"
|
||||||
(setq x (org-property-values key))
|
(org-with-point-at org-babel-current-src-block-location
|
||||||
(pop x)
|
(org-entry-get nil key t))
|
||||||
#+end_src
|
#+end_src
|
||||||
|
|
||||||
#+RESULTS: get_value
|
|
||||||
|
|
||||||
** Table of function arguments
|
** Table of function arguments
|
||||||
|
|
||||||
@ -40,6 +39,7 @@
|
|||||||
f_of_c_d = { '' : ''
|
f_of_c_d = { '' : ''
|
||||||
, 'qmckl_context' : 'integer (c_int64_t)'
|
, 'qmckl_context' : 'integer (c_int64_t)'
|
||||||
, 'qmckl_exit_code' : 'integer (c_int32_t)'
|
, 'qmckl_exit_code' : 'integer (c_int32_t)'
|
||||||
|
, 'bool' : 'logical*8'
|
||||||
, 'int32_t' : 'integer (c_int32_t)'
|
, 'int32_t' : 'integer (c_int32_t)'
|
||||||
, 'int64_t' : 'integer (c_int64_t)'
|
, 'int64_t' : 'integer (c_int64_t)'
|
||||||
, 'uint32_t' : 'integer (c_int32_t)'
|
, 'uint32_t' : 'integer (c_int32_t)'
|
||||||
@ -62,6 +62,7 @@ ctypeid_d = { '' : ''
|
|||||||
, 'real' : 'real(c_float)'
|
, 'real' : 'real(c_float)'
|
||||||
, 'real*8' : 'real(c_double)'
|
, 'real*8' : 'real(c_double)'
|
||||||
, 'character' : 'character(c_char)'
|
, 'character' : 'character(c_char)'
|
||||||
|
, 'character' : 'character(c_char)'
|
||||||
}
|
}
|
||||||
#+END_SRC
|
#+END_SRC
|
||||||
|
|
||||||
@ -188,23 +189,6 @@ results='\n'.join(results)
|
|||||||
return results
|
return results
|
||||||
#+END_SRC
|
#+END_SRC
|
||||||
|
|
||||||
#+RESULTS: generate_c_interface
|
|
||||||
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
|
||||||
integer(c_int32_t) function [] &
|
|
||||||
() &
|
|
||||||
bind(C) result(info)
|
|
||||||
|
|
||||||
use, intrinsic :: iso_c_binding
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
|
|
||||||
integer(c_int32_t), external :: []_f
|
|
||||||
info = []_f &
|
|
||||||
()
|
|
||||||
|
|
||||||
end function []
|
|
||||||
#+end_src
|
|
||||||
|
|
||||||
*** Generates a Fortran interface to the C function
|
*** Generates a Fortran interface to the C function
|
||||||
|
|
||||||
#+NAME: generate_f_interface
|
#+NAME: generate_f_interface
|
||||||
@ -256,8 +240,5 @@ results='\n'.join(results)
|
|||||||
return results
|
return results
|
||||||
#+END_SRC
|
#+END_SRC
|
||||||
|
|
||||||
#+RESULTS: generate_f_interface
|
|
||||||
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
|
|
||||||
#+end_src
|
|
||||||
|
|
||||||
|
|
||||||
|
Loading…
Reference in New Issue
Block a user