1
0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-09-16 17:55:29 +02:00

Accelerated AOs

This commit is contained in:
Anthony Scemama 2022-02-14 19:11:37 +01:00
parent 2cb8c7b63f
commit 685b7201fc
2 changed files with 580 additions and 211 deletions

View File

@ -3387,6 +3387,7 @@ for (j=0 ; j<point_num ; ++j) {
| ~nucl_num~ | ~int64_t~ | in | Number of nuclei |
| ~nucleus_shell_num~ | ~int64_t[nucl_num]~ | in | Number of shells for each nucleus |
| ~nucleus_index~ | ~int64_t[nucl_num]~ | in | Index of the 1st shell of each nucleus |
| ~nucleus_range~ | ~double[nucl_num]~ | in | Range of the nucleus |
| ~shell_prim_index~ | ~int64_t[shell_num]~ | in | Index of the 1st primitive of each shell |
| ~shell_prim_num~ | ~int64_t[shell_num]~ | in | Number of primitives per shell |
| ~coord~ | ~double[3][point_num]~ | in | Coordinates |
@ -3407,6 +3408,7 @@ for (j=0 ; j<point_num ; ++j) {
const int64_t nucl_num,
const int64_t* nucleus_shell_num,
const int64_t* nucleus_index,
const double* nucleus_range,
const int64_t* shell_prim_index,
const int64_t* shell_prim_num,
const double* coord,
@ -3419,9 +3421,9 @@ for (j=0 ; j<point_num ; ++j) {
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_ao_basis_shell_gaussian_vgl_f( &
context, prim_num, shell_num, point_num, nucl_num, &
nucleus_shell_num, nucleus_index, shell_prim_index, &
shell_prim_num, coord, nucl_coord, expo, &
coef_normalized, shell_vgl) &
nucleus_shell_num, nucleus_index, nucleus_range, &
shell_prim_index, shell_prim_num, coord, nucl_coord, &
expo, coef_normalized, shell_vgl) &
result(info)
use qmckl
implicit none
@ -3432,6 +3434,7 @@ integer function qmckl_compute_ao_basis_shell_gaussian_vgl_f( &
integer*8 , intent(in) :: point_num
integer*8 , intent(in) :: nucleus_shell_num(nucl_num)
integer*8 , intent(in) :: nucleus_index(nucl_num)
double precision , intent(in) :: nucleus_range(nucl_num)
integer*8 , intent(in) :: shell_prim_index(shell_num)
integer*8 , intent(in) :: shell_prim_num(shell_num)
double precision , intent(in) :: coord(point_num,3)
@ -3451,13 +3454,9 @@ integer function qmckl_compute_ao_basis_shell_gaussian_vgl_f( &
! TODO : Use numerical precision here
cutoff = -dlog(1.d-12)
do inucl=1,nucl_num
do ipoint = 1, point_num
! C is zero-based, so shift bounds by one
ishell_start = nucleus_index(inucl) + 1
ishell_end = nucleus_index(inucl) + nucleus_shell_num(inucl)
do ipoint = 1, point_num
do inucl=1,nucl_num
x = coord(ipoint,1) - nucl_coord(inucl,1)
y = coord(ipoint,2) - nucl_coord(inucl,2)
@ -3465,6 +3464,14 @@ integer function qmckl_compute_ao_basis_shell_gaussian_vgl_f( &
r2 = x*x + y*y + z*z
if (r2 > cutoff*nucleus_range(inucl)) then
cycle
end if
! C is zero-based, so shift bounds by one
ishell_start = nucleus_index(inucl) + 1
ishell_end = nucleus_index(inucl) + nucleus_shell_num(inucl)
do ishell=ishell_start, ishell_end
shell_vgl(ishell, 1, ipoint) = 0.d0
@ -3523,6 +3530,7 @@ end function qmckl_compute_ao_basis_shell_gaussian_vgl_f
nucl_num, &
nucleus_shell_num, &
nucleus_index, &
nucleus_range, &
shell_prim_index, &
shell_prim_num, &
coord, &
@ -3542,6 +3550,7 @@ end function qmckl_compute_ao_basis_shell_gaussian_vgl_f
integer (c_int64_t) , intent(in) , value :: nucl_num
integer (c_int64_t) , intent(in) :: nucleus_shell_num(nucl_num)
integer (c_int64_t) , intent(in) :: nucleus_index(nucl_num)
real (c_double ) , intent(in) :: nucleus_range(nucl_num)
integer (c_int64_t) , intent(in) :: shell_prim_index(shell_num)
integer (c_int64_t) , intent(in) :: shell_prim_num(shell_num)
real (c_double ) , intent(in) :: coord(point_num,3)
@ -3559,6 +3568,7 @@ end function qmckl_compute_ao_basis_shell_gaussian_vgl_f
nucl_num, &
nucleus_shell_num, &
nucleus_index, &
nucleus_range, &
shell_prim_index, &
shell_prim_num, &
coord, &
@ -3625,6 +3635,7 @@ qmckl_exit_code qmckl_provide_ao_basis_shell_vgl(qmckl_context context)
ctx->nucleus.num,
ctx->ao_basis.nucleus_shell_num,
ctx->ao_basis.nucleus_index,
ctx->ao_basis.nucleus_range,
ctx->ao_basis.shell_prim_index,
ctx->ao_basis.shell_prim_num,
ctx->point.coord.data,
@ -4263,7 +4274,7 @@ integer function qmckl_ao_polynomial_transp_vgl_f (context, &
integer :: a,b,c,d
real*8 :: Y(3)
integer :: lmax_array(3)
real*8 :: pows(-2:lmax,3)
real*8 :: pows(-2:21,3) ! lmax < 22
double precision :: xy, yz, xz
double precision :: da, db, dc, dd
@ -4290,17 +4301,12 @@ integer function qmckl_ao_polynomial_transp_vgl_f (context, &
endif
do i=1,3
Y(i) = X(i) - R(i)
end do
lmax_array(1:3) = lmax
if (lmax == 0) then
VGL(1,1) = 1.d0
VGL(1,2:5) = 0.d0
l(1:3,1) = 0
n=1
else if (lmax > 0) then
if (lmax > 0) then
do i=1,3
Y(i) = X(i) - R(i)
end do
pows(-2:0,1:3) = 1.d0
do i=1,lmax
pows(i,1) = pows(i-1,1) * Y(1)
@ -4327,6 +4333,12 @@ integer function qmckl_ao_polynomial_transp_vgl_f (context, &
VGL(4,4) = 1.d0
n=4
else
VGL(1,1) = 1.d0
VGL(1,2:5) = 0.d0
l(1:3,1) = 0
n=1
return
endif
! l>=2
@ -4535,7 +4547,8 @@ end function test_qmckl_ao_polynomial_vgl
:FRetType: qmckl_exit_code
:END:
#+NAME: qmckl_ao_vgl_args
** Unoptimized version
#+NAME: qmckl_ao_vgl_args_doc
| Variable | Type | In/Out | Description |
|-----------------------+-----------------------------------+--------+----------------------------------------------|
| ~context~ | ~qmckl_context~ | in | Global state |
@ -4552,9 +4565,8 @@ end function test_qmckl_ao_polynomial_vgl
| ~shell_ang_mom~ | ~int32_t[shell_num]~ | in | Angular momentum of each shell |
| ~ao_factor~ | ~double[ao_num]~ | in | Normalization factor of the AOs |
| ~shell_vgl~ | ~double[point_num][5][shell_num]~ | in | Value, gradients and Laplacian of the shells |
| ~ao_vgl~ | ~double[point_num][5][ao_num]~ | out | Value, gradients and Laplacian of the AOs |
| ~ao_vgl~ | ~double[point_num][5][ao_num]~ | out | Value, gradients and Laplacian of the AOs |
** Unoptimized version
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_ao_vgl_doc_f(context, &
ao_num, shell_num, point_num, nucl_num, &
@ -4591,15 +4603,25 @@ integer function qmckl_compute_ao_vgl_doc_f(context, &
integer, external :: qmckl_ao_polynomial_vgl_f
double precision, allocatable :: poly_vgl(:,:)
integer , allocatable :: powers(:,:)
integer , allocatable :: powers(:,:), ao_index(:)
allocate(poly_vgl(5,ao_num), powers(3,ao_num))
allocate(poly_vgl(5,ao_num), powers(3,ao_num), ao_index(ao_num))
! Pre-computed data
do l=0,20
lstart(l) = l*(l+1)*(l+2)/6 +1
end do
k=1
do inucl=1,nucl_num
ishell_start = nucleus_index(inucl) + 1
ishell_end = nucleus_index(inucl) + nucleus_shell_num(inucl)
do ishell = ishell_start, ishell_end
l = shell_ang_mom(ishell)
ao_index(ishell) = k
k = k + lstart(l+1) - lstart(l)
end do
end do
info = QMCKL_SUCCESS
! Don't compute polynomials when the radial part is zero.
@ -4609,7 +4631,6 @@ integer function qmckl_compute_ao_vgl_doc_f(context, &
e_coord(1) = coord(ipoint,1)
e_coord(2) = coord(ipoint,2)
e_coord(3) = coord(ipoint,3)
k=1
do inucl=1,nucl_num
n_coord(1) = nucl_coord(inucl,1)
n_coord(2) = nucl_coord(inucl,2)
@ -4620,7 +4641,7 @@ integer function qmckl_compute_ao_vgl_doc_f(context, &
y = e_coord(2) - n_coord(2)
z = e_coord(3) - n_coord(3)
r2 = x*x + z*z + z*z
r2 = x*x + y*y + z*z
if (r2 > cutoff*nucleus_range(inucl)) then
cycle
@ -4635,6 +4656,7 @@ integer function qmckl_compute_ao_vgl_doc_f(context, &
ishell_start = nucleus_index(inucl) + 1
ishell_end = nucleus_index(inucl) + nucleus_shell_num(inucl)
do ishell = ishell_start, ishell_end
k = ao_index(ishell)
l = shell_ang_mom(ishell)
do il = lstart(l), lstart(l+1)-1
! Value
@ -4680,18 +4702,43 @@ end function qmckl_compute_ao_vgl_doc_f
#+end_src
** HPC version
#+NAME: qmckl_ao_vgl_args_hpc
| Variable | Type | In/Out | Description |
|-----------------------+--------------------------------+--------+----------------------------------------------|
| ~context~ | ~qmckl_context~ | in | Global state |
| ~ao_num~ | ~int64_t~ | in | Number of AOs |
| ~shell_num~ | ~int64_t~ | in | Number of shells |
| ~prim_num~ | ~int64_t~ | in | Number of primitives |
| ~point_num~ | ~int64_t~ | in | Number of points |
| ~nucl_num~ | ~int64_t~ | in | Number of nuclei |
| ~coord~ | ~double[3][point_num]~ | in | Coordinates |
| ~nucl_coord~ | ~double[3][nucl_num]~ | in | Nuclear coordinates |
| ~nucleus_index~ | ~int64_t[nucl_num]~ | in | Index of the 1st shell of each nucleus |
| ~nucleus_shell_num~ | ~int64_t[nucl_num]~ | in | Number of shells per nucleus |
| ~nucleus_range~ | ~double[nucl_num]~ | in | Range beyond which all is zero |
| ~nucleus_max_ang_mom~ | ~int32_t[nucl_num]~ | in | Maximum angular momentum per nucleus |
| ~shell_ang_mom~ | ~int32_t[shell_num]~ | in | Angular momentum of each shell |
| ~shell_prim_index~ | ~int64_t[shell_num]~ | in | Index of the 1st primitive of each shell |
| ~shell_prim_num~ | ~int64_t[shell_num]~ | in | Number of primitives per shell |
| ~ao_factor~ | ~double[ao_num]~ | in | Normalization factor of the AOs |
| ~ao_expo~ | ~double[prim_num]~ | in | Value, gradients and Laplacian of the shells |
| ~coef_normalized~ | ~double[prim_num]~ | in | Value, gradients and Laplacian of the shells |
| ~ao_vgl~ | ~double[point_num][5][ao_num]~ | out | Value, gradients and Laplacian of the AOs |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_ao_vgl_hpc_f(context, &
ao_num, shell_num, point_num, nucl_num, &
integer function qmckl_compute_ao_vgl_hpc_f(context, &
ao_num, shell_num, prim_num, point_num, nucl_num, &
coord, nucl_coord, nucleus_index, nucleus_shell_num, &
nucleus_range, nucleus_max_ang_mom, shell_ang_mom, &
ao_factor, shell_vgl, ao_vgl) &
nucleus_range, nucleus_max_ang_mom, shell_ang_mom, &
shell_prim_index, shell_prim_num, ao_factor, expo, &
coef_normalized, ao_vgl) &
result(info)
use qmckl
implicit none
integer(qmckl_context), intent(in) :: context
integer*8 , intent(in) :: ao_num
integer*8 , intent(in) :: shell_num
integer*8 , intent(in) :: prim_num
integer*8 , intent(in) :: point_num
integer*8 , intent(in) :: nucl_num
double precision , intent(in) :: coord(point_num,3)
@ -4701,8 +4748,11 @@ integer function qmckl_compute_ao_vgl_hpc_f(context, &
double precision , intent(in) :: nucleus_range(nucl_num)
integer , intent(in) :: nucleus_max_ang_mom(nucl_num)
integer , intent(in) :: shell_ang_mom(shell_num)
integer*8 , intent(in) :: shell_prim_index(shell_num)
integer*8 , intent(in) :: shell_prim_num(shell_num)
double precision , intent(in) :: ao_factor(ao_num)
double precision , intent(in) :: shell_vgl(shell_num,5,point_num)
double precision , intent(in) :: expo(prim_num)
double precision , intent(in) :: coef_normalized(prim_num)
double precision , intent(out) :: ao_vgl(ao_num,5,point_num)
double precision :: e_coord(3), n_coord(3)
@ -4711,20 +4761,36 @@ integer function qmckl_compute_ao_vgl_hpc_f(context, &
integer*8 :: ipoint, inucl, ishell
integer*8 :: ishell_start, ishell_end
integer :: lstart(0:20)
double precision :: x, y, z, r2
double precision :: cutoff
double precision :: x, y, z, r2, s1, s2, s3, s4, s5, s6
double precision :: cutoff, v, two_a
integer*8 :: iprim_start , iprim_end, iprim
integer, external :: qmckl_ao_polynomial_transp_vgl_f
double precision, allocatable :: poly_vgl(:,:)
integer , allocatable :: powers(:,:)
integer , allocatable :: powers(:,:), ao_index(:)
allocate(poly_vgl(ao_num,5), powers(3,ao_num))
integer :: nidx, idx, n
double precision, allocatable :: ar2(:), expo_(:), c_(:)
allocate(poly_vgl(ao_num,5), powers(3,ao_num), ao_index(ao_num))
allocate(c_(prim_num), expo_(prim_num), ar2(prim_num))
! Pre-computed data
do l=0,20
lstart(l) = l*(l+1)*(l+2)/6 +1
end do
k=1
do inucl=1,nucl_num
ishell_start = nucleus_index(inucl) + 1
ishell_end = nucleus_index(inucl) + nucleus_shell_num(inucl)
do ishell = ishell_start, ishell_end
l = shell_ang_mom(ishell)
ao_index(ishell) = k
k = k + lstart(l+1) - lstart(l)
end do
end do
info = QMCKL_SUCCESS
! Don't compute polynomials when the radial part is zero.
@ -4734,7 +4800,6 @@ integer function qmckl_compute_ao_vgl_hpc_f(context, &
e_coord(1) = coord(ipoint,1)
e_coord(2) = coord(ipoint,2)
e_coord(3) = coord(ipoint,3)
k=1
do inucl=1,nucl_num
n_coord(1) = nucl_coord(inucl,1)
n_coord(2) = nucl_coord(inucl,2)
@ -4745,7 +4810,7 @@ integer function qmckl_compute_ao_vgl_hpc_f(context, &
y = e_coord(2) - n_coord(2)
z = e_coord(3) - n_coord(3)
r2 = x*x + z*z + z*z
r2 = x*x + y*y + z*z
if (r2 > cutoff*nucleus_range(inucl)) then
cycle
@ -4759,53 +4824,65 @@ integer function qmckl_compute_ao_vgl_hpc_f(context, &
! Loop over shells
ishell_start = nucleus_index(inucl) + 1
ishell_end = nucleus_index(inucl) + nucleus_shell_num(inucl)
do ishell = ishell_start, ishell_end
iprim_start = shell_prim_index(ishell) + 1
iprim_end = shell_prim_index(ishell) + shell_prim_num(ishell)
! /!\ Gaussian fuctions
nidx = 0
do iprim = iprim_start, iprim_end
v = expo(iprim)*r2
if (v > cutoff) then
cycle
end if
nidx = nidx+1
ar2(nidx) = v
c_(nidx) = coef_normalized(iprim)
expo_(nidx) = expo(iprim)
enddo
s1 = 0.d0
s5 = 0.d0
s6 = 0.d0
do idx = 1, nidx
v = c_(idx) * dexp(-ar2(idx))
s1 = s1 + v
s6 = s6 - expo_(idx) * v
s5 = s5 + ar2(idx)
end do
s6 = s6 +s6
s5 = 2.d0*s5 + s6*3.d0
s2 = s6 * x
s3 = s6 * y
s4 = s6 * z
l = shell_ang_mom(ishell)
if (shell_vgl(ishell,1,ipoint) /= 0.d0) then
do il = lstart(l), lstart(l+1)-1
! Value
ao_vgl(k,1,ipoint) = &
poly_vgl(il,1) * shell_vgl(ishell,1,ipoint) * ao_factor(k)
! Grad_x
ao_vgl(k,2,ipoint) = ( &
poly_vgl(il,2) * shell_vgl(ishell,1,ipoint) + &
poly_vgl(il,1) * shell_vgl(ishell,2,ipoint) &
) * ao_factor(k)
! Grad_y
ao_vgl(k,3,ipoint) = ( &
poly_vgl(il,3) * shell_vgl(ishell,1,ipoint) + &
poly_vgl(il,1) * shell_vgl(ishell,3,ipoint) &
) * ao_factor(k)
! Grad_z
ao_vgl(k,4,ipoint) = ( &
poly_vgl(il,4) * shell_vgl(ishell,1,ipoint) + &
poly_vgl(il,1) * shell_vgl(ishell,4,ipoint) &
) * ao_factor(k)
! Lapl_z
ao_vgl(k,5,ipoint) = ( &
poly_vgl(il,5) * shell_vgl(ishell,1,ipoint) + &
poly_vgl(il,1) * shell_vgl(ishell,5,ipoint) + &
2.d0 * ( &
poly_vgl(il,2) * shell_vgl(ishell,2,ipoint) + &
poly_vgl(il,3) * shell_vgl(ishell,3,ipoint) + &
poly_vgl(il,4) * shell_vgl(ishell,4,ipoint) ) &
) * ao_factor(k)
k = k+1
k = ao_index(ishell)
n = lstart(l+1)-lstart(l)
if (nidx > 0) then
idx = lstart(l)
do il = 0,n-1
ao_vgl(k+il,1,ipoint) = poly_vgl(idx+il,1) * s1 * ao_factor(k+il)
ao_vgl(k+il,2,ipoint) = (poly_vgl(idx+il,2) * s1 + poly_vgl(idx+il,1) * s2) * ao_factor(k+il)
ao_vgl(k+il,3,ipoint) = (poly_vgl(idx+il,3) * s1 + poly_vgl(idx+il,1) * s3) * ao_factor(k+il)
ao_vgl(k+il,4,ipoint) = (poly_vgl(idx+il,4) * s1 + poly_vgl(idx+il,1) * s4) * ao_factor(k+il)
ao_vgl(k+il,5,ipoint) = (poly_vgl(idx+il,5) * s1 + &
poly_vgl(idx+il,1) * s5 + 2.d0*( &
poly_vgl(idx+il,2) * s2 + &
poly_vgl(idx+il,3) * s3 + &
poly_vgl(idx+il,4) * s4 )) * ao_factor(k+il)
end do
else
do il = lstart(l), lstart(l+1)-1
ao_vgl(k,1,ipoint) = 0.d0
ao_vgl(k,2,ipoint) = 0.d0
ao_vgl(k,3,ipoint) = 0.d0
ao_vgl(k,4,ipoint) = 0.d0
ao_vgl(k,5,ipoint) = 0.d0
k = k+1
do il = 0, n-1
ao_vgl(k+il,1,ipoint) = 0.d0
ao_vgl(k+il,2,ipoint) = 0.d0
ao_vgl(k+il,3,ipoint) = 0.d0
ao_vgl(k+il,4,ipoint) = 0.d0
ao_vgl(k+il,5,ipoint) = 0.d0
end do
end if
endif
end do
end do
end do
@ -4815,12 +4892,12 @@ end function qmckl_compute_ao_vgl_hpc_f
#+end_src
** Interfaces
# #+CALL: generate_c_header(table=qmckl_ao_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_ao_vgl"))
# #+CALL: generate_c_header(table=qmckl_ao_vgl_args_doc,rettyp=get_value("CRetType"),fname="qmckl_compute_ao_vgl"))
# (Commented because the header needs to go into h_private_func
#+RESULTS:
#+begin_src c :tangle (eval h_private_func) :comments org
qmckl_exit_code qmckl_compute_ao_vgl (
qmckl_exit_code qmckl_compute_ao_vgl_doc (
const qmckl_context context,
const int64_t ao_num,
const int64_t shell_num,
@ -4837,12 +4914,34 @@ end function qmckl_compute_ao_vgl_hpc_f
const double* shell_vgl,
double* const ao_vgl );
#+end_src
#+begin_src c :tangle (eval h_private_func) :comments org
qmckl_exit_code qmckl_compute_ao_vgl_hpc (
const qmckl_context context,
const int64_t ao_num,
const int64_t shell_num,
const int64_t prim_num,
const int64_t point_num,
const int64_t nucl_num,
const double* coord,
const double* nucl_coord,
const int64_t* nucleus_index,
const int64_t* nucleus_shell_num,
const double* nucleus_range,
const int32_t* nucleus_max_ang_mom,
const int32_t* shell_ang_mom,
const int64_t* shell_prim_index,
const int64_t* shell_prim_num,
const double* ao_factor,
const double* expo,
const double* coef_normalized,
double* const ao_vgl );
#+end_src
#+CALL: generate_c_interface(table=qmckl_ao_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_ao_vgl"))
#+CALL: generate_c_interface(table=qmckl_ao_vgl_args_doc,rettyp=get_value("CRetType"),fname="qmckl_compute_ao_vgl_doc"))
#+RESULTS:
#+begin_src f90 :tangle (eval f) :comments org :exports none
integer(c_int32_t) function qmckl_compute_ao_vgl &
integer(c_int32_t) function qmckl_compute_ao_vgl_doc &
(context, &
ao_num, &
shell_num, &
@ -4879,13 +4978,8 @@ end function qmckl_compute_ao_vgl_hpc_f
real (c_double ) , intent(in) :: shell_vgl(shell_num,5,point_num)
real (c_double ) , intent(out) :: ao_vgl(ao_num,5,point_num)
#ifdef HAVE_HPC
integer(c_int32_t), external :: qmckl_compute_ao_vgl_hpc_f
info = qmckl_compute_ao_vgl_hpc_f &
#else
integer(c_int32_t), external :: qmckl_compute_ao_vgl_doc_f
info = qmckl_compute_ao_vgl_doc_f &
#endif
(context, &
ao_num, &
shell_num, &
@ -4902,7 +4996,81 @@ end function qmckl_compute_ao_vgl_hpc_f
shell_vgl, &
ao_vgl)
end function qmckl_compute_ao_vgl
end function qmckl_compute_ao_vgl_doc
#+end_src
#+CALL: generate_c_interface(table=qmckl_ao_vgl_args_hpc,rettyp=get_value("CRetType"),fname="qmckl_compute_ao_vgl_hpc"))
#+RESULTS:
#+begin_src f90 :tangle (eval f) :comments org :exports none
integer(c_int32_t) function qmckl_compute_ao_vgl_hpc &
(context, &
ao_num, &
shell_num, &
prim_num, &
point_num, &
nucl_num, &
coord, &
nucl_coord, &
nucleus_index, &
nucleus_shell_num, &
nucleus_range, &
nucleus_max_ang_mom, &
shell_ang_mom, &
shell_prim_index, &
shell_prim_num, &
ao_factor, &
ao_expo, &
coef_normalized, &
ao_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 :: shell_num
integer (c_int64_t) , intent(in) , value :: prim_num
integer (c_int64_t) , intent(in) , value :: point_num
integer (c_int64_t) , intent(in) , value :: nucl_num
real (c_double ) , intent(in) :: coord(point_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_shell_num(nucl_num)
real (c_double ) , intent(in) :: nucleus_range(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_int64_t) , intent(in) :: shell_prim_index(shell_num)
integer (c_int64_t) , intent(in) :: shell_prim_num(shell_num)
real (c_double ) , intent(in) :: ao_factor(ao_num)
real (c_double ) , intent(in) :: ao_expo(prim_num)
real (c_double ) , intent(in) :: coef_normalized(prim_num)
real (c_double ) , intent(out) :: ao_vgl(ao_num,5,point_num)
integer(c_int32_t), external :: qmckl_compute_ao_vgl_hpc_f
info = qmckl_compute_ao_vgl_hpc_f &
(context, &
ao_num, &
shell_num, &
prim_num, &
point_num, &
nucl_num, &
coord, &
nucl_coord, &
nucleus_index, &
nucleus_shell_num, &
nucleus_range, &
nucleus_max_ang_mom, &
shell_ang_mom, &
shell_prim_index, &
shell_prim_num, &
ao_factor, &
ao_expo, &
coef_normalized, &
ao_vgl)
end function qmckl_compute_ao_vgl_hpc
#+end_src
*** Provide :noexport:
@ -4938,10 +5106,12 @@ qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context)
qmckl_exit_code rc;
/* Provide required data */
#ifndef HAVE_HPC
rc = qmckl_provide_ao_basis_shell_vgl(context);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context, rc, "qmckl_provide_ao_basis_shell_vgl", NULL);
}
#endif
/* Allocate array */
if (ctx->ao_basis.ao_vgl == NULL) {
@ -4958,22 +5128,43 @@ qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context)
}
ctx->ao_basis.ao_vgl = ao_vgl;
}
rc = qmckl_compute_ao_vgl(context,
ctx->ao_basis.ao_num,
ctx->ao_basis.shell_num,
ctx->point.num,
ctx->nucleus.num,
ctx->point.coord.data,
ctx->nucleus.coord.data,
ctx->ao_basis.nucleus_index,
ctx->ao_basis.nucleus_shell_num,
ctx->ao_basis.nucleus_range,
ctx->ao_basis.nucleus_max_ang_mom,
ctx->ao_basis.shell_ang_mom,
ctx->ao_basis.ao_factor,
ctx->ao_basis.shell_vgl,
ctx->ao_basis.ao_vgl);
#ifdef HAVE_HPC
rc = qmckl_compute_ao_vgl_hpc(context,
ctx->ao_basis.ao_num,
ctx->ao_basis.shell_num,
ctx->ao_basis.prim_num,
ctx->point.num,
ctx->nucleus.num,
ctx->point.coord.data,
ctx->nucleus.coord.data,
ctx->ao_basis.nucleus_index,
ctx->ao_basis.nucleus_shell_num,
ctx->ao_basis.nucleus_range,
ctx->ao_basis.nucleus_max_ang_mom,
ctx->ao_basis.shell_ang_mom,
ctx->ao_basis.shell_prim_index,
ctx->ao_basis.shell_prim_num,
ctx->ao_basis.ao_factor,
ctx->ao_basis.exponent,
ctx->ao_basis.coefficient_normalized,
ctx->ao_basis.ao_vgl);
#else
rc = qmckl_compute_ao_vgl_doc(context,
ctx->ao_basis.ao_num,
ctx->ao_basis.shell_num,
ctx->point.num,
ctx->nucleus.num,
ctx->point.coord.data,
ctx->nucleus.coord.data,
ctx->ao_basis.nucleus_index,
ctx->ao_basis.nucleus_shell_num,
ctx->ao_basis.nucleus_range,
ctx->ao_basis.nucleus_max_ang_mom,
ctx->ao_basis.shell_ang_mom,
ctx->ao_basis.ao_factor,
ctx->ao_basis.shell_vgl,
ctx->ao_basis.ao_vgl);
#endif
if (rc != QMCKL_SUCCESS) {
return rc;
}

View File

@ -412,10 +412,10 @@ qmckl_exit_code qmckl_get_jastrow_aord_num (qmckl_context context, int
qmckl_exit_code qmckl_get_jastrow_bord_num (qmckl_context context, int64_t* const bord_num);
qmckl_exit_code qmckl_get_jastrow_cord_num (qmckl_context context, int64_t* const bord_num);
qmckl_exit_code qmckl_get_jastrow_type_nucl_num (qmckl_context context, int64_t* const type_nucl_num);
qmckl_exit_code qmckl_get_jastrow_type_nucl_vector (qmckl_context context, int64_t* const type_nucl_num, int64_t* size_max);
qmckl_exit_code qmckl_get_jastrow_aord_vector (qmckl_context context, double * const aord_vector, int64_t* size_max);
qmckl_exit_code qmckl_get_jastrow_bord_vector (qmckl_context context, double * const bord_vector, int64_t* size_max);
qmckl_exit_code qmckl_get_jastrow_cord_vector (qmckl_context context, double * const cord_vector, int64_t* size_max);
qmckl_exit_code qmckl_get_jastrow_type_nucl_vector (qmckl_context context, int64_t* const type_nucl_num, const int64_t size_max);
qmckl_exit_code qmckl_get_jastrow_aord_vector (qmckl_context context, double * const aord_vector, const int64_t size_max);
qmckl_exit_code qmckl_get_jastrow_bord_vector (qmckl_context context, double * const bord_vector, const int64_t size_max);
qmckl_exit_code qmckl_get_jastrow_cord_vector (qmckl_context context, double * const cord_vector, const int64_t size_max);
#+end_src
Along with these core functions, calculation of the jastrow factor
@ -474,7 +474,7 @@ qmckl_exit_code qmckl_get_jastrow_aord_num (const qmckl_context context, int64_t
}
assert (ctx->jastrow.aord_num > 0);
*aord_num = ctx->jastrow.aord_num;
,*aord_num = ctx->jastrow.aord_num;
return QMCKL_SUCCESS;
}
@ -501,7 +501,7 @@ qmckl_exit_code qmckl_get_jastrow_bord_num (const qmckl_context context, int64_t
}
assert (ctx->jastrow.bord_num > 0);
*bord_num = ctx->jastrow.bord_num;
,*bord_num = ctx->jastrow.bord_num;
return QMCKL_SUCCESS;
}
@ -528,7 +528,7 @@ qmckl_exit_code qmckl_get_jastrow_cord_num (const qmckl_context context, int64_t
}
assert (ctx->jastrow.cord_num > 0);
*cord_num = ctx->jastrow.cord_num;
,*cord_num = ctx->jastrow.cord_num;
return QMCKL_SUCCESS;
}
@ -555,11 +555,15 @@ qmckl_exit_code qmckl_get_jastrow_type_nucl_num (const qmckl_context context, in
}
assert (ctx->jastrow.type_nucl_num > 0);
*type_nucl_num = ctx->jastrow.type_nucl_num;
,*type_nucl_num = ctx->jastrow.type_nucl_num;
return QMCKL_SUCCESS;
}
qmckl_exit_code qmckl_get_jastrow_type_nucl_vector (const qmckl_context context, int64_t * const type_nucl_vector, int64_t* size_max) {
qmckl_exit_code
qmckl_get_jastrow_type_nucl_vector (const qmckl_context context,
int64_t* const type_nucl_vector,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return (char) 0;
@ -582,12 +586,21 @@ qmckl_exit_code qmckl_get_jastrow_type_nucl_vector (const qmckl_context context,
}
assert (ctx->jastrow.type_nucl_vector != NULL);
if (size_max < ctx->jastrow.type_nucl_num) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_jastrow_type_nucl_vector",
"Array too small. Expected jastrow.type_nucl_num");
}
memcpy(type_nucl_vector, ctx->jastrow.type_nucl_vector, ctx->jastrow.type_nucl_num*sizeof(int64_t));
(*size_max) = ctx->jastrow.type_nucl_num;
return QMCKL_SUCCESS;
}
qmckl_exit_code qmckl_get_jastrow_aord_vector (const qmckl_context context, double * const aord_vector, int64_t* size_max) {
qmckl_exit_code
qmckl_get_jastrow_aord_vector (const qmckl_context context,
double * const aord_vector,
const int64_t size_max) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return (char) 0;
@ -611,12 +624,20 @@ qmckl_exit_code qmckl_get_jastrow_aord_vector (const qmckl_context context, doub
assert (ctx->jastrow.aord_vector != NULL);
int64_t sze = (ctx->jastrow.aord_num + 1)*ctx->jastrow.type_nucl_num;
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_jastrow_aord_vector",
"Array too small. Expected (ctx->jastrow.aord_num + 1)*ctx->jastrow.type_nucl_num");
}
memcpy(aord_vector, ctx->jastrow.aord_vector, sze*sizeof(double));
(*size_max) = sze;
return QMCKL_SUCCESS;
}
qmckl_exit_code qmckl_get_jastrow_bord_vector (const qmckl_context context, double * const bord_vector, int64_t* size_max) {
qmckl_exit_code
qmckl_get_jastrow_bord_vector (const qmckl_context context,
double * const bord_vector,
const int64_t size_max) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return (char) 0;
@ -640,12 +661,20 @@ qmckl_exit_code qmckl_get_jastrow_bord_vector (const qmckl_context context, doub
assert (ctx->jastrow.bord_vector != NULL);
int64_t sze=ctx->jastrow.bord_num +1;
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_jastrow_bord_vector",
"Array too small. Expected (ctx->jastrow.bord_num + 1)");
}
memcpy(bord_vector, ctx->jastrow.bord_vector, sze*sizeof(double));
(*size_max) = sze;
return QMCKL_SUCCESS;
}
qmckl_exit_code qmckl_get_jastrow_cord_vector (const qmckl_context context, double * const cord_vector, int64_t* size_max) {
qmckl_exit_code
qmckl_get_jastrow_cord_vector (const qmckl_context context,
double * const cord_vector,
const int64_t size_max) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return (char) 0;
@ -674,8 +703,13 @@ qmckl_exit_code qmckl_get_jastrow_cord_vector (const qmckl_context context, doub
if (rc != QMCKL_SUCCESS) return rc;
int64_t sze=dim_cord_vect * ctx->jastrow.type_nucl_num;
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_jastrow_cord_vector",
"Array too small. Expected dim_cord_vect * jastrow.type_nucl_num");
}
memcpy(cord_vector, ctx->jastrow.cord_vector, sze*sizeof(double));
(*size_max) = sze;
return QMCKL_SUCCESS;
}
@ -690,9 +724,9 @@ qmckl_exit_code qmckl_get_jastrow_cord_vector (const qmckl_context context, doub
qmckl_exit_code qmckl_set_jastrow_ord_num (qmckl_context context, const int64_t aord_num, const int64_t bord_num, const int64_t cord_num);
qmckl_exit_code qmckl_set_jastrow_type_nucl_num (qmckl_context context, const int64_t type_nucl_num);
qmckl_exit_code qmckl_set_jastrow_type_nucl_vector (qmckl_context context, const int64_t* type_nucl_vector, const int64_t nucl_num);
qmckl_exit_code qmckl_set_jastrow_aord_vector (qmckl_context context, const double * aord_vector, int64_t size_max);
qmckl_exit_code qmckl_set_jastrow_bord_vector (qmckl_context context, const double * bord_vector, int64_t size_max);
qmckl_exit_code qmckl_set_jastrow_cord_vector (qmckl_context context, const double * cord_vector, int64_t size_max);
qmckl_exit_code qmckl_set_jastrow_aord_vector (qmckl_context context, const double * aord_vector, const int64_t size_max);
qmckl_exit_code qmckl_set_jastrow_bord_vector (qmckl_context context, const double * bord_vector, const int64_t size_max);
qmckl_exit_code qmckl_set_jastrow_cord_vector (qmckl_context context, const double * cord_vector, const int64_t size_max);
#+end_src
#+NAME:pre2
@ -718,7 +752,12 @@ return QMCKL_SUCCESS;
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_set_jastrow_ord_num(qmckl_context context, const int64_t aord_num, const int64_t bord_num, const int64_t cord_num) {
qmckl_exit_code
qmckl_set_jastrow_ord_num(qmckl_context context,
const int64_t aord_num,
const int64_t bord_num,
const int64_t cord_num)
{
<<pre2>>
if (aord_num <= 0) {
@ -750,7 +789,10 @@ qmckl_exit_code qmckl_set_jastrow_ord_num(qmckl_context context, const int64_t a
<<post2>>
}
qmckl_exit_code qmckl_set_jastrow_type_nucl_num(qmckl_context context, const int64_t type_nucl_num) {
qmckl_exit_code
qmckl_set_jastrow_type_nucl_num(qmckl_context context, const int64_t type_nucl_num)
{
<<pre2>>
if (type_nucl_num <= 0) {
@ -766,7 +808,12 @@ qmckl_exit_code qmckl_set_jastrow_type_nucl_num(qmckl_context context, const int
<<post2>>
}
qmckl_exit_code qmckl_set_jastrow_type_nucl_vector(qmckl_context context, int64_t const * type_nucl_vector, const int64_t nucl_num) {
qmckl_exit_code
qmckl_set_jastrow_type_nucl_vector(qmckl_context context,
int64_t const * type_nucl_vector,
const int64_t nucl_num)
{
<<pre2>>
int32_t mask = 1 << 2;
@ -816,7 +863,12 @@ qmckl_exit_code qmckl_set_jastrow_type_nucl_vector(qmckl_context context, int64_
<<post2>>
}
qmckl_exit_code qmckl_set_jastrow_aord_vector(qmckl_context context, double const * aord_vector, int64_t size_max) {
qmckl_exit_code
qmckl_set_jastrow_aord_vector(qmckl_context context,
double const * aord_vector,
const int64_t size_max)
{
<<pre2>>
int32_t mask = 1 << 3;
@ -849,7 +901,7 @@ qmckl_exit_code qmckl_set_jastrow_aord_vector(qmckl_context context, double cons
return qmckl_failwith( context, rc,
"qmckl_set_ord_vector",
NULL);
}
}
}
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
@ -878,7 +930,12 @@ qmckl_exit_code qmckl_set_jastrow_aord_vector(qmckl_context context, double cons
<<post2>>
}
qmckl_exit_code qmckl_set_jastrow_bord_vector(qmckl_context context, double const * bord_vector, int64_t size_max) {
qmckl_exit_code
qmckl_set_jastrow_bord_vector(qmckl_context context,
double const * bord_vector,
const int64_t size_max)
{
<<pre2>>
int32_t mask = 1 << 4;
@ -936,7 +993,12 @@ qmckl_exit_code qmckl_set_jastrow_bord_vector(qmckl_context context, double cons
<<post2>>
}
qmckl_exit_code qmckl_set_jastrow_cord_vector(qmckl_context context, double const * cord_vector, int64_t size_max) {
qmckl_exit_code
qmckl_set_jastrow_cord_vector(qmckl_context context,
double const * cord_vector,
const int64_t size_max)
{
<<pre2>>
int32_t mask = 1 << 5;
@ -1069,6 +1131,7 @@ double* elec_coord = &(n2_elec_coord[0][0][0]);
const double* nucl_charge = n2_charge;
int64_t nucl_num = n2_nucl_num;
double* nucl_coord = &(n2_nucl_coord[0][0]);
int64_t size_max;
/* Provide Electron data */
@ -1246,11 +1309,17 @@ assert(qmckl_nucleus_provided(context));
*** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_jastrow_asymp_jasb(qmckl_context context, double* const asymp_jasb, int64_t* size_max);
qmckl_exit_code
qmckl_get_jastrow_asymp_jasb(qmckl_context context,
double* const asymp_jasb,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_jastrow_asymp_jasb(qmckl_context context, double* const asymp_jasb, int64_t* size_max)
qmckl_exit_code
qmckl_get_jastrow_asymp_jasb(qmckl_context context,
double* const asymp_jasb,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
@ -1265,8 +1334,13 @@ qmckl_exit_code qmckl_get_jastrow_asymp_jasb(qmckl_context context, double* cons
assert (ctx != NULL);
size_t sze = 2;
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_jastrow_asymp_jasb",
"Array too small. Expected 2");
}
memcpy(asymp_jasb, ctx->jastrow.asymp_jasb, sze * sizeof(double));
(*size_max) = sze;
return QMCKL_SUCCESS;
}
@ -1500,8 +1574,7 @@ assert(rc == QMCKL_SUCCESS);
assert(qmckl_jastrow_provided(context));
double asymp_jasb[2];
int64_t size_max=0;
rc = qmckl_get_jastrow_asymp_jasb(context, asymp_jasb,&size_max);
rc = qmckl_get_jastrow_asymp_jasb(context, asymp_jasb,2);
// calculate asymp_jasb
assert(fabs(asymp_jasb[0]-0.5323750557252571) < 1.e-12);
@ -1521,11 +1594,17 @@ f_{ee} = \sum_{i,j<i} \left\{ \frac{ \eta B_0 C_{ij}}{1 - B_1 C_{ij}} - J_{asym
*** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_jastrow_factor_ee(qmckl_context context, double* const factor_ee, int64_t* size_max);
qmckl_exit_code
qmckl_get_jastrow_factor_ee(qmckl_context context,
double* const factor_ee,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_jastrow_factor_ee(qmckl_context context, double* const factor_ee, int64_t* size_max)
qmckl_exit_code
qmckl_get_jastrow_factor_ee(qmckl_context context,
double* const factor_ee,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
@ -1540,8 +1619,13 @@ qmckl_exit_code qmckl_get_jastrow_factor_ee(qmckl_context context, double* const
assert (ctx != NULL);
int64_t sze=ctx->electron.walk_num;
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_jastrow_factor_ee",
"Array too small. Expected walk_num");
}
memcpy(factor_ee, ctx->jastrow.factor_ee, sze*sizeof(double));
(*size_max) = sze;
return QMCKL_SUCCESS;
}
@ -1804,8 +1888,7 @@ print("factor_ee :",factor_ee)
assert(qmckl_jastrow_provided(context));
double factor_ee[walk_num];
size_max=0;
rc = qmckl_get_jastrow_factor_ee(context, factor_ee, &size_max);
rc = qmckl_get_jastrow_factor_ee(context, factor_ee, walk_num);
// calculate factor_ee
assert(fabs(factor_ee[0]+4.282760865958113) < 1.e-12);
@ -1824,11 +1907,17 @@ assert(fabs(factor_ee[0]+4.282760865958113) < 1.e-12);
*** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_jastrow_factor_ee_deriv_e(qmckl_context context, double* const factor_ee_deriv_e, int64_t* size_max);
qmckl_exit_code
qmckl_get_jastrow_factor_ee_deriv_e(qmckl_context context,
double* const factor_ee_deriv_e,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_jastrow_factor_ee_deriv_e(qmckl_context context, double* const factor_ee_deriv_e, int64_t* size_max)
qmckl_exit_code
qmckl_get_jastrow_factor_ee_deriv_e(qmckl_context context,
double* const factor_ee_deriv_e,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
@ -1843,8 +1932,14 @@ qmckl_exit_code qmckl_get_jastrow_factor_ee_deriv_e(qmckl_context context, doubl
assert (ctx != NULL);
int64_t sze = ctx->electron.walk_num * 4 * ctx->electron.num;
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_jastrow_factor_ee_deriv_e",
"Array too small. Expected 4*walk_num*elec_num");
}
memcpy(factor_ee_deriv_e, ctx->jastrow.factor_ee_deriv_e, sze * sizeof(double));
(*size_max) = sze;
return QMCKL_SUCCESS;
}
@ -2217,8 +2312,7 @@ assert(qmckl_jastrow_provided(context));
// calculate factor_ee_deriv_e
double factor_ee_deriv_e[walk_num][4][elec_num];
size_max=0;
rc = qmckl_get_jastrow_factor_ee_deriv_e(context, &(factor_ee_deriv_e[0][0][0]),&size_max);
rc = qmckl_get_jastrow_factor_ee_deriv_e(context, &(factor_ee_deriv_e[0][0][0]),walk_num*4*elec_num);
// check factor_ee_deriv_e
assert(fabs(factor_ee_deriv_e[0][0][0]-0.16364894652107934) < 1.e-12);
@ -2240,11 +2334,17 @@ f_{en} = \sum_{i,j<i} \left\{ \frac{ A_0 C_{ij}}{1 - A_1 C_{ij}} + \sum^{nord}_{
*** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_jastrow_factor_en(qmckl_context context, double* const factor_en, int64_t* size_max);
qmckl_exit_code
qmckl_get_jastrow_factor_en(qmckl_context context,
double* const factor_en,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_jastrow_factor_en(qmckl_context context, double* const factor_en, int64_t* size_max)
qmckl_exit_code
qmckl_get_jastrow_factor_en(qmckl_context context,
double* const factor_en,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
@ -2259,8 +2359,13 @@ qmckl_exit_code qmckl_get_jastrow_factor_en(qmckl_context context, double* const
assert (ctx != NULL);
int64_t sze=ctx->electron.walk_num;
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_jastrow_factor_en",
"Array too small. Expected walk_num");
}
memcpy(factor_en, ctx->jastrow.factor_en, sze*sizeof(double));
(*size_max)=sze;
return QMCKL_SUCCESS;
}
@ -2517,8 +2622,7 @@ print("factor_en :",factor_en)
assert(qmckl_jastrow_provided(context));
double factor_en[walk_num];
size_max=0;
rc = qmckl_get_jastrow_factor_en(context, factor_en,&size_max);
rc = qmckl_get_jastrow_factor_en(context, factor_en,walk_num);
// calculate factor_en
assert(fabs(factor_en[0]+5.865822569188727) < 1.e-12);
@ -2534,11 +2638,17 @@ assert(fabs(factor_en[0]+5.865822569188727) < 1.e-12);
*** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_jastrow_factor_en_deriv_e(qmckl_context context, double* const factor_en_deriv_e, int64_t* size_max);
qmckl_exit_code
qmckl_get_jastrow_factor_en_deriv_e(qmckl_context context,
double* const factor_en_deriv_e,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_jastrow_factor_en_deriv_e(qmckl_context context, double* const factor_en_deriv_e, int64_t* size_max)
qmckl_exit_code
qmckl_get_jastrow_factor_en_deriv_e(qmckl_context context,
double* const factor_en_deriv_e,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
@ -2553,8 +2663,13 @@ qmckl_exit_code qmckl_get_jastrow_factor_en_deriv_e(qmckl_context context, doubl
assert (ctx != NULL);
int64_t sze = ctx->electron.walk_num * 4 * ctx->electron.num;
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_jastrow_factor_en_deriv_e",
"Array too small. Expected 4*walk_num*elec_num");
}
memcpy(factor_en_deriv_e, ctx->jastrow.factor_en_deriv_e, sze*sizeof(double));
(*size_max) = sze;
return QMCKL_SUCCESS;
}
@ -2921,8 +3036,7 @@ assert(qmckl_jastrow_provided(context));
// calculate factor_en_deriv_e
double factor_en_deriv_e[walk_num][4][elec_num];
size_max=0;
rc = qmckl_get_jastrow_factor_en_deriv_e(context, &(factor_en_deriv_e[0][0][0]),&size_max);
rc = qmckl_get_jastrow_factor_en_deriv_e(context, &(factor_en_deriv_e[0][0][0]),walk_num*4*elec_num);
// check factor_en_deriv_e
assert(fabs(factor_en_deriv_e[0][0][0]-0.11609919541763383) < 1.e-12);
@ -2946,11 +3060,17 @@ assert(fabs(factor_en_deriv_e[0][3][0]+0.9667363412285741 ) < 1.e-12);
*** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_jastrow_een_rescaled_e(qmckl_context context, double* const distance_rescaled, int64_t* size_max);
qmckl_exit_code
qmckl_get_jastrow_een_rescaled_e(qmckl_context context,
double* const distance_rescaled,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_jastrow_een_rescaled_e(qmckl_context context, double* const distance_rescaled, int64_t* size_max)
qmckl_exit_code
qmckl_get_jastrow_een_rescaled_e(qmckl_context context,
double* const distance_rescaled,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
@ -2965,8 +3085,13 @@ qmckl_exit_code qmckl_get_jastrow_een_rescaled_e(qmckl_context context, double*
assert (ctx != NULL);
size_t sze = ctx->electron.num * ctx->electron.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1);
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_jastrow_factor_een_rescaled_e",
"Array too small. Expected ctx->electron.num * ctx->electron.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1)");
}
memcpy(distance_rescaled, ctx->jastrow.een_rescaled_e, sze * sizeof(double));
(*size_max) = sze;
return QMCKL_SUCCESS;
}
@ -3242,8 +3367,7 @@ assert(qmckl_electron_provided(context));
double een_rescaled_e[walk_num][(cord_num + 1)][elec_num][elec_num];
size_max=0;
rc = qmckl_get_jastrow_een_rescaled_e(context, &(een_rescaled_e[0][0][0][0]),&size_max);
rc = qmckl_get_jastrow_een_rescaled_e(context, &(een_rescaled_e[0][0][0][0]),elec_num*elec_num*(cord_num+1)*walk_num);
// value of (0,2,1)
assert(fabs(een_rescaled_e[0][1][0][2]-0.08084493981483197) < 1.e-12);
@ -3268,11 +3392,17 @@ assert(fabs(een_rescaled_e[0][2][1][5]-0.3424402276009091) < 1.e-12);
*** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_jastrow_een_rescaled_e_deriv_e(qmckl_context context, double* const distance_rescaled, int64_t* size_max);
qmckl_exit_code
qmckl_get_jastrow_een_rescaled_e_deriv_e(qmckl_context context,
double* const distance_rescaled,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_jastrow_een_rescaled_e_deriv_e(qmckl_context context, double* const distance_rescaled, int64_t* size_max)
qmckl_exit_code
qmckl_get_jastrow_een_rescaled_e_deriv_e(qmckl_context context,
double* const distance_rescaled,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
@ -3287,8 +3417,13 @@ qmckl_exit_code qmckl_get_jastrow_een_rescaled_e_deriv_e(qmckl_context context,
assert (ctx != NULL);
size_t sze = ctx->electron.num * 4 * ctx->electron.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1);
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_jastrow_factor_een_deriv_e",
"Array too small. Expected ctx->electron.num * 4 * ctx->electron.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1)");
}
memcpy(distance_rescaled, ctx->jastrow.een_rescaled_e_deriv_e, sze * sizeof(double));
(*size_max) = sze;
return QMCKL_SUCCESS;
}
@ -3609,8 +3744,9 @@ for l in range(0,cord_num+1):
#+begin_src c :tangle (eval c_test)
//assert(qmckl_electron_provided(context));
double een_rescaled_e_deriv_e[walk_num][(cord_num + 1)][elec_num][4][elec_num];
size_max=0;
rc = qmckl_get_jastrow_een_rescaled_e_deriv_e(context, &(een_rescaled_e_deriv_e[0][0][0][0][0]),&size_max);
size_max=walk_num*(cord_num + 1)*elec_num*4*elec_num;
rc = qmckl_get_jastrow_een_rescaled_e_deriv_e(context,
&(een_rescaled_e_deriv_e[0][0][0][0][0]),size_max);
// value of (0,0,0,2,1)
assert(fabs(een_rescaled_e_deriv_e[0][1][0][0][2] + 0.05991352796887283 ) < 1.e-12);
@ -3635,11 +3771,17 @@ assert(fabs(een_rescaled_e_deriv_e[0][2][1][0][5] + 0.5880599146214673 ) < 1.
*** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_jastrow_een_rescaled_n(qmckl_context context, double* const distance_rescaled, int64_t* size_max);
qmckl_exit_code
qmckl_get_jastrow_een_rescaled_n(qmckl_context context,
double* const distance_rescaled,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_jastrow_een_rescaled_n(qmckl_context context, double* const distance_rescaled, int64_t* size_max)
qmckl_exit_code
qmckl_get_jastrow_een_rescaled_n(qmckl_context context,
double* const distance_rescaled,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
@ -3654,8 +3796,13 @@ qmckl_exit_code qmckl_get_jastrow_een_rescaled_n(qmckl_context context, double*
assert (ctx != NULL);
size_t sze = ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1);
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_jastrow_factor_een_deriv_e",
"Array too small. Expected ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1)");
}
memcpy(distance_rescaled, ctx->jastrow.een_rescaled_n, sze * sizeof(double));
(*size_max)=sze;
return QMCKL_SUCCESS;
}
@ -3914,8 +4061,8 @@ print(" een_rescaled_n[1, 5, 2] = ",een_rescaled_n[1, 5, 2])
assert(qmckl_electron_provided(context));
double een_rescaled_n[walk_num][(cord_num + 1)][nucl_num][elec_num];
size_max=0;
rc = qmckl_get_jastrow_een_rescaled_n(context, &(een_rescaled_n[0][0][0][0]),&size_max);
size_max=walk_num*(cord_num + 1)*nucl_num*elec_num;
rc = qmckl_get_jastrow_een_rescaled_n(context, &(een_rescaled_n[0][0][0][0]),size_max);
// value of (0,2,1)
assert(fabs(een_rescaled_n[0][1][0][2]-0.10612983920006765) < 1.e-12);
@ -3936,11 +4083,17 @@ assert(fabs(een_rescaled_n[0][2][1][5]-0.01343938025140174) < 1.e-12);
*** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_jastrow_een_rescaled_n_deriv_e(qmckl_context context, double* const distance_rescaled, int64_t* size_max);
qmckl_exit_code
qmckl_get_jastrow_een_rescaled_n_deriv_e(qmckl_context context,
double* const distance_rescaled,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_jastrow_een_rescaled_n_deriv_e(qmckl_context context, double* const distance_rescaled, int64_t* size_max)
qmckl_exit_code
qmckl_get_jastrow_een_rescaled_n_deriv_e(qmckl_context context,
double* const distance_rescaled,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
@ -3955,8 +4108,13 @@ qmckl_exit_code qmckl_get_jastrow_een_rescaled_n_deriv_e(qmckl_context context,
assert (ctx != NULL);
size_t sze = ctx->electron.num * 4 * ctx->nucleus.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1);
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_jastrow_factor_een_deriv_e",
"Array too small. Expected ctx->electron.num * 4 * ctx->nucleus.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1)");
}
memcpy(distance_rescaled, ctx->jastrow.een_rescaled_n_deriv_e, sze * sizeof(double));
(*size_max)=sze;
return QMCKL_SUCCESS;
}
@ -4287,8 +4445,8 @@ print(" een_rescaled_n_deriv_e[2, 1, 6, 2] = ",een_rescaled_n_deriv_e[5, 0, 1, 2
assert(qmckl_electron_provided(context));
double een_rescaled_n_deriv_e[walk_num][(cord_num + 1)][nucl_num][4][elec_num];
size_max=0;
rc = qmckl_get_jastrow_een_rescaled_n_deriv_e(context, &(een_rescaled_n_deriv_e[0][0][0][0][0]),&size_max);
size_max=walk_num*(cord_num + 1)*nucl_num*4*elec_num;
rc = qmckl_get_jastrow_een_rescaled_n_deriv_e(context, &(een_rescaled_n_deriv_e[0][0][0][0][0]),size_max);
// value of (0,2,1)
assert(fabs(een_rescaled_n_deriv_e[0][1][0][0][2]+0.07633444246999128 ) < 1.e-12);
@ -5330,11 +5488,17 @@ assert(fabs(dtmp_c[0][1][0][0][0][0] - 0.237440520852232) < 1e-12);
*** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_jastrow_factor_een(qmckl_context context, double* const factor_een, int64_t* size_max);
qmckl_exit_code
qmckl_get_jastrow_factor_een(qmckl_context context,
double* const factor_een,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_jastrow_factor_een(qmckl_context context, double* const factor_een, int64_t* size_max)
qmckl_exit_code
qmckl_get_jastrow_factor_een(qmckl_context context,
double* const factor_een,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
@ -5348,9 +5512,14 @@ qmckl_exit_code qmckl_get_jastrow_factor_een(qmckl_context context, double* cons
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
int64_t sze = ctx->electron.walk_num * ctx->electron.num;
int64_t sze = ctx->electron.walk_num;
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_jastrow_factor_een",
"Array too small. Expected walk_num");
}
memcpy(factor_een, ctx->jastrow.factor_een, sze*sizeof(double));
(*size_max)=sze;
return QMCKL_SUCCESS;
}
@ -5804,8 +5973,7 @@ print("factor_een:",factor_een)
assert(qmckl_jastrow_provided(context));
double factor_een[walk_num];
size_max=0;
rc = qmckl_get_jastrow_factor_een(context, &(factor_een[0]),&size_max);
rc = qmckl_get_jastrow_factor_een(context, &(factor_een[0]),walk_num);
assert(fabs(factor_een[0] + 0.37407972141304213) < 1e-12);
#+end_src
@ -5819,11 +5987,17 @@ assert(fabs(factor_een[0] + 0.37407972141304213) < 1e-12);
*** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_jastrow_factor_een_deriv_e(qmckl_context context, double* const factor_een_deriv_e, int64_t* size_max);
qmckl_exit_code
qmckl_get_jastrow_factor_een_deriv_e(qmckl_context context,
double* const factor_een_deriv_e,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_jastrow_factor_een_deriv_e(qmckl_context context, double* const factor_een_deriv_e, int64_t* size_max)
qmckl_exit_code
qmckl_get_jastrow_factor_een_deriv_e(qmckl_context context,
double* const factor_een_deriv_e,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
@ -5837,9 +6011,14 @@ qmckl_exit_code qmckl_get_jastrow_factor_een_deriv_e(qmckl_context context, doub
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
int64_t sze = ctx->electron.walk_num * ctx->electron.num;
int64_t sze = ctx->electron.walk_num * 4 * ctx->electron.num;
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_jastrow_factor_een_deriv_e",
"Array too small. Expected 4*walk_num*elec_num");
}
memcpy(factor_een_deriv_e, ctx->jastrow.factor_een_deriv_e, sze*sizeof(double));
(*size_max)=sze;
return QMCKL_SUCCESS;
}
@ -6384,11 +6563,10 @@ print("factor_een:",factor_een)
/* Check if Jastrow is properly initialized */
assert(qmckl_jastrow_provided(context));
double factor_een_deriv_e[walk_num][elec_num];
size_max=0;
rc = qmckl_get_jastrow_factor_een_deriv_e(context, &(factor_een_deriv_e[0][0]),&size_max);
double factor_een_deriv_e[4][walk_num][elec_num];
rc = qmckl_get_jastrow_factor_een_deriv_e(context, &(factor_een_deriv_e[0][0][0]),4*walk_num*elec_num);
assert(fabs(factor_een_deriv_e[0][0] + 0.0005481671107226865) < 1e-12);
assert(fabs(factor_een_deriv_e[0][0][0] + 0.0005481671107226865) < 1e-12);
#+end_src
* End of files :noexport: