From 574ce713d25c519580fcc2e11047f5c35e8c4d16 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 23 Sep 2021 11:16:54 +0200 Subject: [PATCH] Removed walk_num in AOs --- org/qmckl_ao.org | 698 +++++++++++++++++++++-------------------------- 1 file changed, 317 insertions(+), 381 deletions(-) diff --git a/org/qmckl_ao.org b/org/qmckl_ao.org index e913c0e..9ce7a3a 100644 --- a/org/qmckl_ao.org +++ b/org/qmckl_ao.org @@ -118,23 +118,23 @@ int main() { Computed data: - |--------------------------+--------------------------------------+-----------------------------------------------------------------------------------------------| - | ~coefficient_normalized~ | ~[prim_num]~ | Normalized primitive coefficients | - | ~nucleus_prim_index~ | ~[nucl_num]~ | Index of the first primitive for each nucleus | - | ~nucleus_max_ang_mom~ | ~[nucl_num]~ | Maximum angular momentum for each nucleus | - | ~nucleus_range~ | ~[nucl_num]~ | Distance beyond which all the AOs are zero | - |--------------------------+--------------------------------------+-----------------------------------------------------------------------------------------------| - | ~primitive_vgl~ | ~[5][walk_num][elec_num][prim_num]~ | Value, gradients, Laplacian of the primitives at electron positions | - | ~primitive_vgl_date~ | ~uint64_t~ | Late modification date of Value, gradients, Laplacian of the primitives at electron positions | - | ~shell_vgl~ | ~[5][walk_num][elec_num][shell_num]~ | Value, gradients, Laplacian of the primitives at electron positions | - | ~shell_vgl_date~ | ~uint64_t~ | Late modification date of Value, gradients, Laplacian of the AOs at electron positions | - | ~ao_vgl~ | ~[5][walk_num][elec_num][ao_num]~ | Value, gradients, Laplacian of the primitives at electron positions | - | ~ao_vgl_date~ | ~uint64_t~ | Late modification date of Value, gradients, Laplacian of the AOs at electron positions | - |--------------------------+--------------------------------------+-----------------------------------------------------------------------------------------------| - | ~nucl_shell_index~ | ~[nucl_num]~ | Index of the first shell for each nucleus | - | ~exponent_sorted~ | ~[prim_num]~ | Array of exponents for sorted primitives | - | ~coeff_norm_sorted~ | ~[prim_num]~ | Array of normalized coefficients for sorted primitives | - | ~prim_factor_sorted~ | ~[prim_num]~ | Normalization factors of the sorted primtives | + |--------------------------+----------------------------+-----------------------------------------------------------------------------------------------| + | ~coefficient_normalized~ | ~[prim_num]~ | Normalized primitive coefficients | + | ~nucleus_prim_index~ | ~[nucl_num]~ | Index of the first primitive for each nucleus | + | ~nucleus_max_ang_mom~ | ~[nucl_num]~ | Maximum angular momentum for each nucleus | + | ~nucleus_range~ | ~[nucl_num]~ | Distance beyond which all the AOs are zero | + |--------------------------+----------------------------+-----------------------------------------------------------------------------------------------| + | ~primitive_vgl~ | ~[5][elec_num][prim_num]~ | Value, gradients, Laplacian of the primitives at electron positions | + | ~primitive_vgl_date~ | ~uint64_t~ | Late modification date of Value, gradients, Laplacian of the primitives at electron positions | + | ~shell_vgl~ | ~[5][elec_num][shell_num]~ | Value, gradients, Laplacian of the primitives at electron positions | + | ~shell_vgl_date~ | ~uint64_t~ | Late modification date of Value, gradients, Laplacian of the AOs at electron positions | + | ~ao_vgl~ | ~[5][elec_num][ao_num]~ | Value, gradients, Laplacian of the primitives at electron positions | + | ~ao_vgl_date~ | ~uint64_t~ | Late modification date of Value, gradients, Laplacian of the AOs at electron positions | + |--------------------------+----------------------------+-----------------------------------------------------------------------------------------------| + | ~nucl_shell_index~ | ~[nucl_num]~ | Index of the first shell for each nucleus | + | ~exponent_sorted~ | ~[prim_num]~ | Array of exponents for sorted primitives | + | ~coeff_norm_sorted~ | ~[prim_num]~ | Array of normalized coefficients for sorted primitives | + | ~prim_factor_sorted~ | ~[prim_num]~ | Normalization factors of the sorted primtives | For H_2 with the following basis set, @@ -1837,7 +1837,7 @@ qmckl_exit_code qmckl_get_ao_basis_primitive_vgl(qmckl_context context, double* qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - size_t sze = ctx->ao_basis.prim_num * 5 * ctx->electron.num * ctx->electron.walk_num; + size_t sze = ctx->ao_basis.prim_num * 5 * ctx->electron.num; memcpy(primitive_vgl, ctx->ao_basis.primitive_vgl, sze * sizeof(double)); return QMCKL_SUCCESS; @@ -1876,7 +1876,7 @@ qmckl_exit_code qmckl_provide_ao_basis_primitive_vgl(qmckl_context context) qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = ctx->ao_basis.prim_num * 5 * ctx->electron.num * - ctx->electron.walk_num * sizeof(double); + sizeof(double); double* primitive_vgl = (double*) qmckl_malloc(context, mem_info); if (primitive_vgl == NULL) { @@ -1894,7 +1894,6 @@ qmckl_exit_code qmckl_provide_ao_basis_primitive_vgl(qmckl_context context) ctx->ao_basis.prim_num, ctx->electron.num, ctx->nucleus.num, - ctx->electron.walk_num, ctx->ao_basis.nucleus_prim_index, ctx->electron.coord_new, ctx->nucleus.coord, @@ -1925,20 +1924,19 @@ qmckl_exit_code qmckl_provide_ao_basis_primitive_vgl(qmckl_context context) :END: #+NAME: qmckl_ao_basis_primitive_gaussian_vgl_args - | qmckl_context | context | in | Global state | - | int64_t | prim_num | in | Number of primitives | - | int64_t | elec_num | in | Number of electrons | - | int64_t | nucl_num | in | Number of nuclei | - | int64_t | walk_num | in | Number of walkers | - | int64_t | nucleus_prim_index[nucl_num] | in | Index of the 1st primitive of each nucleus | - | double | elec_coord[walk_num][3][elec_num] | in | Electron coordinates | - | double | nucl_coord[3][elec_num] | in | Nuclear coordinates | - | double | expo[prim_num] | in | Exponents of the primitives | - | double | primitive_vgl[5][walk_num][elec_num][prim_num] | out | Value, gradients and Laplacian of the primitives | + | qmckl_context | context | in | Global state | + | int64_t | prim_num | in | Number of primitives | + | int64_t | elec_num | in | Number of electrons | + | int64_t | nucl_num | in | Number of nuclei | + | int64_t | nucleus_prim_index[nucl_num] | in | Index of the 1st primitive of each nucleus | + | double | elec_coord[3][elec_num] | in | Electron coordinates | + | double | nucl_coord[3][elec_num] | in | Nuclear coordinates | + | double | expo[prim_num] | in | Exponents of the primitives | + | double | primitive_vgl[5][elec_num][prim_num] | out | Value, gradients and Laplacian of the primitives | #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_ao_basis_primitive_gaussian_vgl_f(context, & - prim_num, elec_num, nucl_num, walk_num, & + prim_num, elec_num, nucl_num, & nucleus_prim_index, elec_coord, nucl_coord, expo, primitive_vgl) & result(info) use qmckl @@ -1947,14 +1945,13 @@ integer function qmckl_compute_ao_basis_primitive_gaussian_vgl_f(context, & integer*8 , intent(in) :: prim_num integer*8 , intent(in) :: nucl_num integer*8 , intent(in) :: elec_num - integer*8 , intent(in) :: walk_num integer*8 , intent(in) :: nucleus_prim_index(nucl_num+1) - double precision , intent(in) :: elec_coord(elec_num,3,walk_num) + double precision , intent(in) :: elec_coord(elec_num,3) double precision , intent(in) :: nucl_coord(nucl_num,3) double precision , intent(in) :: expo(prim_num) - double precision , intent(out) :: primitive_vgl(prim_num,elec_num,walk_num,5) + double precision , intent(out) :: primitive_vgl(prim_num,elec_num,5) - integer*8 :: inucl, iprim, iwalk, ielec + integer*8 :: inucl, iprim, ielec double precision :: x, y, z, two_a, ar2, r2, v, cutoff info = QMCKL_SUCCESS @@ -1965,30 +1962,28 @@ integer function qmckl_compute_ao_basis_primitive_gaussian_vgl_f(context, & do inucl=1,nucl_num ! C is zero-based, so shift bounds by one do iprim = nucleus_prim_index(inucl)+1, nucleus_prim_index(inucl+1) - do iwalk = 1, walk_num - do ielec = 1, elec_num - x = elec_coord(ielec,1,iwalk) - nucl_coord(inucl,1) - y = elec_coord(ielec,2,iwalk) - nucl_coord(inucl,2) - z = elec_coord(ielec,3,iwalk) - nucl_coord(inucl,3) + do ielec = 1, elec_num + x = elec_coord(ielec,1) - nucl_coord(inucl,1) + y = elec_coord(ielec,2) - nucl_coord(inucl,2) + z = elec_coord(ielec,3) - nucl_coord(inucl,3) - r2 = x*x + y*y + z*z - ar2 = expo(iprim)*r2 - if (ar2 > cutoff) cycle - - v = dexp(-ar2) - two_a = -2.d0 * expo(iprim) * v + r2 = x*x + y*y + z*z + ar2 = expo(iprim)*r2 + if (ar2 > cutoff) cycle - primitive_vgl(iprim, ielec, iwalk, 1) = v - primitive_vgl(iprim, ielec, iwalk, 2) = two_a * x - primitive_vgl(iprim, ielec, iwalk, 3) = two_a * y - primitive_vgl(iprim, ielec, iwalk, 4) = two_a * z - primitive_vgl(iprim, ielec, iwalk, 5) = two_a * (3.d0 - 2.d0*ar2) + v = dexp(-ar2) + two_a = -2.d0 * expo(iprim) * v + + primitive_vgl(iprim, ielec, 1) = v + primitive_vgl(iprim, ielec, 2) = two_a * x + primitive_vgl(iprim, ielec, 3) = two_a * y + primitive_vgl(iprim, ielec, 4) = two_a * z + primitive_vgl(iprim, ielec, 5) = two_a * (3.d0 - 2.d0*ar2) - end do end do end do end do - + end function qmckl_compute_ao_basis_primitive_gaussian_vgl_f #+end_src @@ -1998,7 +1993,6 @@ qmckl_exit_code qmckl_compute_ao_basis_primitive_gaussian_vgl( const int64_t prim_num, const int64_t elec_num, const int64_t nucl_num, - const int64_t walk_num, const int64_t* nucleus_prim_index, const double* elec_coord, const double* nucl_coord, @@ -2011,7 +2005,15 @@ qmckl_exit_code qmckl_compute_ao_basis_primitive_gaussian_vgl( #+RESULTS: #+begin_src f90 :tangle (eval f) :comments org :exports none integer(c_int32_t) function qmckl_compute_ao_basis_primitive_gaussian_vgl & - (context, prim_num, elec_num, nucl_num, walk_num, nucleus_prim_index, elec_coord, nucl_coord, expo, primitive_vgl) & + (context, & + prim_num, & + elec_num, & + nucl_num, & + nucleus_prim_index, & + elec_coord, & + nucl_coord, & + expo, & + primitive_vgl) & bind(C) result(info) use, intrinsic :: iso_c_binding @@ -2021,16 +2023,23 @@ qmckl_exit_code qmckl_compute_ao_basis_primitive_gaussian_vgl( integer (c_int64_t) , intent(in) , value :: prim_num integer (c_int64_t) , intent(in) , value :: elec_num integer (c_int64_t) , intent(in) , value :: nucl_num - integer (c_int64_t) , intent(in) , value :: walk_num integer (c_int64_t) , intent(in) :: nucleus_prim_index(nucl_num) - real (c_double ) , intent(in) :: elec_coord(elec_num,3,walk_num) + real (c_double ) , intent(in) :: elec_coord(elec_num,3) real (c_double ) , intent(in) :: nucl_coord(elec_num,3) real (c_double ) , intent(in) :: expo(prim_num) - real (c_double ) , intent(out) :: primitive_vgl(elec_num,walk_num,5,prim_num) + real (c_double ) , intent(out) :: primitive_vgl(prim_num,elec_num,5) integer(c_int32_t), external :: qmckl_compute_ao_basis_primitive_gaussian_vgl_f info = qmckl_compute_ao_basis_primitive_gaussian_vgl_f & - (context, prim_num, elec_num, nucl_num, walk_num, nucleus_prim_index, elec_coord, nucl_coord, expo, primitive_vgl) + (context, & + prim_num, & + elec_num, & + nucl_num, & + nucleus_prim_index, & + elec_coord, & + nucl_coord, & + expo, & + primitive_vgl) end function qmckl_compute_ao_basis_primitive_gaussian_vgl #+end_src @@ -2059,40 +2068,25 @@ def lf(a,x,y): return d2f(a,x,y,1) + d2f(a,x,y,2) + d2f(a,x,y,3) elec_26_w1 = np.array( [ 1.49050402641, 2.90106987953, -1.05920815468 ] ) -elec_15_w2 = np.array( [ -2.20180344582,-1.9113150239, 2.2193744778600002 ] ) nucl_1 = np.array( [ 1.096243353458458e+00, 8.907054016973815e-01, 7.777092280258892e-01 ] ) nucl_2 = np.array( [ 1.168459237342663e+00, 1.125660720053393e+00, 2.833370314829343e+00 ] ) -#double prim_vgl[prim_num][5][walk_num][elec_num]; +#double prim_vgl[prim_num][5][elec_num]; a = 0.9059; x = elec_26_w1 ; y = nucl_1 -print ( "[7][0][0][26] : %e"% f(a,x,y)) -print ( "[7][1][0][26] : %e"% df(a,x,y,1)) -print ( "[7][2][0][26] : %e"% df(a,x,y,2)) -print ( "[7][3][0][26] : %e"% df(a,x,y,3)) -print ( "[7][4][0][26] : %e"% lf(a,x,y)) - -a = 0.32578; x = elec_15_w2 ; y = nucl_2 -print ( "[39][0][1][15] : %e"% f(a,x,y)) -print ( "[39][1][1][15] : %e"% df(a,x,y,1)) -print ( "[39][2][1][15] : %e"% df(a,x,y,2)) -print ( "[39][3][1][15] : %e"% df(a,x,y,3)) -print ( "[39][4][1][15] : %e"% lf(a,x,y)) +print ( "[7][0][26] : %e"% f(a,x,y)) +print ( "[7][1][26] : %e"% df(a,x,y,1)) +print ( "[7][2][26] : %e"% df(a,x,y,2)) +print ( "[7][3][26] : %e"% df(a,x,y,3)) +print ( "[7][4][26] : %e"% lf(a,x,y)) #+end_src #+RESULTS: - #+begin_example - [7][0][0][26] : 1.050157e-03 - [7][1][0][26] : -7.501497e-04 - [7][2][0][26] : -3.825069e-03 - [7][3][0][26] : 3.495056e-03 - [7][4][0][26] : 2.040013e-02 - [39][0][1][15] : 1.083038e-03 - [39][1][1][15] : 2.378275e-03 - [39][2][1][15] : 2.143086e-03 - [39][3][1][15] : 4.332750e-04 - [39][4][1][15] : 7.514605e-03 - #+end_example + : [7][0][26] : 1.050157e-03 + : [7][1][26] : -7.501497e-04 + : [7][2][26] : -3.825069e-03 + : [7][3][26] : 3.495056e-03 + : [7][4][26] : 2.040013e-02 *** Test @@ -2119,22 +2113,16 @@ assert(rc == QMCKL_SUCCESS); -double prim_vgl[5][walk_num][elec_num][prim_num]; +double prim_vgl[5][elec_num][prim_num]; -rc = qmckl_get_ao_basis_primitive_vgl(context, &(prim_vgl[0][0][0][0])); +rc = qmckl_get_ao_basis_primitive_vgl(context, &(prim_vgl[0][0][0])); assert (rc == QMCKL_SUCCESS); -assert( fabs(prim_vgl[0][0][26][7] - ( 1.0501570432064878E-003)) < 1.e-14 ); -assert( fabs(prim_vgl[1][0][26][7] - (-7.5014974095310560E-004)) < 1.e-14 ); -assert( fabs(prim_vgl[2][0][26][7] - (-3.8250692897610380E-003)) < 1.e-14 ); -assert( fabs(prim_vgl[3][0][26][7] - ( 3.4950559194080275E-003)) < 1.e-14 ); -assert( fabs(prim_vgl[4][0][26][7] - ( 2.0392163767356572E-002)) < 1.e-14 ); - -assert( fabs(prim_vgl[0][1][15][39] - ( 1.0825844173157661E-003)) < 1.e-14 ); -assert( fabs(prim_vgl[1][1][15][39] - ( 2.3774237611651531E-003)) < 1.e-14 ); -assert( fabs(prim_vgl[2][1][15][39] - ( 2.1423191526963063E-003)) < 1.e-14 ); -assert( fabs(prim_vgl[3][1][15][39] - ( 4.3312003523048492E-004)) < 1.e-14 ); -assert( fabs(prim_vgl[4][1][15][39] - ( 7.5174404780004771E-003)) < 1.e-14 ); +assert( fabs(prim_vgl[0][26][7] - ( 1.0501570432064878E-003)) < 1.e-14 ); +assert( fabs(prim_vgl[1][26][7] - (-7.5014974095310560E-004)) < 1.e-14 ); +assert( fabs(prim_vgl[2][26][7] - (-3.8250692897610380E-003)) < 1.e-14 ); +assert( fabs(prim_vgl[3][26][7] - ( 3.4950559194080275E-003)) < 1.e-14 ); +assert( fabs(prim_vgl[4][26][7] - ( 2.0392163767356572E-002)) < 1.e-14 ); } @@ -2144,26 +2132,23 @@ assert( fabs(prim_vgl[4][1][15][39] - ( 7.5174404780004771E-003)) < 1.e-14 ); *** Ideas for improvement #+begin_src c -// m : walkers // j : electrons // l : primitives k=0; -for (m=0 ; mao_basis.shell_num * 5 * ctx->electron.num * ctx->electron.walk_num; + size_t sze = ctx->ao_basis.shell_num * 5 * ctx->electron.num; memcpy(shell_vgl, ctx->ao_basis.shell_vgl, sze * sizeof(double)); return QMCKL_SUCCESS; @@ -2256,8 +2241,7 @@ qmckl_exit_code qmckl_provide_ao_basis_shell_vgl(qmckl_context context) if (ctx->ao_basis.shell_vgl == NULL) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; - mem_info.size = ctx->ao_basis.shell_num * 5 * ctx->electron.num * - ctx->electron.walk_num * sizeof(double); + mem_info.size = ctx->ao_basis.shell_num * 5 * ctx->electron.num * sizeof(double); double* shell_vgl = (double*) qmckl_malloc(context, mem_info); if (shell_vgl == NULL) { @@ -2276,7 +2260,6 @@ qmckl_exit_code qmckl_provide_ao_basis_shell_vgl(qmckl_context context) ctx->ao_basis.shell_num, ctx->electron.num, ctx->nucleus.num, - ctx->electron.walk_num, ctx->ao_basis.nucleus_shell_num, ctx->ao_basis.nucleus_index, ctx->ao_basis.shell_prim_index, @@ -2311,25 +2294,24 @@ qmckl_exit_code qmckl_provide_ao_basis_shell_vgl(qmckl_context context) :END: #+NAME: qmckl_ao_basis_shell_gaussian_vgl_args - | ~qmckl_context~ | ~context~ | in | Global state | - | ~int64_t~ | ~prim_num~ | in | Number of primitives | - | ~int64_t~ | ~shell_num~ | in | Number of shells | - | ~int64_t~ | ~elec_num~ | in | Number of electrons | - | ~int64_t~ | ~nucl_num~ | in | Number of nuclei | - | ~int64_t~ | ~walk_num~ | in | Number of walkers | - | ~int64_t~ | ~nucleus_shell_num[nucl_num]~ | in | Number of shells for each nucleus | - | ~int64_t~ | ~nucleus_index[nucl_num]~ | in | Index of the 1st shell of each nucleus | - | ~int64_t~ | ~shell_prim_index[shell_num]~ | in | Index of the 1st primitive of each shell | - | ~int64_t~ | ~shell_prim_num[shell_num]~ | in | Number of primitives per shell | - | ~double~ | ~elec_coord[walk_num][3][elec_num]~ | in | Electron coordinates | - | ~double~ | ~nucl_coord[3][elec_num]~ | in | Nuclear coordinates | - | ~double~ | ~expo[prim_num]~ | in | Exponents of the primitives | - | ~double~ | ~coef_normalized[prim_num]~ | in | Coefficients of the primitives | - | ~double~ | ~shell_vgl[5][walk_num][elec_num][shell_num]~ | out | Value, gradients and Laplacian of the shells | + | ~qmckl_context~ | ~context~ | in | Global state | + | ~int64_t~ | ~prim_num~ | in | Number of primitives | + | ~int64_t~ | ~shell_num~ | in | Number of shells | + | ~int64_t~ | ~elec_num~ | in | Number of electrons | + | ~int64_t~ | ~nucl_num~ | in | Number of nuclei | + | ~int64_t~ | ~nucleus_shell_num[nucl_num]~ | in | Number of shells for each nucleus | + | ~int64_t~ | ~nucleus_index[nucl_num]~ | in | Index of the 1st shell of each nucleus | + | ~int64_t~ | ~shell_prim_index[shell_num]~ | in | Index of the 1st primitive of each shell | + | ~int64_t~ | ~shell_prim_num[shell_num]~ | in | Number of primitives per shell | + | ~double~ | ~elec_coord[3][elec_num]~ | in | Electron coordinates | + | ~double~ | ~nucl_coord[3][elec_num]~ | in | Nuclear coordinates | + | ~double~ | ~expo[prim_num]~ | in | Exponents of the primitives | + | ~double~ | ~coef_normalized[prim_num]~ | in | Coefficients of the primitives | + | ~double~ | ~shell_vgl[5][elec_num][shell_num]~ | out | Value, gradients and Laplacian of the shells | #+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, elec_num, nucl_num, walk_num, & + prim_num, shell_num, elec_num, nucl_num, & nucleus_shell_num, nucleus_index, shell_prim_index, shell_prim_num, & elec_coord, nucl_coord, expo, coef_normalized, shell_vgl) & result(info) @@ -2340,79 +2322,76 @@ integer function qmckl_compute_ao_basis_shell_gaussian_vgl_f(context, & integer*8 , intent(in) :: shell_num integer*8 , intent(in) :: nucl_num integer*8 , intent(in) :: elec_num - integer*8 , intent(in) :: walk_num integer*8 , intent(in) :: nucleus_shell_num(nucl_num) integer*8 , intent(in) :: nucleus_index(nucl_num) integer*8 , intent(in) :: shell_prim_index(shell_num) integer*8 , intent(in) :: shell_prim_num(shell_num) - double precision , intent(in) :: elec_coord(elec_num,3,walk_num) + double precision , intent(in) :: elec_coord(elec_num,3) double precision , intent(in) :: nucl_coord(nucl_num,3) double precision , intent(in) :: expo(prim_num) double precision , intent(in) :: coef_normalized(prim_num) - double precision , intent(out) :: shell_vgl(shell_num,elec_num,walk_num,5) + double precision , intent(out) :: shell_vgl(shell_num,elec_num,5) - integer*8 :: inucl, iprim, iwalk, ielec, ishell + integer*8 :: inucl, iprim, ielec, ishell double precision :: x, y, z, two_a, ar2, r2, v, cutoff info = QMCKL_SUCCESS - + ! Don't compute exponentials when the result will be almost zero. ! TODO : Use numerical precision here cutoff = -dlog(1.d-15) do inucl=1,nucl_num - do iwalk = 1, walk_num - do ielec = 1, elec_num + do ielec = 1, elec_num - x = elec_coord(ielec,1,iwalk) - nucl_coord(inucl,1) - y = elec_coord(ielec,2,iwalk) - nucl_coord(inucl,2) - z = elec_coord(ielec,3,iwalk) - nucl_coord(inucl,3) + x = elec_coord(ielec,1) - nucl_coord(inucl,1) + y = elec_coord(ielec,2) - nucl_coord(inucl,2) + z = elec_coord(ielec,3) - nucl_coord(inucl,3) - r2 = x*x + y*y + z*z + r2 = x*x + y*y + z*z - do ishell=nucleus_index(inucl)+1, nucleus_index(inucl)+nucleus_shell_num(inucl) - ! C is zero-based, so shift bounds by one + do ishell=nucleus_index(inucl)+1, nucleus_index(inucl)+nucleus_shell_num(inucl) + ! C is zero-based, so shift bounds by one - shell_vgl(ishell, ielec, iwalk, 1) = 0.d0 - shell_vgl(ishell, ielec, iwalk, 2) = 0.d0 - shell_vgl(ishell, ielec, iwalk, 3) = 0.d0 - shell_vgl(ishell, ielec, iwalk, 4) = 0.d0 - shell_vgl(ishell, ielec, iwalk, 5) = 0.d0 + shell_vgl(ishell, ielec, 1) = 0.d0 + shell_vgl(ishell, ielec, 2) = 0.d0 + shell_vgl(ishell, ielec, 3) = 0.d0 + shell_vgl(ishell, ielec, 4) = 0.d0 + shell_vgl(ishell, ielec, 5) = 0.d0 - do iprim = shell_prim_index(ishell)+1, shell_prim_index(ishell)+shell_prim_num(ishell) + do iprim = shell_prim_index(ishell)+1, shell_prim_index(ishell)+shell_prim_num(ishell) - ar2 = expo(iprim)*r2 - if (ar2 > cutoff) then - cycle - end if + ar2 = expo(iprim)*r2 + if (ar2 > cutoff) then + cycle + end if - v = coef_normalized(iprim) * dexp(-ar2) - two_a = -2.d0 * expo(iprim) * v + v = coef_normalized(iprim) * dexp(-ar2) + two_a = -2.d0 * expo(iprim) * v - shell_vgl(ishell, ielec, iwalk, 1) = & - shell_vgl(ishell, ielec, iwalk, 1) + v + shell_vgl(ishell, ielec, 1) = & + shell_vgl(ishell, ielec, 1) + v - shell_vgl(ishell, ielec, iwalk, 2) = & - shell_vgl(ishell, ielec, iwalk, 2) + two_a * x + shell_vgl(ishell, ielec, 2) = & + shell_vgl(ishell, ielec, 2) + two_a * x - shell_vgl(ishell, ielec, iwalk, 3) = & - shell_vgl(ishell, ielec, iwalk, 3) + two_a * y + shell_vgl(ishell, ielec, 3) = & + shell_vgl(ishell, ielec, 3) + two_a * y - shell_vgl(ishell, ielec, iwalk, 4) = & - shell_vgl(ishell, ielec, iwalk, 4) + two_a * z + shell_vgl(ishell, ielec, 4) = & + shell_vgl(ishell, ielec, 4) + two_a * z - shell_vgl(ishell, ielec, iwalk, 5) = & - shell_vgl(ishell, ielec, iwalk, 5) + two_a * (3.d0 - 2.d0*ar2) - - end do + shell_vgl(ishell, ielec, 5) = & + shell_vgl(ishell, ielec, 5) + two_a * (3.d0 - 2.d0*ar2) end do + end do - end do + end do - + end function qmckl_compute_ao_basis_shell_gaussian_vgl_f #+end_src @@ -2426,7 +2405,6 @@ end function qmckl_compute_ao_basis_shell_gaussian_vgl_f const int64_t shell_num, const int64_t elec_num, const int64_t nucl_num, - const int64_t walk_num, const int64_t* nucleus_shell_num, const int64_t* nucleus_index, const int64_t* shell_prim_index, @@ -2448,7 +2426,6 @@ end function qmckl_compute_ao_basis_shell_gaussian_vgl_f shell_num, & elec_num, & nucl_num, & - walk_num, & nucleus_shell_num, & nucleus_index, & shell_prim_index, & @@ -2468,16 +2445,15 @@ end function qmckl_compute_ao_basis_shell_gaussian_vgl_f integer (c_int64_t) , intent(in) , value :: shell_num integer (c_int64_t) , intent(in) , value :: elec_num integer (c_int64_t) , intent(in) , value :: nucl_num - integer (c_int64_t) , intent(in) , value :: walk_num integer (c_int64_t) , intent(in) :: nucleus_shell_num(nucl_num) integer (c_int64_t) , intent(in) :: nucleus_index(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) :: elec_coord(elec_num,3,walk_num) + real (c_double ) , intent(in) :: elec_coord(elec_num,3) real (c_double ) , intent(in) :: nucl_coord(elec_num,3) real (c_double ) , intent(in) :: expo(prim_num) real (c_double ) , intent(in) :: coef_normalized(prim_num) - real (c_double ) , intent(out) :: shell_vgl(shell_num,elec_num,walk_num,5) + real (c_double ) , intent(out) :: shell_vgl(shell_num,elec_num,5) integer(c_int32_t), external :: qmckl_compute_ao_basis_shell_gaussian_vgl_f info = qmckl_compute_ao_basis_shell_gaussian_vgl_f & @@ -2486,7 +2462,6 @@ end function qmckl_compute_ao_basis_shell_gaussian_vgl_f shell_num, & elec_num, & nucl_num, & - walk_num, & nucleus_shell_num, & nucleus_index, & shell_prim_index, & @@ -2528,7 +2503,7 @@ elec_15_w2 = np.array( [ -2.20180344582,-1.9113150239, 2.2193744778600002 ] ) nucl_1 = np.array( [ 1.096243353458458e+00, 8.907054016973815e-01, 7.777092280258892e-01 ] ) nucl_2 = np.array( [ 1.168459237342663e+00, 1.125660720053393e+00, 2.833370314829343e+00 ] ) -#double prim_vgl[prim_num][5][walk_num][elec_num]; +#double prim_vgl[prim_num][5][elec_num]; x = elec_26_w1 ; y = nucl_1 a = [( 8.236000E+03, -1.130000E-04 * 6.1616545431994848e+02 ), ( 1.235000E+03, -8.780000E-04 * 1.4847738511079908e+02 ), @@ -2541,40 +2516,20 @@ a = [( 8.236000E+03, -1.130000E-04 * 6.1616545431994848e+02 ), ( 3.643000E-01, 5.986840E-01 * 3.3419848027174592e-01 ), ( 1.285000E-01, 3.953890E-01 * 1.5296336817449557e-01 )] -print ( "[1][0][0][26] : %25.15e"% f(a,x,y)) -print ( "[1][1][0][26] : %25.15e"% df(a,x,y,1)) -print ( "[1][2][0][26] : %25.15e"% df(a,x,y,2)) -print ( "[1][3][0][26] : %25.15e"% df(a,x,y,3)) -print ( "[1][4][0][26] : %25.15e"% lf(a,x,y)) - -x = elec_15_w2 ; y = nucl_2 -a = [(3.387000E+01, 6.068000E-03 *1.0006253235944540e+01), - (5.095000E+00, 4.530800E-02 *2.4169531573445120e+00), - (1.159000E+00, 2.028220E-01 *7.9610924849766440e-01), - (3.258000E-01, 5.039030E-01 *3.0734305383061117e-01), - (1.027000E-01, 3.834210E-01 *1.2929684417481876e-01)] - -print ( "[0][1][15][14] : %25.15e"% f(a,x,y)) -print ( "[1][1][15][14] : %25.15e"% df(a,x,y,1)) -print ( "[2][1][15][14] : %25.15e"% df(a,x,y,2)) -print ( "[3][1][15][14] : %25.15e"% df(a,x,y,3)) -print ( "[4][1][15][14] : %25.15e"% lf(a,x,y)) +print ( "[1][0][26] : %25.15e"% f(a,x,y)) +print ( "[1][1][26] : %25.15e"% df(a,x,y,1)) +print ( "[1][2][26] : %25.15e"% df(a,x,y,2)) +print ( "[1][3][26] : %25.15e"% df(a,x,y,3)) +print ( "[1][4][26] : %25.15e"% lf(a,x,y)) #+end_src #+RESULTS: - #+begin_example - [1][0][0][26] : 3.564393437193867e-02 - [1][1][0][26] : -6.030177988891605e-03 - [1][2][0][26] : -3.074832579871845e-02 - [1][3][0][26] : 2.809546963133958e-02 - [1][4][0][26] : 1.903338597841753e-02 - [0][1][15][14] : 5.928089771361000e-03 - [1][1][15][14] : 4.355862298893037e-03 - [2][1][15][14] : 3.925108924950765e-03 - [3][1][15][14] : 7.935527764416084e-04 - [4][1][15][14] : 2.697495005143935e-03 - #+end_example + : [1][0][26] : 3.564393437193867e-02 + : [1][1][26] : -6.030177988891605e-03 + : [1][2][26] : -3.074832579871845e-02 + : [1][3][26] : 2.809546963133958e-02 + : [1][4][26] : 1.903338597841753e-02 *** Test @@ -2600,34 +2555,23 @@ rc = qmckl_set_electron_coord (context, 'N', elec_coord); assert(rc == QMCKL_SUCCESS); -double shell_vgl[5][walk_num][elec_num][shell_num]; +double shell_vgl[5][elec_num][shell_num]; -rc = qmckl_get_ao_basis_shell_vgl(context, &(shell_vgl[0][0][0][0])); +rc = qmckl_get_ao_basis_shell_vgl(context, &(shell_vgl[0][0][0])); assert (rc == QMCKL_SUCCESS); -printf(" shell_vgl[1][0][0][26] %25.15e\n", shell_vgl[0][0][26][1]); -printf(" shell_vgl[1][1][0][26] %25.15e\n", shell_vgl[1][0][26][1]); -printf(" shell_vgl[1][2][0][26] %25.15e\n", shell_vgl[2][0][26][1]); -printf(" shell_vgl[1][3][0][26] %25.15e\n", shell_vgl[3][0][26][1]); -printf(" shell_vgl[1][4][0][26] %25.15e\n", shell_vgl[4][0][26][1]); +printf(" shell_vgl[1][0][26] %25.15e\n", shell_vgl[0][26][1]); +printf(" shell_vgl[1][1][26] %25.15e\n", shell_vgl[1][26][1]); +printf(" shell_vgl[1][2][26] %25.15e\n", shell_vgl[2][26][1]); +printf(" shell_vgl[1][3][26] %25.15e\n", shell_vgl[3][26][1]); +printf(" shell_vgl[1][4][26] %25.15e\n", shell_vgl[4][26][1]); -printf(" shell_vgl[14][0][1][15] %25.15e\n", shell_vgl[0][1][15][14]); -printf(" shell_vgl[14][1][1][15] %25.15e\n", shell_vgl[1][1][15][14]); -printf(" shell_vgl[14][2][1][15] %25.15e\n", shell_vgl[2][1][15][14]); -printf(" shell_vgl[14][3][1][15] %25.15e\n", shell_vgl[3][1][15][14]); -printf(" shell_vgl[14][4][1][15] %25.15e\n", shell_vgl[4][1][15][14]); - -assert( fabs(shell_vgl[0][0][26][1] - ( 3.564393437193868e-02)) < 1.e-14 ); -assert( fabs(shell_vgl[1][0][26][1] - (-6.030177987072189e-03)) < 1.e-14 ); -assert( fabs(shell_vgl[2][0][26][1] - (-3.074832579537582e-02)) < 1.e-14 ); -assert( fabs(shell_vgl[3][0][26][1] - ( 2.809546963519935e-02)) < 1.e-14 ); -assert( fabs(shell_vgl[4][0][26][1] - ( 1.896046117183968e-02)) < 1.e-14 ); +assert( fabs(shell_vgl[0][26][1] - ( 3.564393437193868e-02)) < 1.e-14 ); +assert( fabs(shell_vgl[1][26][1] - (-6.030177987072189e-03)) < 1.e-14 ); +assert( fabs(shell_vgl[2][26][1] - (-3.074832579537582e-02)) < 1.e-14 ); +assert( fabs(shell_vgl[3][26][1] - ( 2.809546963519935e-02)) < 1.e-14 ); +assert( fabs(shell_vgl[4][26][1] - ( 1.896046117183968e-02)) < 1.e-14 ); -assert( fabs(shell_vgl[0][1][15][14] - ( 5.928089771361000e-03)) < 1.e-14 ); -assert( fabs(shell_vgl[1][1][15][14] - ( 4.355862296021654e-03)) < 1.e-14 ); -assert( fabs(shell_vgl[2][1][15][14] - ( 3.925108924923650e-03)) < 1.e-14 ); -assert( fabs(shell_vgl[3][1][15][14] - ( 7.935527784022099e-04)) < 1.e-14 ); -assert( fabs(shell_vgl[4][1][15][14] - ( 2.708246573703548e-03)) < 1.e-14 ); } #+end_src @@ -3243,7 +3187,7 @@ qmckl_exit_code qmckl_get_ao_vgl(qmckl_context context, double* const ao_vgl) { qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - size_t sze = ctx->ao_basis.ao_num * 5 * ctx->electron.num * ctx->electron.walk_num; + size_t sze = ctx->ao_basis.ao_num * 5 * ctx->electron.num; memcpy(ao_vgl, ctx->ao_basis.ao_vgl, sze * sizeof(double)); return QMCKL_SUCCESS; @@ -3310,8 +3254,7 @@ qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context) if (ctx->ao_basis.ao_vgl == NULL) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; - mem_info.size = ctx->ao_basis.ao_num * 5 * ctx->electron.num * - ctx->electron.walk_num * sizeof(double); + mem_info.size = ctx->ao_basis.ao_num * 5 * ctx->electron.num * sizeof(double); double* ao_vgl = (double*) qmckl_malloc(context, mem_info); if (ao_vgl == NULL) { @@ -3328,7 +3271,6 @@ qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context) ctx->ao_basis.shell_num, ctx->electron.num, ctx->nucleus.num, - ctx->electron.walk_num, ctx->electron.coord_new, ctx->nucleus.coord, ctx->ao_basis.nucleus_index, @@ -3358,26 +3300,25 @@ qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context) :END: #+NAME: qmckl_ao_vgl_args - | ~qmckl_context~ | ~context~ | in | Global state | - | ~int64_t~ | ~ao_num~ | in | Number of AOs | - | ~int64_t~ | ~shell_num~ | in | Number of shells | - | ~int64_t~ | ~elec_num~ | in | Number of electrons | - | ~int64_t~ | ~nucl_num~ | in | Number of nuclei | - | ~int64_t~ | ~walk_num~ | in | Number of walkers | - | ~double~ | ~elec_coord[walk_num][3][elec_num]~ | in | Electron coordinates | - | ~double~ | ~nucl_coord[3][nucl_num]~ | in | Nuclear coordinates | - | ~int64_t~ | ~nucleus_index[nucl_num]~ | in | Index of the 1st shell of each nucleus | - | ~int64_t~ | ~nucleus_shell_num[nucl_num]~ | in | Number of shells per nucleus | - | ~double~ | ~nucleus_range[nucl_num]~ | in | Range beyond which all is zero | - | ~int32_t~ | ~nucleus_max_ang_mom[nucl_num]~ | in | Maximum angular momentum per nucleus | - | ~int32_t~ | ~shell_ang_mom[shell_num]~ | in | Angular momentum of each shell | - | ~double~ | ~ao_factor[ao_num]~ | in | Normalization factor of the AOs | - | ~double~ | ~shell_vgl[5][walk_num][elec_num][shell_num]~ | in | Value, gradients and Laplacian of the shells | - | ~double~ | ~ao_vgl[5][walk_num][elec_num][ao_num]~ | out | Value, gradients and Laplacian of the AOs | + | ~qmckl_context~ | ~context~ | in | Global state | + | ~int64_t~ | ~ao_num~ | in | Number of AOs | + | ~int64_t~ | ~shell_num~ | in | Number of shells | + | ~int64_t~ | ~elec_num~ | in | Number of electrons | + | ~int64_t~ | ~nucl_num~ | in | Number of nuclei | + | ~double~ | ~elec_coord[3][elec_num]~ | in | Electron coordinates | + | ~double~ | ~nucl_coord[3][nucl_num]~ | in | Nuclear coordinates | + | ~int64_t~ | ~nucleus_index[nucl_num]~ | in | Index of the 1st shell of each nucleus | + | ~int64_t~ | ~nucleus_shell_num[nucl_num]~ | in | Number of shells per nucleus | + | ~double~ | ~nucleus_range[nucl_num]~ | in | Range beyond which all is zero | + | ~int32_t~ | ~nucleus_max_ang_mom[nucl_num]~ | in | Maximum angular momentum per nucleus | + | ~int32_t~ | ~shell_ang_mom[shell_num]~ | in | Angular momentum of each shell | + | ~double~ | ~ao_factor[ao_num]~ | in | Normalization factor of the AOs | + | ~double~ | ~shell_vgl[5][elec_num][shell_num]~ | in | Value, gradients and Laplacian of the shells | + | ~double~ | ~ao_vgl[5][elec_num][ao_num]~ | out | Value, gradients and Laplacian of the AOs | #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_ao_vgl_f(context, & - ao_num, shell_num, elec_num, nucl_num, walk_num, & + ao_num, shell_num, elec_num, nucl_num, & elec_coord, nucl_coord, nucleus_index, nucleus_shell_num, & nucleus_range, nucleus_max_ang_mom, shell_ang_mom, & ao_factor, shell_vgl, ao_vgl) & @@ -3389,8 +3330,7 @@ integer function qmckl_compute_ao_vgl_f(context, & integer*8 , intent(in) :: shell_num integer*8 , intent(in) :: elec_num integer*8 , intent(in) :: nucl_num - integer*8 , intent(in) :: walk_num - double precision , intent(in) :: elec_coord(elec_num,3,walk_num) + double precision , intent(in) :: elec_coord(elec_num,3) double precision , intent(in) :: nucl_coord(nucl_num,3) integer*8 , intent(in) :: nucleus_index(nucl_num) integer*8 , intent(in) :: nucleus_shell_num(nucl_num) @@ -3398,13 +3338,13 @@ integer function qmckl_compute_ao_vgl_f(context, & integer , intent(in) :: nucleus_max_ang_mom(nucl_num) integer , intent(in) :: shell_ang_mom(shell_num) double precision , intent(in) :: ao_factor(ao_num) - double precision , intent(in) :: shell_vgl(shell_num,elec_num,walk_num,5) - double precision , intent(out) :: ao_vgl(ao_num,elec_num,walk_num,5) + double precision , intent(in) :: shell_vgl(shell_num,elec_num,5) + double precision , intent(out) :: ao_vgl(ao_num,elec_num,5) double precision :: e_coord(3), n_coord(3) integer*8 :: n_poly integer :: l, il, k - integer*8 :: ielec, inucl, ishell, iwalk + integer*8 :: ielec, inucl, ishell integer :: lstart(0:20) double precision :: x, y, z, r2 double precision :: cutoff @@ -3421,76 +3361,74 @@ integer function qmckl_compute_ao_vgl_f(context, & end do info = QMCKL_SUCCESS - + ! Don't compute polynomials when the radial part is zero. ! TODO : Use numerical precision here cutoff = -dlog(1.d-15) - do iwalk = 1,walk_num - do ielec = 1, elec_num - e_coord(1) = elec_coord(ielec,1,iwalk) - e_coord(2) = elec_coord(ielec,2,iwalk) - e_coord(3) = elec_coord(ielec,3,iwalk) - k=1 - do inucl=1,nucl_num - n_coord(1) = nucl_coord(inucl,1) - n_coord(2) = nucl_coord(inucl,2) - n_coord(3) = nucl_coord(inucl,3) + do ielec = 1, elec_num + e_coord(1) = elec_coord(ielec,1) + e_coord(2) = elec_coord(ielec,2) + e_coord(3) = elec_coord(ielec,3) + k=1 + do inucl=1,nucl_num + n_coord(1) = nucl_coord(inucl,1) + n_coord(2) = nucl_coord(inucl,2) + n_coord(3) = nucl_coord(inucl,3) - ! Test if the electron is in the range of the nucleus - x = e_coord(1) - n_coord(1) - y = e_coord(2) - n_coord(2) - z = e_coord(3) - n_coord(3) + ! Test if the electron is in the range of the nucleus + x = e_coord(1) - n_coord(1) + y = e_coord(2) - n_coord(2) + z = e_coord(3) - n_coord(3) - r2 = x*x + z*z + z*z + r2 = x*x + z*z + z*z - if (r2 > cutoff*nucleus_range(inucl)) then - cycle - end if + if (r2 > cutoff*nucleus_range(inucl)) then + cycle + end if - ! Compute polynomials - info = qmckl_ao_polynomial_vgl_f(context, e_coord, n_coord, & - nucleus_max_ang_mom(inucl), n_poly, powers, 3_8, & - poly_vgl, 5_8) + ! Compute polynomials + info = qmckl_ao_polynomial_vgl_f(context, e_coord, n_coord, & + nucleus_max_ang_mom(inucl), n_poly, powers, 3_8, & + poly_vgl, 5_8) - ! Loop over shells - do ishell = nucleus_index(inucl)+1, nucleus_index(inucl)+nucleus_shell_num(inucl) - l = shell_ang_mom(ishell) - do il = lstart(l), lstart(l+1)-1 - ! Value - ao_vgl(k,ielec,iwalk,1) = & - poly_vgl(1,il) * shell_vgl(ishell,ielec,iwalk,1) * ao_factor(k) + ! Loop over shells + do ishell = nucleus_index(inucl)+1, nucleus_index(inucl)+nucleus_shell_num(inucl) + l = shell_ang_mom(ishell) + do il = lstart(l), lstart(l+1)-1 + ! Value + ao_vgl(k,ielec,1) = & + poly_vgl(1,il) * shell_vgl(ishell,ielec,1) * ao_factor(k) - ! Grad_x - ao_vgl(k,ielec,iwalk,2) = ( & - poly_vgl(2,il) * shell_vgl(ishell,ielec,iwalk,1) + & - poly_vgl(1,il) * shell_vgl(ishell,ielec,iwalk,2) & - ) * ao_factor(k) + ! Grad_x + ao_vgl(k,ielec,2) = ( & + poly_vgl(2,il) * shell_vgl(ishell,ielec,1) + & + poly_vgl(1,il) * shell_vgl(ishell,ielec,2) & + ) * ao_factor(k) - ! Grad_y - ao_vgl(k,ielec,iwalk,3) = ( & - poly_vgl(3,il) * shell_vgl(ishell,ielec,iwalk,1) + & - poly_vgl(1,il) * shell_vgl(ishell,ielec,iwalk,3) & - ) * ao_factor(k) + ! Grad_y + ao_vgl(k,ielec,3) = ( & + poly_vgl(3,il) * shell_vgl(ishell,ielec,1) + & + poly_vgl(1,il) * shell_vgl(ishell,ielec,3) & + ) * ao_factor(k) - ! Grad_z - ao_vgl(k,ielec,iwalk,4) = ( & - poly_vgl(4,il) * shell_vgl(ishell,ielec,iwalk,1) + & - poly_vgl(1,il) * shell_vgl(ishell,ielec,iwalk,4) & - ) * ao_factor(k) + ! Grad_z + ao_vgl(k,ielec,4) = ( & + poly_vgl(4,il) * shell_vgl(ishell,ielec,1) + & + poly_vgl(1,il) * shell_vgl(ishell,ielec,4) & + ) * ao_factor(k) - ! Lapl_z - ao_vgl(k,ielec,iwalk,5) = ( & - poly_vgl(5,il) * shell_vgl(ishell,ielec,iwalk,1) + & - poly_vgl(1,il) * shell_vgl(ishell,ielec,iwalk,5) + & - 2.d0 * ( & - poly_vgl(2,il) * shell_vgl(ishell,ielec,iwalk,2) + & - poly_vgl(3,il) * shell_vgl(ishell,ielec,iwalk,3) + & - poly_vgl(4,il) * shell_vgl(ishell,ielec,iwalk,4) ) & - ) * ao_factor(k) - - k = k+1 - end do + ! Lapl_z + ao_vgl(k,ielec,5) = ( & + poly_vgl(5,il) * shell_vgl(ishell,ielec,1) + & + poly_vgl(1,il) * shell_vgl(ishell,ielec,5) + & + 2.d0 * ( & + poly_vgl(2,il) * shell_vgl(ishell,ielec,2) + & + poly_vgl(3,il) * shell_vgl(ishell,ielec,3) + & + poly_vgl(4,il) * shell_vgl(ishell,ielec,4) ) & + ) * ao_factor(k) + + k = k+1 end do end do end do @@ -3510,7 +3448,6 @@ end function qmckl_compute_ao_vgl_f const int64_t shell_num, const int64_t elec_num, const int64_t nucl_num, - const int64_t walk_num, const double* elec_coord, const double* nucl_coord, const int64_t* nucleus_index, @@ -3533,7 +3470,6 @@ end function qmckl_compute_ao_vgl_f shell_num, & elec_num, & nucl_num, & - walk_num, & elec_coord, & nucl_coord, & nucleus_index, & @@ -3554,8 +3490,7 @@ end function qmckl_compute_ao_vgl_f integer (c_int64_t) , intent(in) , value :: shell_num integer (c_int64_t) , intent(in) , value :: elec_num integer (c_int64_t) , intent(in) , value :: nucl_num - integer (c_int64_t) , intent(in) , value :: walk_num - real (c_double ) , intent(in) :: elec_coord(elec_num,3,walk_num) + real (c_double ) , intent(in) :: elec_coord(elec_num,3) real (c_double ) , intent(in) :: nucl_coord(nucl_num,3) integer (c_int64_t) , intent(in) :: nucleus_index(nucl_num) integer (c_int64_t) , intent(in) :: nucleus_shell_num(nucl_num) @@ -3563,8 +3498,8 @@ end function qmckl_compute_ao_vgl_f integer (c_int32_t) , intent(in) :: nucleus_max_ang_mom(nucl_num) integer (c_int32_t) , intent(in) :: shell_ang_mom(shell_num) real (c_double ) , intent(in) :: ao_factor(ao_num) - real (c_double ) , intent(in) :: shell_vgl(shell_num,elec_num,walk_num,5) - real (c_double ) , intent(out) :: ao_vgl(ao_num,elec_num,walk_num,5) + real (c_double ) , intent(in) :: shell_vgl(shell_num,elec_num,5) + real (c_double ) , intent(out) :: ao_vgl(ao_num,elec_num,5) integer(c_int32_t), external :: qmckl_compute_ao_vgl_f info = qmckl_compute_ao_vgl_f & @@ -3573,7 +3508,6 @@ end function qmckl_compute_ao_vgl_f shell_num, & elec_num, & nucl_num, & - walk_num, & elec_coord, & nucl_coord, & nucleus_index, & @@ -3590,6 +3524,7 @@ end function qmckl_compute_ao_vgl_f #+begin_src python :results output :exports none import numpy as np +from math import sqrt def f(a,x,y): return np.sum( [c * np.exp( -b*(np.linalg.norm(x-y))**2) for b,c in a] ) @@ -3615,51 +3550,52 @@ elec_26_w1 = np.array( [ 1.49050402641, 2.90106987953, -1.05920815468 ] ) elec_15_w2 = np.array( [ -2.20180344582,-1.9113150239, 2.2193744778600002 ] ) nucl_1 = np.array( [ -2.302574592081335e+00, -3.542027060505035e-01, -5.334129934317614e-02] ) -#double prim_vgl[prim_num][5][walk_num][elec_num]; +#double ao_vgl[prim_num][5][elec_num]; x = elec_26_w1 ; y = nucl_1 -a = [( 403.830000, 0.001473 * 5.9876577632594533e+04, 1.0), - ( 121.170000, 0.012672 * 7.2836806319891484e+03, 1.7320508075688774e+00), - ( 46.345000, 0.058045 * 1.3549226646722386e+03, 1.7320508075688774e+00), - ( 19.721000, 0.170510 * 3.0376315094739988e+02, 1.0), - ( 8.862400, 0.318596 * 7.4924579607137730e+01, 1.7320508075688774e+00), - ( 3.996200, 0.384502 * 1.8590543353806009e+01, 1.0), - ( 1.763600, 0.273774 * 4.4423176930919421e+00, 1.0), - ( 0.706190, 0.074397 * 8.9541051939952665e-01, 1.7320508075688774e+00)] +a = [( 403.830000, 0.001473 * 5.9876577632594533e+04), + ( 121.170000, 0.012672 * 7.2836806319891484e+03), + ( 46.345000, 0.058045 * 1.3549226646722386e+03), + ( 19.721000, 0.170510 * 3.0376315094739988e+02), + ( 8.862400, 0.318596 * 7.4924579607137730e+01), + ( 3.996200, 0.384502 * 1.8590543353806009e+01), + ( 1.763600, 0.273774 * 4.4423176930919421e+00), + ( 0.706190, 0.074397 * 8.9541051939952665e-01)] -print ( "[0][0][26][219] : %25.15e"%(f(a,x,y) * (x[0] - y[0])**2) ) -print ( "[1][0][26][219] : %25.15e"%(df(a,x,y,1)* (x[0] - y[0]) * (x[1] - y[1]) + 2.*f(a,x,y) * (x[0] - y[0])) ) +norm = sqrt(3.) +print ( "[0][26][219] : %25.15e"%(f(a,x,y) * (x[0] - y[0])**2) ) +print ( "[1][26][219] : %25.15e"%(df(a,x,y,1)* (x[0] - y[0]) * (x[1] - y[1]) + 2.*f(a,x,y) * (x[0] - y[0])) ) -print ( "[0][0][26][220] : %25.15e"%(f(a,x,y) * (x[0] - y[0]) * (x[1] - y[1])) ) -print ( "[1][0][26][220] : %25.15e"%(df(a,x,y,1)* (x[0] - y[0]) * (x[1] - y[1]) + f(a,x,y) * (x[1] - y[1])) ) +print ( "[0][26][220] : %25.15e"%(norm*f(a,x,y) * (x[0] - y[0]) * (x[1] - y[1]) )) +print ( "[1][26][220] : %25.15e"%(norm*df(a,x,y,1)* (x[0] - y[0]) * (x[1] - y[1]) + norm*f(a,x,y) * (x[1] - y[1])) ) -print ( "[0][0][26][221] : %25.15e"%(f(a,x,y) * (x[0] - y[0]) * (x[2] - y[2])) ) -print ( "[1][0][26][221] : %25.15e"%(df(a,x,y,1)* (x[0] - y[0]) * (x[2] - y[2]) + f(a,x,y) * (x[2] - y[2])) ) +print ( "[0][26][221] : %25.15e"%(norm*f(a,x,y) * (x[0] - y[0]) * (x[2] - y[2])) ) +print ( "[1][26][221] : %25.15e"%(norm*df(a,x,y,1)* (x[0] - y[0]) * (x[2] - y[2]) + norm*f(a,x,y) * (x[2] - y[2])) ) -print ( "[0][0][26][222] : %25.15e"%(f(a,x,y) * (x[1] - y[1]) * (x[1] - y[1])) ) -print ( "[1][0][26][222] : %25.15e"%(df(a,x,y,1)* (x[1] - y[1]) * (x[1] - y[1])) ) +print ( "[0][26][222] : %25.15e"%(f(a,x,y) * (x[1] - y[1]) * (x[1] - y[1])) ) +print ( "[1][26][222] : %25.15e"%(df(a,x,y,1)* (x[1] - y[1]) * (x[1] - y[1])) ) -print ( "[0][0][26][223] : %25.15e"%(f(a,x,y) * (x[1] - y[1]) * (x[2] - y[2])) ) -print ( "[1][0][26][223] : %25.15e"%(df(a,x,y,1)* (x[1] - y[1]) * (x[2] - y[2])) ) +print ( "[0][26][223] : %25.15e"%(norm*f(a,x,y) * (x[1] - y[1]) * (x[2] - y[2])) ) +print ( "[1][26][223] : %25.15e"%(norm*df(a,x,y,1)* (x[1] - y[1]) * (x[2] - y[2])) ) -print ( "[0][0][26][224] : %25.15e"%(f(a,x,y) * (x[2] - y[2]) * (x[2] - y[2])) ) -print ( "[1][0][26][224] : %25.15e"%(df(a,x,y,1)* (x[2] - y[2]) * (x[2] - y[2])) ) +print ( "[0][26][224] : %25.15e"%(f(a,x,y) * (x[2] - y[2]) * (x[2] - y[2])) ) +print ( "[1][26][224] : %25.15e"%(df(a,x,y,1)* (x[2] - y[2]) * (x[2] - y[2])) ) #+end_src #+RESULTS: #+begin_example - [0][0][26][219] : 1.020302912653649e-08 - [1][0][26][219] : -4.153046808203204e-08 - [0][0][26][220] : 8.756380857379661e-09 - [1][0][26][220] : -4.460176677299534e-08 - [0][0][26][221] : -2.705688401075445e-09 - [1][0][26][221] : 1.378177639720419e-08 - [0][0][26][222] : 7.514847283937212e-09 - [1][0][26][222] : -4.025905373647693e-08 - [0][0][26][223] : -2.322059246071533e-09 - [1][0][26][223] : 1.243989457599443e-08 - [0][0][26][224] : 7.175074806631758e-10 - [1][0][26][224] : -3.843880138733679e-09 + [0][26][219] : 1.020302912653649e-08 + [1][26][219] : -4.153046808203204e-08 + [0][26][220] : 1.516649653540510e-08 + [1][26][220] : -7.725252615816528e-08 + [0][26][221] : -4.686389780112468e-09 + [1][26][221] : 2.387073693851122e-08 + [0][26][222] : 7.514847283937212e-09 + [1][26][222] : -4.025905373647693e-08 + [0][26][223] : -4.021924592380977e-09 + [1][26][223] : 2.154652944642284e-08 + [0][26][224] : 7.175074806631758e-10 + [1][26][224] : -3.843880138733679e-09 #+end_example *** Test @@ -3687,38 +3623,38 @@ rc = qmckl_set_electron_coord (context, 'N', elec_coord); assert(rc == QMCKL_SUCCESS); -double ao_vgl[5][walk_num][elec_num][ao_num]; +double ao_vgl[5][elec_num][ao_num]; -rc = qmckl_get_ao_vgl(context, &(ao_vgl[0][0][0][0])); +rc = qmckl_get_ao_vgl(context, &(ao_vgl[0][0][0])); assert (rc == QMCKL_SUCCESS); printf("\n"); -printf(" ao_vgl ao_vgl[0][0][26][219] %25.15e\n", ao_vgl[0][0][26][219]); -printf(" ao_vgl ao_vgl[1][0][26][219] %25.15e\n", ao_vgl[1][0][26][219]); -printf(" ao_vgl ao_vgl[0][0][26][220] %25.15e\n", ao_vgl[0][0][26][220]); -printf(" ao_vgl ao_vgl[1][0][26][220] %25.15e\n", ao_vgl[1][0][26][220]); -printf(" ao_vgl ao_vgl[0][0][26][221] %25.15e\n", ao_vgl[0][0][26][221]); -printf(" ao_vgl ao_vgl[1][0][26][221] %25.15e\n", ao_vgl[1][0][26][221]); -printf(" ao_vgl ao_vgl[0][0][26][222] %25.15e\n", ao_vgl[0][0][26][222]); -printf(" ao_vgl ao_vgl[1][0][26][222] %25.15e\n", ao_vgl[1][0][26][222]); -printf(" ao_vgl ao_vgl[0][0][26][223] %25.15e\n", ao_vgl[0][0][26][223]); -printf(" ao_vgl ao_vgl[1][0][26][223] %25.15e\n", ao_vgl[1][0][26][223]); -printf(" ao_vgl ao_vgl[0][0][26][224] %25.15e\n", ao_vgl[0][0][26][224]); -printf(" ao_vgl ao_vgl[1][0][26][224] %25.15e\n", ao_vgl[1][0][26][224]); +printf(" ao_vgl ao_vgl[0][26][219] %25.15e\n", ao_vgl[0][26][219]); +printf(" ao_vgl ao_vgl[1][26][219] %25.15e\n", ao_vgl[1][26][219]); +printf(" ao_vgl ao_vgl[0][26][220] %25.15e\n", ao_vgl[0][26][220]); +printf(" ao_vgl ao_vgl[1][26][220] %25.15e\n", ao_vgl[1][26][220]); +printf(" ao_vgl ao_vgl[0][26][221] %25.15e\n", ao_vgl[0][26][221]); +printf(" ao_vgl ao_vgl[1][26][221] %25.15e\n", ao_vgl[1][26][221]); +printf(" ao_vgl ao_vgl[0][26][222] %25.15e\n", ao_vgl[0][26][222]); +printf(" ao_vgl ao_vgl[1][26][222] %25.15e\n", ao_vgl[1][26][222]); +printf(" ao_vgl ao_vgl[0][26][223] %25.15e\n", ao_vgl[0][26][223]); +printf(" ao_vgl ao_vgl[1][26][223] %25.15e\n", ao_vgl[1][26][223]); +printf(" ao_vgl ao_vgl[0][26][224] %25.15e\n", ao_vgl[0][26][224]); +printf(" ao_vgl ao_vgl[1][26][224] %25.15e\n", ao_vgl[1][26][224]); printf("\n"); -assert( fabs(ao_vgl[0][0][26][219] - ( 1.020298798341620e-08)) < 1.e-14 ); -assert( fabs(ao_vgl[1][0][26][219] - ( -4.928035238010602e-08)) < 1.e-14 ); -assert( fabs(ao_vgl[0][0][26][220] - ( 1.516643537739178e-08)) < 1.e-14 ); -assert( fabs(ao_vgl[1][0][26][220] - ( -7.725221462603871e-08)) < 1.e-14 ); -assert( fabs(ao_vgl[0][0][26][221] - ( -4.686370882518819e-09)) < 1.e-14 ); -assert( fabs(ao_vgl[1][0][26][221] - ( 2.387064067626827e-08)) < 1.e-14 ); -assert( fabs(ao_vgl[0][0][26][222] - ( 7.514816980753531e-09)) < 1.e-14 ); -assert( fabs(ao_vgl[1][0][26][222] - ( -4.025889138635182e-08)) < 1.e-14 ); -assert( fabs(ao_vgl[0][0][26][223] - ( -4.021908374204471e-09)) < 1.e-14 ); -assert( fabs(ao_vgl[1][0][26][223] - ( 2.154644255710413e-08)) < 1.e-14 ); -assert( fabs(ao_vgl[0][0][26][224] - ( 7.175045873560788e-10)) < 1.e-14 ); -assert( fabs(ao_vgl[1][0][26][224] - ( -3.843864637762753e-09)) < 1.e-14 ); +assert( fabs(ao_vgl[0][26][219] - ( 1.020298798341620e-08)) < 1.e-14 ); +assert( fabs(ao_vgl[1][26][219] - (-4.928035238010602e-08)) < 1.e-14 ); +assert( fabs(ao_vgl[0][26][220] - ( 1.516643537739178e-08)) < 1.e-14 ); +assert( fabs(ao_vgl[1][26][220] - (-7.725221462603871e-08)) < 1.e-14 ); +assert( fabs(ao_vgl[0][26][221] - (-4.686370882518819e-09)) < 1.e-14 ); +assert( fabs(ao_vgl[1][26][221] - ( 2.387064067626827e-08)) < 1.e-14 ); +assert( fabs(ao_vgl[0][26][222] - ( 7.514816980753531e-09)) < 1.e-14 ); +assert( fabs(ao_vgl[1][26][222] - (-4.025889138635182e-08)) < 1.e-14 ); +assert( fabs(ao_vgl[0][26][223] - (-4.021908374204471e-09)) < 1.e-14 ); +assert( fabs(ao_vgl[1][26][223] - ( 2.154644255710413e-08)) < 1.e-14 ); +assert( fabs(ao_vgl[0][26][224] - ( 7.175045873560788e-10)) < 1.e-14 ); +assert( fabs(ao_vgl[1][26][224] - (-3.843864637762753e-09)) < 1.e-14 ); } #+end_src