diff --git a/org/qmckl_ao.org b/org/qmckl_ao.org index d40e478..8783884 100644 --- a/org/qmckl_ao.org +++ b/org/qmckl_ao.org @@ -3387,6 +3387,7 @@ for (j=0 ; j 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) @@ -4710,21 +4760,37 @@ integer function qmckl_compute_ao_vgl_hpc_f(context, & integer :: l, il, k integer*8 :: ipoint, inucl, ishell integer*8 :: ishell_start, ishell_end - integer :: lstart(0:20) - double precision :: x, y, z, r2 - double precision :: cutoff + integer :: lstart(0:20) + 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,68 +4824,80 @@ 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 - + deallocate(poly_vgl, powers) 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; } diff --git a/org/qmckl_jastrow.org b/org/qmckl_jastrow.org index abeac12..d745dcf 100644 --- a/org/qmckl_jastrow.org +++ b/org/qmckl_jastrow.org @@ -9,10 +9,10 @@ nuclear ($\mathbf{R}$) coordinates. Its defined as $\exp(J(\mathbf{r},\mathbf{R}))$, where \[ - J(\mathbf{r},\mathbf{R}) = J_{\text{eN}}(\mathbf{r},\mathbf{R}) + J_{\text{ee}}(\mathbf{r}) + J_{\text{eeN}}(\mathbf{r},\mathbf{R}) + J(\mathbf{r},\mathbf{R}) = J_{\text{eN}}(\mathbf{r},\mathbf{R}) + J_{\text{ee}}(\mathbf{r}) + J_{\text{eeN}}(\mathbf{r},\mathbf{R}) \] - In the following, we us the notations $r_{ij} = |\mathbf{r}_i - \mathbf{r}_j|$ and + In the following, we us the notations $r_{ij} = |\mathbf{r}_i - \mathbf{r}_j|$ and $R_{i\alpha} = |\mathbf{r}_i - \mathbf{R}_\alpha|$. $J_{\text{eN}}$ contains electron-nucleus terms: @@ -23,9 +23,9 @@ \sum_{p=2}^{N_\text{ord}^a} a_{p+1}\, [f(R_{i\alpha})]^p - J_{eN}^\infty \] - $J_{\text{ee}}$ contains electron-electron terms: + $J_{\text{ee}}$ contains electron-electron terms: \[ - J_{\text{ee}}(\mathbf{r}) = + J_{\text{ee}}(\mathbf{r}) = \sum_{i=1}^{N_\text{elec}} \sum_{j=1}^{i-1} \frac{b_1\, f(r_{ij})}{1+b_2\, f(r_{ij})} + \sum_{p=2}^{N_\text{ord}^b} a_{p+1}\, [f(r_{ij})]^p - J_{ee}^\infty @@ -43,7 +43,7 @@ \sum_{l=0}^{p-k-2\delta_{k,0}} c_{lkp\alpha} \left[ g({r}_{ij}) \right]^k \left[ \left[ g({R}_{i\alpha}) \right]^l + \left[ g({R}_{j\alpha}) \right]^l \right] - \left[ g({R}_{i\,\alpha}) \, g({R}_{j\alpha}) \right]^{(p-k-l)/2} + \left[ g({R}_{i\,\alpha}) \, g({R}_{j\alpha}) \right]^{(p-k-l)/2} \] $c_{lkp\alpha}$ are non-zero only when $p-k-l$ is even. @@ -55,10 +55,10 @@ g(r) = e^{-\kappa\, r}. \] - The terms $J_{\text{ee}}^\infty$ and $J_{\text{eN}}^\infty$ are shifts to ensure that + The terms $J_{\text{ee}}^\infty$ and $J_{\text{eN}}^\infty$ are shifts to ensure that $J_{\text{ee}}$ and $J_{\text{eN}}$ have an asymptotic value of zero. - - + + * Headers :noexport: #+begin_src elisp :noexport :results none (org-babel-lob-ingest "../tools/lib.org") @@ -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; @@ -668,14 +697,19 @@ qmckl_exit_code qmckl_get_jastrow_cord_vector (const qmckl_context context, doub } assert (ctx->jastrow.cord_vector != NULL); - + int64_t dim_cord_vect; qmckl_exit_code rc = qmckl_get_jastrow_dim_cord_vect(context, &dim_cord_vect); 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) +{ <> if (aord_num <= 0) { @@ -750,7 +789,10 @@ qmckl_exit_code qmckl_set_jastrow_ord_num(qmckl_context context, const int64_t a <> } -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) +{ <> if (type_nucl_num <= 0) { @@ -766,7 +808,12 @@ qmckl_exit_code qmckl_set_jastrow_type_nucl_num(qmckl_context context, const int <> } -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) +{ <> int32_t mask = 1 << 2; @@ -816,7 +863,12 @@ qmckl_exit_code qmckl_set_jastrow_type_nucl_vector(qmckl_context context, int64_ <> } -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) +{ <> 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 <> } -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) +{ <> int32_t mask = 1 << 4; @@ -936,7 +993,12 @@ qmckl_exit_code qmckl_set_jastrow_bord_vector(qmckl_context context, double cons <> } -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) +{ <> 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; } @@ -1409,7 +1483,7 @@ qmckl_exit_code qmckl_compute_asymp_jasb ( } asym_one = bord_vector[0] * kappa_inv / (1.0 + bord_vector[1] * kappa_inv); - asymp_jasb[0] = asym_one; + asymp_jasb[0] = asym_one; asymp_jasb[1] = 0.5 * asym_one; for (int i = 0 ; i <= 1; ++i) { @@ -1418,7 +1492,7 @@ qmckl_exit_code qmckl_compute_asymp_jasb ( x = x * kappa_inv; asymp_jasb[i] = asymp_jasb[i] + bord_vector[p + 1] * x; } - } + } return QMCKL_SUCCESS; } @@ -1433,7 +1507,7 @@ qmckl_exit_code qmckl_compute_asymp_jasb ( const int64_t bord_num, const double* bord_vector, const double rescale_factor_kappa_ee, - double* const asymp_jasb ); + double* const asymp_jasb ); #+end_src @@ -1457,7 +1531,7 @@ print("asym_one : ", asym_one) print("asymp_jasb[0] : ", asymp_jasb[0]) print("asymp_jasb[1] : ", asymp_jasb[1]) #+end_src - + #+RESULTS: asymp_jasb : asym_one : 0.43340325572525706 : asymp_jasb[0] : 0.5323750557252571 @@ -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,jelectron.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; } @@ -2002,7 +2097,7 @@ integer function qmckl_compute_factor_ee_deriv_e_f(context, walk_num, elec_num, dx(1) = ee_distance_rescaled_deriv_e(1, i, j, nw) dx(2) = ee_distance_rescaled_deriv_e(2, i, j, nw) dx(3) = ee_distance_rescaled_deriv_e(3, i, j, nw) - dx(4) = ee_distance_rescaled_deriv_e(4, i, j, nw) + dx(4) = ee_distance_rescaled_deriv_e(4, i, j, nw) if((i .LE. up_num .AND. j .LE. up_num ) .OR. & (i .GT. up_num .AND. j .GT. up_num)) then @@ -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,jelectron.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; } @@ -3178,7 +3303,7 @@ end function qmckl_compute_een_rescaled_e_f #+end_src *** Test - + #+begin_src python :results output :exports none :noweb yes import numpy as np @@ -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; } @@ -4215,9 +4373,9 @@ end function qmckl_compute_factor_een_rescaled_n_deriv_e_f end function qmckl_compute_factor_een_rescaled_n_deriv_e #+end_src - + *** Test - + #+begin_src python :results output :exports none :noweb yes import numpy as np @@ -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; } @@ -6327,7 +6506,7 @@ end function qmckl_compute_factor_een_deriv_e_f end function qmckl_compute_factor_een_deriv_e #+end_src - + *** Test #+begin_src python :results output :exports none :noweb yes import numpy as np @@ -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: