From 0c9a50a681f2d97371a9e5551c3e24f3624efff1 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 26 Jan 2022 17:06:51 +0100 Subject: [PATCH] Introduced points in electron --- org/qmckl_ao.org | 374 ++++++++++++++++++++--------------------- org/qmckl_context.org | 2 +- org/qmckl_electron.org | 322 ++++++++++++++++------------------- org/qmckl_jastrow.org | 4 +- org/qmckl_point.org | 94 +++++------ 5 files changed, 373 insertions(+), 423 deletions(-) diff --git a/org/qmckl_ao.org b/org/qmckl_ao.org index 594cfa7..8a6fbf6 100644 --- a/org/qmckl_ao.org +++ b/org/qmckl_ao.org @@ -279,11 +279,11 @@ typedef struct qmckl_ao_basis_struct { int32_t * nucleus_max_ang_mom; double * nucleus_range; double * primitive_vgl; - int64_t primitive_vgl_date; + uint64_t primitive_vgl_date; double * shell_vgl; - int64_t shell_vgl_date; + uint64_t shell_vgl_date; double * ao_vgl; - int64_t ao_vgl_date; + uint64_t ao_vgl_date; int32_t uninitialized; bool provided; @@ -2457,15 +2457,15 @@ for (int64_t i=0 ; i < ao_num ; ++i) { The following data is computed as described in the next sections: - |----------------------+----------------------------------+-----------------------------------------------------------------------------------------------| - | Variable | Type | Description | - |----------------------+----------------------------------+-----------------------------------------------------------------------------------------------| - | ~primitive_vgl~ | ~double[5][elec_num][prim_num]~ | Value, gradients, Laplacian of the primitives at electron positions | - | ~primitive_vgl_date~ | ~uint64_t~ | Last modification date of Value, gradients, Laplacian of the primitives at electron positions | - | ~shell_vgl~ | ~double[5][elec_num][shell_num]~ | Value, gradients, Laplacian of the primitives at electron positions | - | ~shell_vgl_date~ | ~uint64_t~ | Last modification date of Value, gradients, Laplacian of the AOs at electron positions | - | ~ao_vgl~ | ~double[5][elec_num][ao_num]~ | Value, gradients, Laplacian of the primitives at electron positions | - | ~ao_vgl_date~ | ~uint64_t~ | Last modification date of Value, gradients, Laplacian of the AOs at electron positions | + |----------------------+-----------------------------------+----------------------------------------------------------------------------------------------| + | Variable | Type | Description | + |----------------------+-----------------------------------+----------------------------------------------------------------------------------------------| + | ~primitive_vgl~ | ~double[5][point_num][prim_num]~ | Value, gradients, Laplacian of the primitives at current positions | + | ~primitive_vgl_date~ | ~uint64_t~ | Last modification date of Value, gradients, Laplacian of the primitives at current positions | + | ~shell_vgl~ | ~double[5][point_num][shell_num]~ | Value, gradients, Laplacian of the primitives at current positions | + | ~shell_vgl_date~ | ~uint64_t~ | Last modification date of Value, gradients, Laplacian of the AOs at current positions | + | ~ao_vgl~ | ~double[5][point_num][ao_num]~ | Value, gradients, Laplacian of the primitives at current positions | + | ~ao_vgl_date~ | ~uint64_t~ | Last modification date of Value, gradients, Laplacian of the AOs at current positions | *** After initialization @@ -2653,7 +2653,7 @@ qmckl_get_ao_basis_primitive_vgl (qmckl_context context, qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - int64_t sze = ctx->ao_basis.prim_num * 5 * ctx->electron.num; + int64_t sze = ctx->ao_basis.prim_num * 5 * ctx->point.num; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, @@ -2714,7 +2714,7 @@ qmckl_get_ao_basis_shell_vgl (qmckl_context context, qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - int64_t sze = ctx->ao_basis.shell_num * 5 * ctx->electron.num; + int64_t sze = ctx->ao_basis.shell_num * 5 * ctx->point.num; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, @@ -2777,7 +2777,7 @@ qmckl_get_ao_basis_ao_vgl (qmckl_context context, qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - int64_t sze = ctx->ao_basis.ao_num * 5 * ctx->electron.num; + int64_t sze = ctx->ao_basis.ao_num * 5 * ctx->point.num; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, @@ -3045,17 +3045,17 @@ assert(0 == test_qmckl_ao_gaussian_vgl(context)); :END: #+NAME: qmckl_ao_basis_primitive_gaussian_vgl_args - | Variable | Type | In/Out | Description | - |----------------------+---------------------------------+--------+--------------------------------------------------| - | ~context~ | ~qmckl_context~ | in | Global state | - | ~prim_num~ | ~int64_t~ | in | Number of primitives | - | ~elec_num~ | ~int64_t~ | in | Number of electrons | - | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | - | ~nucleus_prim_index~ | ~int64_t[nucl_num]~ | in | Index of the 1st primitive of each nucleus | - | ~elec_coord~ | ~double[3][elec_num]~ | in | Electron coordinates | - | ~nucl_coord~ | ~double[3][elec_num]~ | in | Nuclear coordinates | - | ~expo~ | ~double[prim_num]~ | in | Exponents of the primitives | - | ~primitive_vgl~ | ~double[5][elec_num][prim_num]~ | out | Value, gradients and Laplacian of the primitives | + | Variable | Type | In/Out | Description | + |----------------------+----------------------------------+--------+--------------------------------------------------| + | ~context~ | ~qmckl_context~ | in | Global state | + | ~prim_num~ | ~int64_t~ | in | Number of primitives | + | ~point_num~ | ~int64_t~ | in | Number of points considered | + | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | + | ~nucleus_prim_index~ | ~int64_t[nucl_num]~ | in | Index of the 1st primitive of each nucleus | + | ~coord~ | ~double[3][point_num]~ | in | Coordinates | + | ~nucl_coord~ | ~double[3][nucl_num]~ | in | Nuclear coordinates | + | ~expo~ | ~double[prim_num]~ | in | Exponents of the primitives | + | ~primitive_vgl~ | ~double[5][point_num][prim_num]~ | out | Value, gradients and Laplacian of the primitives | #+CALL: generate_c_header(table=qmckl_ao_basis_primitive_gaussian_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_ao_basis_primitive_gaussian_vgl")) @@ -3064,20 +3064,20 @@ assert(0 == test_qmckl_ao_gaussian_vgl(context)); qmckl_exit_code qmckl_compute_ao_basis_primitive_gaussian_vgl ( const qmckl_context context, const int64_t prim_num, - const int64_t elec_num, + const int64_t point_num, const int64_t nucl_num, const int64_t* nucleus_prim_index, - const double* elec_coord, + const double* coord, const double* nucl_coord, const double* expo, - double* const primitive_vgl ); + double* const primitive_vgl ); #+end_src #+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, & - nucleus_prim_index, elec_coord, nucl_coord, & + context, prim_num, point_num, nucl_num, & + nucleus_prim_index, coord, nucl_coord, & expo, primitive_vgl) & result(info) @@ -3086,14 +3086,14 @@ integer function qmckl_compute_ao_basis_primitive_gaussian_vgl_f( & integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: prim_num integer*8 , intent(in) :: nucl_num - integer*8 , intent(in) :: elec_num + integer*8 , intent(in) :: point_num integer*8 , intent(in) :: nucleus_prim_index(nucl_num+1) - double precision , intent(in) :: elec_coord(elec_num,3) + double precision , intent(in) :: coord(point_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,5) + double precision , intent(out) :: primitive_vgl(prim_num,point_num,5) - integer*8 :: inucl, iprim, ielec + integer*8 :: inucl, iprim, ipoint double precision :: x, y, z, two_a, ar2, r2, v, cutoff info = QMCKL_SUCCESS @@ -3104,10 +3104,10 @@ integer function qmckl_compute_ao_basis_primitive_gaussian_vgl_f( & 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 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) + do ipoint = 1, point_num + x = coord(ipoint,1) - nucl_coord(inucl,1) + y = coord(ipoint,2) - nucl_coord(inucl,2) + z = coord(ipoint,3) - nucl_coord(inucl,3) r2 = x*x + y*y + z*z ar2 = expo(iprim)*r2 @@ -3116,11 +3116,11 @@ integer function qmckl_compute_ao_basis_primitive_gaussian_vgl_f( & 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) + primitive_vgl(iprim, ipoint, 1) = v + primitive_vgl(iprim, ipoint, 2) = two_a * x + primitive_vgl(iprim, ipoint, 3) = two_a * y + primitive_vgl(iprim, ipoint, 4) = two_a * z + primitive_vgl(iprim, ipoint, 5) = two_a * (3.d0 - 2.d0*ar2) end do end do @@ -3136,10 +3136,10 @@ end function qmckl_compute_ao_basis_primitive_gaussian_vgl_f integer(c_int32_t) function qmckl_compute_ao_basis_primitive_gaussian_vgl & (context, & prim_num, & - elec_num, & + point_num, & nucl_num, & nucleus_prim_index, & - elec_coord, & + coord, & nucl_coord, & expo, & primitive_vgl) & @@ -3150,22 +3150,22 @@ end function qmckl_compute_ao_basis_primitive_gaussian_vgl_f integer (c_int64_t) , intent(in) , value :: context 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 :: point_num integer (c_int64_t) , intent(in) , value :: nucl_num integer (c_int64_t) , intent(in) :: nucleus_prim_index(nucl_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) :: coord(point_num,3) + real (c_double ) , intent(in) :: nucl_coord(nucl_num,3) real (c_double ) , intent(in) :: expo(prim_num) - real (c_double ) , intent(out) :: primitive_vgl(prim_num,elec_num,5) + real (c_double ) , intent(out) :: primitive_vgl(prim_num,point_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, & + point_num, & nucl_num, & nucleus_prim_index, & - elec_coord, & + coord, & nucl_coord, & expo, & primitive_vgl) @@ -3201,13 +3201,13 @@ qmckl_exit_code qmckl_provide_ao_basis_primitive_vgl(qmckl_context context) } /* Compute if necessary */ - if (ctx->electron.coord_new_date > ctx->ao_basis.primitive_vgl_date) { + if (ctx->point.date > ctx->ao_basis.primitive_vgl_date) { /* Allocate array */ if (ctx->ao_basis.primitive_vgl == NULL) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; - mem_info.size = ctx->ao_basis.prim_num * 5 * ctx->electron.num * + mem_info.size = ctx->ao_basis.prim_num * 5 * ctx->point.num * sizeof(double); double* primitive_vgl = (double*) qmckl_malloc(context, mem_info); @@ -3224,10 +3224,10 @@ qmckl_exit_code qmckl_provide_ao_basis_primitive_vgl(qmckl_context context) if (ctx->ao_basis.type == 'G') { rc = qmckl_compute_ao_basis_primitive_gaussian_vgl(context, ctx->ao_basis.prim_num, - ctx->electron.num, + ctx->point.num, ctx->nucleus.num, ctx->ao_basis.nucleus_prim_index, - ctx->electron.coord_new, + ctx->point.coord.data, ctx->nucleus.coord.data, ctx->ao_basis.exponent, ctx->ao_basis.primitive_vgl); @@ -3312,15 +3312,15 @@ print ( "[7][4][26] : %e"% lf(a,x,y)) assert(qmckl_electron_provided(context)); - rc = qmckl_set_electron_coord (context, 'N', elec_coord, walk_num*3*elec_num); + rc = qmckl_set_electron_coord (context, 'N', elec_coord, walk_num*elec_num*3); assert(rc == QMCKL_SUCCESS); - double prim_vgl[5][elec_num][prim_num]; + double prim_vgl[5][elec_num*walk_num][prim_num]; rc = qmckl_get_ao_basis_primitive_vgl(context, &(prim_vgl[0][0][0]), - (int64_t) 5*elec_num*prim_num ); + (int64_t) 5*elec_num*walk_num*prim_num ); assert (rc == QMCKL_SUCCESS); assert( fabs(prim_vgl[0][26][7] - ( 1.0501570432064878E-003)) < 1.e-14 ); @@ -3341,10 +3341,10 @@ print ( "[7][4][26] : %e"% lf(a,x,y)) // l : primitives k=0; -for (j=0 ; jelectron.provided)) { - return qmckl_failwith( context, - QMCKL_NOT_PROVIDED, - "qmckl_electron", - NULL); - } - /* Compute if necessary */ - if (ctx->electron.coord_new_date > ctx->ao_basis.shell_vgl_date) { + if (ctx->point.date > ctx->ao_basis.shell_vgl_date) { /* Allocate array */ 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 * sizeof(double); + mem_info.size = ctx->ao_basis.shell_num * 5 * ctx->point.num * sizeof(double); double* shell_vgl = (double*) qmckl_malloc(context, mem_info); if (shell_vgl == NULL) { @@ -3622,13 +3615,13 @@ qmckl_exit_code qmckl_provide_ao_basis_shell_vgl(qmckl_context context) rc = qmckl_compute_ao_basis_shell_gaussian_vgl(context, ctx->ao_basis.prim_num, ctx->ao_basis.shell_num, - ctx->electron.num, + ctx->point.num, ctx->nucleus.num, ctx->ao_basis.nucleus_shell_num, ctx->ao_basis.nucleus_index, ctx->ao_basis.shell_prim_index, ctx->ao_basis.shell_prim_num, - ctx->electron.coord_new, + ctx->point.coord.data, ctx->nucleus.coord.data, ctx->ao_basis.exponent, ctx->ao_basis.coefficient_normalized, @@ -3680,7 +3673,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][elec_num]; +#double prim_vgl[prim_num][5][point_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 ), @@ -3726,7 +3719,7 @@ print ( "[1][4][26] : %25.15e"% lf(a,x,y)) assert(qmckl_electron_provided(context)); - rc = qmckl_set_electron_coord (context, 'N', elec_coord, walk_num*3*elec_num); + rc = qmckl_set_electron_coord (context, 'N', elec_coord, walk_num*elec_num*3); assert(rc == QMCKL_SUCCESS); @@ -4340,28 +4333,28 @@ end function test_qmckl_ao_polynomial_vgl :END: #+NAME: qmckl_ao_vgl_args - | 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 | - | ~elec_num~ | ~int64_t~ | in | Number of electrons | - | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | - | ~elec_coord~ | ~double[3][elec_num]~ | in | Electron 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 | - | ~ao_factor~ | ~double[ao_num]~ | in | Normalization factor of the AOs | - | ~shell_vgl~ | ~double[5][elec_num][shell_num]~ | in | Value, gradients and Laplacian of the shells | - | ~ao_vgl~ | ~double[5][elec_num][ao_num]~ | out | Value, gradients and Laplacian of the AOs | + | 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 | + | ~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 | + | ~ao_factor~ | ~double[ao_num]~ | in | Normalization factor of the AOs | + | ~shell_vgl~ | ~double[5][point_num][shell_num]~ | in | Value, gradients and Laplacian of the shells | + | ~ao_vgl~ | ~double[5][point_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, & - elec_coord, nucl_coord, nucleus_index, nucleus_shell_num, & + ao_num, shell_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) & result(info) @@ -4370,9 +4363,9 @@ integer function qmckl_compute_ao_vgl_f(context, & integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: ao_num integer*8 , intent(in) :: shell_num - integer*8 , intent(in) :: elec_num + integer*8 , intent(in) :: point_num integer*8 , intent(in) :: nucl_num - double precision , intent(in) :: elec_coord(elec_num,3) + double precision , intent(in) :: coord(point_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) @@ -4380,13 +4373,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,5) - double precision , intent(out) :: ao_vgl(ao_num,elec_num,5) + double precision , intent(in) :: shell_vgl(shell_num,point_num,5) + double precision , intent(out) :: ao_vgl(ao_num,point_num,5) double precision :: e_coord(3), n_coord(3) integer*8 :: n_poly integer :: l, il, k - integer*8 :: ielec, inucl, ishell + integer*8 :: ipoint, inucl, ishell integer*8 :: ishell_start, ishell_end integer :: lstart(0:20) double precision :: x, y, z, r2 @@ -4409,17 +4402,17 @@ integer function qmckl_compute_ao_vgl_f(context, & ! TODO : Use numerical precision here cutoff = -dlog(1.d-15) - 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) + do ipoint = 1, point_num + 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) n_coord(3) = nucl_coord(inucl,3) - ! Test if the electron is in the range of the nucleus + ! Test if the point 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) @@ -4442,35 +4435,35 @@ integer function qmckl_compute_ao_vgl_f(context, & 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) + ao_vgl(k,ipoint,1) = & + poly_vgl(1,il) * shell_vgl(ishell,ipoint,1) * 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_vgl(k,ipoint,2) = ( & + poly_vgl(2,il) * shell_vgl(ishell,ipoint,1) + & + poly_vgl(1,il) * shell_vgl(ishell,ipoint,2) & ) * 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_vgl(k,ipoint,3) = ( & + poly_vgl(3,il) * shell_vgl(ishell,ipoint,1) + & + poly_vgl(1,il) * shell_vgl(ishell,ipoint,3) & ) * 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_vgl(k,ipoint,4) = ( & + poly_vgl(4,il) * shell_vgl(ishell,ipoint,1) + & + poly_vgl(1,il) * shell_vgl(ishell,ipoint,4) & ) * ao_factor(k) ! 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) + & + ao_vgl(k,ipoint,5) = ( & + poly_vgl(5,il) * shell_vgl(ishell,ipoint,1) + & + poly_vgl(1,il) * shell_vgl(ishell,ipoint,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) ) & + poly_vgl(2,il) * shell_vgl(ishell,ipoint,2) + & + poly_vgl(3,il) * shell_vgl(ishell,ipoint,3) + & + poly_vgl(4,il) * shell_vgl(ishell,ipoint,4) ) & ) * ao_factor(k) k = k+1 @@ -4492,9 +4485,9 @@ end function qmckl_compute_ao_vgl_f const qmckl_context context, const int64_t ao_num, const int64_t shell_num, - const int64_t elec_num, + const int64_t point_num, const int64_t nucl_num, - const double* elec_coord, + const double* coord, const double* nucl_coord, const int64_t* nucleus_index, const int64_t* nucleus_shell_num, @@ -4514,9 +4507,9 @@ end function qmckl_compute_ao_vgl_f (context, & ao_num, & shell_num, & - elec_num, & + point_num, & nucl_num, & - elec_coord, & + coord, & nucl_coord, & nucleus_index, & nucleus_shell_num, & @@ -4534,9 +4527,9 @@ end function qmckl_compute_ao_vgl_f 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 :: elec_num + integer (c_int64_t) , intent(in) , value :: point_num integer (c_int64_t) , intent(in) , value :: nucl_num - real (c_double ) , intent(in) :: elec_coord(elec_num,3) + 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) @@ -4544,17 +4537,17 @@ 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,5) - real (c_double ) , intent(out) :: ao_vgl(ao_num,elec_num,5) + real (c_double ) , intent(in) :: shell_vgl(shell_num,point_num,5) + real (c_double ) , intent(out) :: ao_vgl(ao_num,point_num,5) integer(c_int32_t), external :: qmckl_compute_ao_vgl_f info = qmckl_compute_ao_vgl_f & (context, & ao_num, & shell_num, & - elec_num, & + point_num, & nucl_num, & - elec_coord, & + coord, & nucl_coord, & nucleus_index, & nucleus_shell_num, & @@ -4595,15 +4588,8 @@ qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context) NULL); } - if(!(ctx->electron.provided)) { - return qmckl_failwith( context, - QMCKL_NOT_PROVIDED, - "qmckl_electron", - NULL); - } - /* Compute if necessary */ - if (ctx->electron.coord_new_date > ctx->ao_basis.ao_vgl_date) { + if (ctx->point.date > ctx->ao_basis.ao_vgl_date) { qmckl_exit_code rc; @@ -4617,7 +4603,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 * sizeof(double); + mem_info.size = ctx->ao_basis.ao_num * 5 * ctx->point.num * sizeof(double); double* ao_vgl = (double*) qmckl_malloc(context, mem_info); if (ao_vgl == NULL) { @@ -4632,9 +4618,9 @@ qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context) rc = qmckl_compute_ao_vgl(context, ctx->ao_basis.ao_num, ctx->ao_basis.shell_num, - ctx->electron.num, + ctx->point.num, ctx->nucleus.num, - ctx->electron.coord_new, + ctx->point.coord.data, ctx->nucleus.coord.data, ctx->ao_basis.nucleus_index, ctx->ao_basis.nucleus_shell_num, @@ -4752,7 +4738,7 @@ assert (rc == QMCKL_SUCCESS); assert(qmckl_electron_provided(context)); -rc = qmckl_set_electron_coord (context, 'N', elec_coord, walk_num*3*elec_num); +rc = qmckl_set_electron_coord (context, 'N', elec_coord, walk_num*elec_num*3); assert(rc == QMCKL_SUCCESS); diff --git a/org/qmckl_context.org b/org/qmckl_context.org index 5dadf94..9c54c25 100644 --- a/org/qmckl_context.org +++ b/org/qmckl_context.org @@ -124,7 +124,7 @@ typedef struct qmckl_context_struct { uint64_t date; /* Points */ - qmckl_point_struct *point; + qmckl_point_struct point; /* -- Molecular system -- */ qmckl_nucleus_struct nucleus; diff --git a/org/qmckl_electron.org b/org/qmckl_electron.org index c09a9ca..9f79bcc 100644 --- a/org/qmckl_electron.org +++ b/org/qmckl_electron.org @@ -69,22 +69,22 @@ int main() { The following data stored in the context: - | Variable | Type | Description | - |---------------------------+----------------------------+-------------------------------------------| - | ~uninitialized~ | ~int32_t~ | Keeps bit set for uninitialized data | - | ~num~ | ~int64_t~ | Total number of electrons | - | ~up_num~ | ~int64_t~ | Number of up-spin electrons | - | ~down_num~ | ~int64_t~ | Number of down-spin electrons | - | ~walk_num~ | ~int64_t~ | Number of walkers | - | ~rescale_factor_kappa_ee~ | ~double~ | The distance scaling factor | - | ~rescale_factor_kappa_en~ | ~double~ | The distance scaling factor | - | ~provided~ | ~bool~ | If true, ~electron~ is valid | - | ~coord_new~ | ~double[3][walk_num][num]~ | New set of electron coordinates | - | ~coord_old~ | ~double[3][walk_num][num]~ | Old set of electron coordinates | - | ~coord_new_date~ | ~uint64_t~ | Last modification date of the coordinates | + | Variable | Type | Description | + |---------------------------+----------------+---------------------------------------------------------------| + | ~uninitialized~ | ~int32_t~ | Keeps bit set for uninitialized data | + | ~num~ | ~int64_t~ | Total number of electrons | + | ~up_num~ | ~int64_t~ | Number of up-spin electrons | + | ~down_num~ | ~int64_t~ | Number of down-spin electrons | + | ~walk_num~ | ~int64_t~ | Number of walkers | + | ~rescale_factor_kappa_ee~ | ~double~ | The distance scaling factor | + | ~rescale_factor_kappa_en~ | ~double~ | The distance scaling factor | + | ~provided~ | ~bool~ | If true, ~electron~ is valid | + | ~coord_new~ | ~qmckl_matrix~ | Current set of electron coordinates. Pointer to ~ctx->points~ | + | ~coord_old~ | ~qmckl_matrix~ | Old set of electron coordinates | + | ~coord_new_date~ | ~uint64_t~ | Last modification date of the coordinates | Computed data: - + | Variable | Type | Description | |-------------------------------------+--------------------------------------+----------------------------------------------------------------------| | ~ee_distance~ | ~double[walk_num][num][num]~ | Electron-electron distances | @@ -108,33 +108,33 @@ int main() { #+begin_src c :comments org :tangle (eval h_private_type) typedef struct qmckl_electron_struct { - int64_t num; - int64_t up_num; - int64_t down_num; - int64_t walk_num; - double rescale_factor_kappa_ee; - double rescale_factor_kappa_en; - int64_t coord_new_date; - int64_t ee_distance_date; - int64_t en_distance_date; - int64_t ee_pot_date; - int64_t en_pot_date; - int64_t ee_distance_rescaled_date; - int64_t ee_distance_rescaled_deriv_e_date; - int64_t en_distance_rescaled_date; - int64_t en_distance_rescaled_deriv_e_date; - double* coord_new; - double* coord_old; - double* ee_distance; - double* en_distance; - double* ee_pot; - double* en_pot; - double* ee_distance_rescaled; - double* ee_distance_rescaled_deriv_e; - double* en_distance_rescaled; - double* en_distance_rescaled_deriv_e; - int32_t uninitialized; - bool provided; + int64_t num; + int64_t up_num; + int64_t down_num; + int64_t walk_num; + double rescale_factor_kappa_ee; + double rescale_factor_kappa_en; + int64_t coord_new_date; + int64_t ee_distance_date; + int64_t en_distance_date; + int64_t ee_pot_date; + int64_t en_pot_date; + int64_t ee_distance_rescaled_date; + int64_t ee_distance_rescaled_deriv_e_date; + int64_t en_distance_rescaled_date; + int64_t en_distance_rescaled_deriv_e_date; + qmckl_matrix coord_new; + qmckl_matrix coord_old; + double* ee_distance; + double* en_distance; + double* ee_pot; + double* en_pot; + double* ee_distance_rescaled; + double* ee_distance_rescaled_deriv_e; + double* en_distance_rescaled; + double* en_distance_rescaled_deriv_e; + int32_t uninitialized; + bool provided; } qmckl_electron_struct; #+end_src @@ -146,10 +146,10 @@ typedef struct qmckl_electron_struct { Some values are initialized by default, and are not concerned by this mechanism. - #+begin_src c :comments org :tangle (eval h_private_func) + #+begin_src c :comments org :tangle (eval h_private_func) qmckl_exit_code qmckl_init_electron(qmckl_context context); #+end_src - + #+begin_src c :comments org :tangle (eval c) qmckl_exit_code qmckl_init_electron(qmckl_context context) { @@ -169,8 +169,8 @@ qmckl_exit_code qmckl_init_electron(qmckl_context context) { return QMCKL_SUCCESS; } #+end_src - - + + #+begin_src c :comments org :tangle (eval h_func) bool qmckl_electron_provided (const qmckl_context context); #+end_src @@ -204,7 +204,7 @@ if ( (ctx->electron.uninitialized & mask) != 0) { return NULL; } #+end_src - + *** Number of electrons #+begin_src c :comments org :tangle (eval h_func) :exports none @@ -395,7 +395,7 @@ qmckl_get_electron_rescale_factor_en (const qmckl_context context, double* const *** Electron coordinates Returns the current electron coordinates. The pointer is assumed - to point on a memory block of size ~size_max~ \ge ~3 * elec_num * walk_num~. + to point on a memory block of size ~size_max~ \ge ~3 * elec_num * walk_num~. The order of the indices is: | | Normal | Transposed | @@ -412,6 +412,10 @@ qmckl_get_electron_coord (const qmckl_context context, const int64_t size_max); #+end_src + As the ~coord_new~ attribute is a pointer equal to ~points~, + returning the current electron coordinates is equivalent to + returning the current points. + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_electron_coord (const qmckl_context context, @@ -419,10 +423,6 @@ qmckl_get_electron_coord (const qmckl_context context, double* const coord, const int64_t size_max) { - if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { - return QMCKL_INVALID_CONTEXT; - } - if (transp != 'N' && transp != 'T') { return qmckl_failwith( context, QMCKL_INVALID_ARG_2, @@ -437,46 +437,32 @@ qmckl_get_electron_coord (const qmckl_context context, "coord is a null pointer"); } - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; - assert (ctx != NULL); - - if ( !(ctx->electron.provided) ) { - return qmckl_failwith( context, - QMCKL_NOT_PROVIDED, - "qmckl_get_electron_coord", - "electron data is not provided"); - } - - int64_t elec_num = ctx->electron.num; - int64_t walk_num = ctx->electron.walk_num; - - if (size_max < elec_num*walk_num*3) { + if (size_max <= 0) { return qmckl_failwith( context, QMCKL_INVALID_ARG_4, "qmckl_get_electron_coord", - "Array too small"); + "size_max should be > 0"); } - assert (ctx->electron.coord_new != NULL); - - double* ptr1 = ctx->electron.coord_new; - double* ptr2 = coord; - - if (transp == 'N') { - - for (int64_t i=0 ; ielectron.provided) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_get_electron_coord", + NULL); + } + + assert (ctx->point.coord.data == ctx->electron.coord_new.data); + assert (ctx->point.coord.size[0] == ctx->electron.coord_new.size[0]); + assert (ctx->point.coord.size[1] == ctx->electron.coord_new.size[1]); + + return qmckl_get_point(context, transp, coord, size_max); } #+end_src @@ -512,20 +498,23 @@ ctx->electron.uninitialized &= ~mask; ctx->electron.provided = (ctx->electron.uninitialized == 0); if (ctx->electron.provided) { - if (ctx->electron.coord_new != NULL) { - qmckl_free(context, ctx->electron.coord_new); - ctx->electron.coord_new = NULL; + if (ctx->electron.coord_new.data != NULL) { + const qmckl_exit_code rc = qmckl_matrix_free(context, ctx->electron.coord_new); + assert (rc == QMCKL_SUCCESS); } - if (ctx->electron.coord_old != NULL) { - qmckl_free(context, ctx->electron.coord_old); - ctx->electron.coord_old = NULL; + if (ctx->electron.coord_old.data != NULL) { + const qmckl_exit_code rc = qmckl_matrix_free(context, ctx->electron.coord_old); + assert (rc == QMCKL_SUCCESS); } - qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; - mem_info.size = ctx->electron.num * ctx->electron.walk_num * 3 * sizeof(double); + const int64_t walk_num = ctx->electron.walk_num; + const int64_t elec_num = ctx->electron.num; - double* coord_new = (double*) qmckl_malloc(context, mem_info); - if (coord_new == NULL) { + const int64_t point_num = walk_num * elec_num; + + qmckl_matrix coord_new = qmckl_matrix_alloc(context, point_num, 3); + + if (coord_new.data == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_set_electron_num", @@ -533,8 +522,9 @@ if (ctx->electron.provided) { } ctx->electron.coord_new = coord_new; - double* coord_old = (double*) qmckl_malloc(context, mem_info); - if (coord_old == NULL) { + qmckl_matrix coord_old = qmckl_matrix_alloc(context, point_num, 3); + + if (coord_old.data == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_set_electron_num", @@ -542,6 +532,8 @@ if (ctx->electron.provided) { } ctx->electron.coord_old = coord_old; + ctx->point.num = point_num; + ctx->point.coord = coord_new; } return QMCKL_SUCCESS; @@ -648,9 +640,9 @@ interface import implicit none - integer (c_int64_t) , intent(in) , value :: context + integer (c_int64_t) , intent(in) , value :: context integer (c_int64_t) , intent(in) , value :: alpha - integer (c_int64_t) , intent(in) , value :: beta + integer (c_int64_t) , intent(in) , value :: beta end function end interface @@ -660,7 +652,7 @@ interface import implicit none - integer (c_int64_t) , intent(in) , value :: context + integer (c_int64_t) , intent(in) , value :: context integer (c_int64_t) , intent(in) , value :: walk_num end function end interface @@ -672,6 +664,9 @@ end interface overwritten. This can be done only when the data relative to electrons have been set. + ~size_max~ should be equal to ~elec_num * walk_num * 3~, to be symmetric + with ~qmckl_get_electron_coord~. + Important: changing the electron coordinates increments the date in the context. @@ -699,10 +694,7 @@ qmckl_set_electron_coord(qmckl_context context, "coord is a null pointer"); } - int64_t elec_num; - qmckl_exit_code rc; - rc = qmckl_get_electron_num(context, &elec_num); - if (rc != QMCKL_SUCCESS) return rc; + const int64_t elec_num = ctx->electron.num; if (elec_num == 0L) { return qmckl_failwith( context, @@ -711,9 +703,9 @@ qmckl_set_electron_coord(qmckl_context context, "elec_num is not set"); } - int64_t walk_num; - rc = qmckl_get_electron_walk_num(context, &walk_num); - if (rc != QMCKL_SUCCESS) return rc; + const int64_t walk_num = ctx->electron.walk_num; + + const int64_t point_num = ctx->point.num ; if (walk_num == 0L) { return qmckl_failwith( context, @@ -722,40 +714,18 @@ qmckl_set_electron_coord(qmckl_context context, "walk_num is not set"); } - if (size_max < walk_num*3*elec_num) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_4, - "qmckl_set_electron_coord", - "Array too small"); - } - - /* If num and walk_num are set, the arrays should be allocated */ - assert (ctx->electron.coord_old != NULL); - assert (ctx->electron.coord_new != NULL); - - /* Increment the date of the context */ - ctx->date += 1UL; + assert(point_num == elec_num * walk_num); /* Swap pointers */ - double * swap; - swap = ctx->electron.coord_old; - ctx->electron.coord_old = ctx->electron.coord_new; - ctx->electron.coord_new = swap; + ctx->point.coord = ctx->electron.coord_old; + ctx->electron.coord_old = ctx->electron.coord_new ; - double* ptr1 = ctx->electron.coord_new; - if (transp == 'N') { + qmckl_exit_code rc; + rc = qmckl_set_point(context, transp, coord, size_max/3); + assert (rc == QMCKL_SUCCESS); - for (int64_t i=0 ; ielectron.coord_new = ctx->point.coord ; - memcpy(ptr1, coord, 3*elec_num*walk_num*sizeof(double)); - - } ctx->electron.coord_new_date = ctx->date; return QMCKL_SUCCESS; @@ -770,8 +740,8 @@ interface import implicit none - integer (c_int64_t) , intent(in) , value :: context - character , intent(in) , value :: transp + integer (c_int64_t) , intent(in) , value :: context + character , intent(in) , value :: transp double precision , intent(in) :: coord(*) integer (c_int64_t) , intent(in) , value :: size_max end function @@ -880,6 +850,7 @@ double elec_coord2[walk_num*3*elec_num]; rc = qmckl_get_electron_coord (context, 'N', elec_coord2, walk_num*3*elec_num); assert(rc == QMCKL_SUCCESS); for (int64_t i=0 ; i<3*elec_num*walk_num ; ++i) { + printf("%f %f\n", elec_coord[i], elec_coord2[i]); assert( elec_coord[i] == elec_coord2[i] ); } @@ -894,7 +865,7 @@ for (int64_t i=0 ; i<3*elec_num*walk_num ; ++i) { the dependencies are more recent than the date of the data to compute. If it is the case, then the data is recomputed and the current date is stored. - + ** Electron-electron distances *** Get @@ -981,7 +952,7 @@ qmckl_exit_code qmckl_provide_ee_distance(qmckl_context context) qmckl_compute_ee_distance(context, ctx->electron.num, ctx->electron.walk_num, - ctx->electron.coord_new, + ctx->electron.coord_new.data, ctx->electron.ee_distance); if (rc != QMCKL_SUCCESS) { return rc; @@ -1088,7 +1059,7 @@ qmckl_exit_code qmckl_compute_ee_distance ( #+end_src *** Test - + #+begin_src python :results output :exports none import numpy as np @@ -1223,7 +1194,7 @@ qmckl_exit_code qmckl_provide_ee_distance_rescaled(qmckl_context context) ctx->electron.num, ctx->electron.rescale_factor_kappa_en, ctx->electron.walk_num, - ctx->electron.coord_new, + ctx->electron.coord_new.data, ctx->electron.ee_distance_rescaled); if (rc != QMCKL_SUCCESS) { return rc; @@ -1389,12 +1360,12 @@ assert(fabs(ee_distance_rescaled[elec_num*elec_num+1]-0.9985724058042633) < 1.e- #+end_src -** Electron-electron rescaled distance gradients and laplacian with respect to electron coords +** Electron-electron rescaled distance gradients and laplacian with respect to electron coords - The rescaled distances which is given as $R = (1 - \exp{-\kappa r})/\kappa$ + The rescaled distances which is given as $R = (1 - \exp{-\kappa r})/\kappa$ needs to be perturbed with respect to the electorn coordinates. - This data is stored in the ~ee_distance_rescaled_deriv_e~ tensor. The - The first three elements of this three index tensor ~[4][num][num]~ gives the + This data is stored in the ~ee_distance_rescaled_deriv_e~ tensor. The + The first three elements of this three index tensor ~[4][num][num]~ gives the derivatives in the x, y, and z directions $dx, dy, dz$ and the last index gives the Laplacian $\partial x^2 + \partial y^2 + \partial z^2$. @@ -1427,7 +1398,7 @@ qmckl_exit_code qmckl_get_electron_ee_distance_rescaled_deriv_e(qmckl_context co #+end_src *** Provide :noexport: - + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_provide_ee_distance_rescaled_deriv_e(qmckl_context context); #+end_src @@ -1469,7 +1440,7 @@ qmckl_exit_code qmckl_provide_ee_distance_rescaled_deriv_e(qmckl_context context ctx->electron.num, ctx->electron.rescale_factor_kappa_en, ctx->electron.walk_num, - ctx->electron.coord_new, + ctx->electron.coord_new.data, ctx->electron.ee_distance_rescaled_deriv_e); if (rc != QMCKL_SUCCESS) { return rc; @@ -1617,7 +1588,7 @@ rc = qmckl_get_electron_ee_distance_rescaled_deriv_e(context, ee_distance_rescal #+end_src ** Electron-electron potential - + ~ee_pot~ calculates the ~ee~ potential energy. \[ @@ -1654,7 +1625,7 @@ qmckl_exit_code qmckl_get_electron_ee_potential(qmckl_context context, double* c return QMCKL_SUCCESS; } #+end_src - + *** Provide :noexport: #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none @@ -1783,7 +1754,7 @@ end function qmckl_compute_ee_potential_f const int64_t elec_num, const int64_t walk_num, const double* ee_distance, - double* const ee_pot ); + double* const ee_pot ); #+end_src #+CALL: generate_c_interface(table=qmckl_ee_potential_args,rettyp=get_value("CRetType"),fname=get_value("Name")) @@ -1910,7 +1881,7 @@ qmckl_exit_code qmckl_provide_en_distance(qmckl_context context) ctx->electron.num, ctx->nucleus.num, ctx->electron.walk_num, - ctx->electron.coord_new, + ctx->electron.coord_new.data, ctx->nucleus.coord.data, ctx->electron.en_distance); if (rc != QMCKL_SUCCESS) { @@ -2101,7 +2072,7 @@ assert(fabs(en_distance[1][0][1] - 3.1804527583077356) < 1.e-12); ** Electron-nucleus rescaled distances - ~en_distance_rescaled~ stores the matrix of the rescaled distances between + ~en_distance_rescaled~ stores the matrix of the rescaled distances between electrons and nuclei. \[ @@ -2187,7 +2158,7 @@ qmckl_exit_code qmckl_provide_en_distance_rescaled(qmckl_context context) ctx->nucleus.num, ctx->electron.rescale_factor_kappa_en, ctx->electron.walk_num, - ctx->electron.coord_new, + ctx->electron.coord_new.data, ctx->nucleus.coord.data, ctx->electron.en_distance_rescaled); if (rc != QMCKL_SUCCESS) { @@ -2374,7 +2345,6 @@ assert (rc == QMCKL_SUCCESS); assert(fabs(en_distance_rescaled[0][0][0] - 0.9994721712909764) < 1.e-12); // (1,2,1) -printf("%f\n%f\n", en_distance_rescaled[0][1][0] , 0.9998448354439821); assert(fabs(en_distance_rescaled[0][1][0] - 0.9998448354439821) < 1.e-12); // (2,1,1) @@ -2393,10 +2363,10 @@ assert(fabs(en_distance_rescaled[1][0][1] - 0.9584331688679852) < 1.e-12); ** Electron-nucleus rescaled distance gradients and laplacian with respect to electron coords - The rescaled distances which is given as $R = (1 - \exp{-\kappa r})/\kappa$ + The rescaled distances which is given as $R = (1 - \exp{-\kappa r})/\kappa$ needs to be perturbed with respect to the nuclear coordinates. - This data is stored in the ~en_distance_rescaled_deriv_e~ tensor. The - The first three elements of this three index tensor ~[4][nucl_num][elec_num]~ gives the + This data is stored in the ~en_distance_rescaled_deriv_e~ tensor. The + The first three elements of this three index tensor ~[4][nucl_num][elec_num]~ gives the derivatives in the x, y, and z directions $dx, dy, dz$ and the last index gives the Laplacian $\partial x^2 + \partial y^2 + \partial z^2$. @@ -2476,7 +2446,7 @@ qmckl_exit_code qmckl_provide_en_distance_rescaled_deriv_e(qmckl_context context ctx->nucleus.num, ctx->electron.rescale_factor_kappa_en, ctx->electron.walk_num, - ctx->electron.coord_new, + ctx->electron.coord_new.data, ctx->nucleus.coord.data, ctx->electron.en_distance_rescaled_deriv_e); if (rc != QMCKL_SUCCESS) { @@ -2608,7 +2578,7 @@ qmckl_exit_code qmckl_compute_en_distance_rescaled_deriv_e ( #+end_src *** Test - + #+begin_src python :results output :exports none import numpy as np @@ -2671,7 +2641,7 @@ assert (rc == QMCKL_SUCCESS); where \(\mathcal{V}_{en}\) is the ~en~ potential, \[r_{iA}\] the ~en~ distance and \[Z_A\] is the nuclear charge. - + *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes @@ -2699,7 +2669,7 @@ qmckl_exit_code qmckl_get_electron_en_potential(qmckl_context context, double* c return QMCKL_SUCCESS; } #+end_src - + *** Provide :noexport: #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none @@ -2837,7 +2807,7 @@ end function qmckl_compute_en_potential_f const int64_t walk_num, const double* charge, const double* en_distance, - double* const en_pot ); + double* const en_pot ); #+end_src #+CALL: generate_c_interface(table=qmckl_en_potential_args,rettyp=get_value("CRetType"),fname=get_value("Name")) @@ -2873,7 +2843,7 @@ double en_pot[walk_num]; rc = qmckl_get_electron_en_potential(context, &(en_pot[0])); assert (rc == QMCKL_SUCCESS); #+end_src - + ** Generate initial coordinates *** Compute :noexport: @@ -2887,11 +2857,11 @@ subroutine draw_init_points integer :: iwalk logical, allocatable :: do_elec(:) integer :: acc_num - + real, allocatable :: xmin(:,:) - + integer :: i, j, k, l, kk - + real :: norm allocate (do_elec(elec_num), xmin(3,elec_num)) xmin = -huge(1.) @@ -2905,7 +2875,7 @@ subroutine draw_init_points call rinfo( irp_here, 'Norm : ', norm ) call rinfo( irp_here, 'mo_scale: ' , mo_scale ) mo_coef_transp = mo_coef_transp/norm - + double precision :: qmc_ranf real :: mo_max do i=1,elec_alpha_num @@ -2922,7 +2892,7 @@ subroutine draw_init_points xmin(2,i) = nucl_coord(l,2) xmin(3,i) = nucl_coord(l,3) enddo - + call iinfo(irp_here, 'Det num = ', det_num ) do k=1,elec_beta_num i = k+elec_alpha_num @@ -2939,10 +2909,10 @@ subroutine draw_init_points xmin(2,i) = nucl_coord(l,2) xmin(3,i) = nucl_coord(l,3) enddo - + call rinfo( irp_here, 'time step =', time_step ) do iwalk=1,walk_num - print *, 'Generating initial positions for walker', iwalk + print *, 'Generating initial positions for walker', iwalk acc_num = 0 do_elec = .True. integer :: iter @@ -2964,7 +2934,7 @@ subroutine draw_init_points TOUCH elec_coord re_compute = minval(nucl_elec_dist(1:nucl_num,1:elec_num)) enddo - + do i=1,elec_alpha_num if (do_elec(i)) then if ( mo_value_transp(i,i)**2 >= qmc_ranf()) then @@ -2973,7 +2943,7 @@ subroutine draw_init_points endif endif enddo - + do i=1,elec_beta_num if (do_elec(i+elec_alpha_num)) then if ( mo_value_transp(i,i+elec_alpha_num)**2 >= qmc_ranf()) then @@ -2982,9 +2952,9 @@ subroutine draw_init_points endif endif enddo - + enddo - + do l=1,3 do i=1,elec_num+1 elec_coord_full(i,l,iwalk) = elec_coord(i,l) @@ -2997,7 +2967,7 @@ subroutine draw_init_points endif SOFT_TOUCH elec_coord elec_coord_full deallocate (do_elec, xmin) - + end # end_src diff --git a/org/qmckl_jastrow.org b/org/qmckl_jastrow.org index 011c55c..ad6826f 100644 --- a/org/qmckl_jastrow.org +++ b/org/qmckl_jastrow.org @@ -3260,7 +3260,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_e_deriv_e(qmckl_context context) ctx->electron.num, ctx->jastrow.cord_num, ctx->electron.rescale_factor_kappa_ee, - ctx->electron.coord_new, + ctx->electron.coord_new.data, ctx->electron.ee_distance, ctx->jastrow.een_rescaled_e, ctx->jastrow.een_rescaled_e_deriv_e); @@ -3924,7 +3924,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_n_deriv_e(qmckl_context context) ctx->nucleus.num, ctx->jastrow.cord_num, ctx->electron.rescale_factor_kappa_en, - ctx->electron.coord_new, + ctx->electron.coord_new.data, ctx->nucleus.coord.data, ctx->electron.en_distance, ctx->jastrow.een_rescaled_n, diff --git a/org/qmckl_point.org b/org/qmckl_point.org index 7212ffe..d0f4614 100644 --- a/org/qmckl_point.org +++ b/org/qmckl_point.org @@ -77,10 +77,11 @@ int main() { The following data stored in the context: - | Variable | Type | Description | - |----------+----------------+------------------------| - | ~num~ | ~int64_t~ | Total number of points | - | ~coord~ | ~qmckl_matrix~ | ~num~ \times 3 matrix3 | + | Variable | Type | Description | + |----------+----------------+-------------------------------------------| + | ~num~ | ~int64_t~ | Total number of points | + | ~date~ | ~uint64_t~ | Last modification date of the coordinates | + | ~coord~ | ~qmckl_matrix~ | ~num~ \times 3 matrix3 | We consider that the matrix is stored 'transposed' and 'normal' corresponds to the 3 \times ~num~ matrix. @@ -90,15 +91,16 @@ int main() { #+begin_src c :comments org :tangle (eval h_private_type) typedef struct qmckl_point_struct { int64_t num; + uint64_t date; qmckl_matrix coord; } qmckl_point_struct; #+end_src - #+begin_src c :comments org :tangle (eval h_private_func) + #+begin_src c :comments org :tangle (eval h_private_func) qmckl_exit_code qmckl_init_point(qmckl_context context); #+end_src - + #+begin_src c :comments org :tangle (eval c) qmckl_exit_code qmckl_init_point(qmckl_context context) { @@ -109,16 +111,7 @@ qmckl_exit_code qmckl_init_point(qmckl_context context) { qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; - mem_info.size = sizeof(qmckl_point_struct); - ctx->point = (qmckl_point_struct*) qmckl_malloc(context, mem_info); - if (ctx->point == NULL) { - return qmckl_failwith( context, - QMCKL_ALLOCATION_FAILED, - "qmckl_init_point", - NULL); - } - memset(ctx->point, 0, sizeof(qmckl_point_struct)); + memset(&(ctx->point), 0, sizeof(qmckl_point_struct)); return QMCKL_SUCCESS; } @@ -139,7 +132,7 @@ qmckl_exit_code qmckl_get_point_num (const qmckl_context context, int64_t* const #+end_src Returns the number of points stored in the context. - + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_point_num (const qmckl_context context, int64_t* const num) { @@ -157,10 +150,9 @@ qmckl_get_point_num (const qmckl_context context, int64_t* const num) { qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - assert (ctx->point != NULL); - assert (ctx->point->num > (int64_t) 0); - ,*num = ctx->point->num; + assert (ctx->point.num >= (int64_t) 0); + ,*num = ctx->point.num; return QMCKL_SUCCESS; } #+end_src @@ -172,14 +164,14 @@ interface import implicit none - integer (c_int64_t) , intent(in) , value :: context - integer (c_int64_t) , intent(out) :: num + integer (c_int64_t) , intent(in) , value :: context + integer (c_int64_t) , intent(out) :: num end function end interface #+end_src *** Point coordinates - + #+begin_src c :comments org :tangle (eval h_func) :exports none qmckl_exit_code qmckl_get_point(const qmckl_context context, const char transp, @@ -212,12 +204,11 @@ qmckl_get_point(const qmckl_context context, qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - assert (ctx->point != NULL); - int64_t point_num = ctx->point->num; + int64_t point_num = ctx->point.num; if (point_num == 0) return QMCKL_NOT_PROVIDED; - assert (ctx->point->coord.data != NULL); + assert (ctx->point.coord.data != NULL); if (size_max < 3*point_num) { return qmckl_failwith( context, @@ -229,16 +220,16 @@ qmckl_get_point(const qmckl_context context, qmckl_exit_code rc; if (transp == 'N') { qmckl_matrix At = qmckl_matrix_alloc( context, 3, point_num); - rc = qmckl_transpose( context, ctx->point->coord, At ); + rc = qmckl_transpose( context, ctx->point.coord, At ); if (rc != QMCKL_SUCCESS) return rc; rc = qmckl_double_of_matrix( context, At, coord, size_max); if (rc != QMCKL_SUCCESS) return rc; rc = qmckl_matrix_free(context, At); } else { - rc = qmckl_double_of_matrix( context, ctx->point->coord, coord, size_max); + rc = qmckl_double_of_matrix( context, ctx->point.coord, coord, size_max); } if (rc != QMCKL_SUCCESS) return rc; - + return rc; } @@ -259,23 +250,23 @@ interface end interface #+end_src - + ** Initialization functions When the data is set in the context, if the arrays are large enough, we overwrite the data contained in them. - + To set the data relative to the points in the context, one of the - following functions need to be called. - + following functions need to be called. + #+begin_src c :comments org :tangle (eval h_func) qmckl_exit_code qmckl_set_point (qmckl_context context, const char transp, const double* coord, - const int64_t size_max); + const int64_t num); #+end_src - Copy a sequence of (x,y,z) into the context. + Copy a sequence of ~num~ points $(x,y,z)$ into the context. #+begin_src c :comments org :tangle (eval c) :noweb yes qmckl_exit_code @@ -284,7 +275,7 @@ qmckl_set_point (qmckl_context context, const double* coord, const int64_t num) { - + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } @@ -295,28 +286,27 @@ qmckl_set_point (qmckl_context context, "qmckl_set_point", "transp should be 'N' or 'T'"); } - + if (coord == NULL) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_set_point", "coord is a NULL pointer"); } - + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - assert (ctx->point != NULL); qmckl_exit_code rc; - if (ctx->point->num < num) { + if (ctx->point.num < num) { - if (ctx->point->coord.data != NULL) { - rc = qmckl_matrix_free(context, ctx->point->coord); - return rc; + if (ctx->point.coord.data != NULL) { + rc = qmckl_matrix_free(context, ctx->point.coord); + assert (rc == QMCKL_SUCCESS); } - ctx->point->coord = qmckl_matrix_alloc(context, num, 3); - if (ctx->point->coord.data == NULL) { + ctx->point.coord = qmckl_matrix_alloc(context, num, 3); + if (ctx->point.coord.data == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_set_point", @@ -325,18 +315,22 @@ qmckl_set_point (qmckl_context context, }; - ctx->point->num = num; + ctx->point.num = num; if (transp == 'T') { - memcpy(ctx->point->coord.data, coord, 3*num*sizeof(double)); + memcpy(ctx->point.coord.data, coord, 3*num*sizeof(double)); } else { for (int64_t i=0 ; ipoint->coord, i, 0) = coord[3*i ]; - qmckl_mat(ctx->point->coord, i, 1) = coord[3*i+1]; - qmckl_mat(ctx->point->coord, i, 2) = coord[3*i+2]; + qmckl_mat(ctx->point.coord, i, 0) = coord[3*i ]; + qmckl_mat(ctx->point.coord, i, 1) = coord[3*i+1]; + qmckl_mat(ctx->point.coord, i, 2) = coord[3*i+2]; } } + /* Increment the date of the context */ + ctx->date += 1UL; + ctx->point.date = ctx->date; + return QMCKL_SUCCESS; }