diff --git a/org/qmckl_ao.org b/org/qmckl_ao.org index 5ba081b..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; @@ -2258,10 +2258,10 @@ qmckl_exit_code rc; rc = qmckl_set_nucleus_num (context, nucl_num); assert(rc == QMCKL_SUCCESS); -rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0])); +rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), 3*nucl_num); assert(rc == QMCKL_SUCCESS); -rc = qmckl_set_nucleus_charge(context, nucl_charge); +rc = qmckl_set_nucleus_charge(context, nucl_charge, nucl_num); assert(rc == QMCKL_SUCCESS); assert(qmckl_nucleus_provided(context)); @@ -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,11 +3224,11 @@ 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->nucleus.coord, + ctx->point.coord.data, + ctx->nucleus.coord.data, ctx->ao_basis.exponent, ctx->ao_basis.primitive_vgl); } else { @@ -3296,7 +3296,7 @@ print ( "[7][4][26] : %e"% lf(a,x,y)) #+begin_src c :tangle (eval c_test) :exports none { -#define walk_num chbrclf_walk_num +#define walk_num 1 // chbrclf_walk_num #define elec_num chbrclf_elec_num #define prim_num chbrclf_prim_num @@ -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,14 +3615,14 @@ 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->nucleus.coord, + ctx->point.coord.data, + ctx->nucleus.coord.data, ctx->ao_basis.exponent, ctx->ao_basis.coefficient_normalized, ctx->ao_basis.shell_vgl); @@ -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 ), @@ -3710,7 +3703,7 @@ print ( "[1][4][26] : %25.15e"% lf(a,x,y)) #+begin_src c :tangle (eval c_test) :exports none { -#define walk_num chbrclf_walk_num +#define walk_num 1 // chbrclf_walk_num #define elec_num chbrclf_elec_num #define shell_num chbrclf_shell_num @@ -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,10 +4618,10 @@ 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->nucleus.coord, + ctx->point.coord.data, + ctx->nucleus.coord.data, ctx->ao_basis.nucleus_index, ctx->ao_basis.nucleus_shell_num, ctx->ao_basis.nucleus_range, @@ -4735,7 +4721,7 @@ print ( "[1][26][224] : %25.15e"%(df(a,x,y,1)* (x[2] - y[2]) * (x[2] - y[2])) ) #+begin_src c :tangle (eval c_test) :exports none { -#define walk_num chbrclf_walk_num +#define walk_num 1 // chbrclf_walk_num #define elec_num chbrclf_elec_num #define shell_num chbrclf_shell_num #define ao_num chbrclf_ao_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_blas.org b/org/qmckl_blas.org index 4a03a04..998e6ef 100644 --- a/org/qmckl_blas.org +++ b/org/qmckl_blas.org @@ -38,8 +38,9 @@ #include "qmckl.h" #include "qmckl_context_private_type.h" #include "qmckl_memory_private_type.h" -#include "qmckl_memory_private_func.h" #include "qmckl_blas_private_type.h" + +#include "qmckl_memory_private_func.h" #include "qmckl_blas_private_func.h" #+end_src @@ -379,23 +380,18 @@ qmckl_tensor_of_vector(const qmckl_vector vector, #+begin_src c :comments org :tangle (eval h_private_func) qmckl_vector -qmckl_vector_of_matrix(const qmckl_matrix matrix, - const int64_t size); +qmckl_vector_of_matrix(const qmckl_matrix matrix); #+end_src Reshapes a matrix into a vector. #+begin_src c :comments org :tangle (eval c) qmckl_vector -qmckl_vector_of_matrix(const qmckl_matrix matrix, - const int64_t size) +qmckl_vector_of_matrix(const qmckl_matrix matrix) { - /* Always true */ - assert (matrix.size[0] * matrix.size[1] == size); - qmckl_vector result; - result.size = size; + result.size = matrix.size[0] * matrix.size[1]; result.data = matrix.data; return result; @@ -438,27 +434,23 @@ qmckl_tensor_of_matrix(const qmckl_matrix matrix, #+begin_src c :comments org :tangle (eval h_private_func) qmckl_vector -qmckl_vector_of_tensor(const qmckl_tensor tensor, - const int64_t size); +qmckl_vector_of_tensor(const qmckl_tensor tensor); #+end_src Reshapes a tensor into a vector. #+begin_src c :comments org :tangle (eval c) qmckl_vector -qmckl_vector_of_tensor(const qmckl_tensor tensor, - const int64_t size) +qmckl_vector_of_tensor(const qmckl_tensor tensor) { - /* Always true */ - int64_t prod_size = (int64_t) 1; - for (int64_t i=0 ; i (int64_t) 0); + assert (target != NULL); + assert (size_max > (int64_t) 0); + assert (size_max >= vector.size); + for (int64_t i=0 ; isize[0], matrix->size[1]); + return rc; +} + #+end_src + + + #+begin_src c :comments org :tangle (eval h_private_func) +qmckl_exit_code +qmckl_tensor_of_double(const qmckl_context context, + const double* target, + const int64_t size_max, + qmckl_tensor* tensor); + #+end_src + + Converts a matrix to a ~double*~. + + #+begin_src c :comments org :tangle (eval c) +qmckl_exit_code +qmckl_tensor_of_double(const qmckl_context context, + const double* target, + const int64_t size_max, + qmckl_tensor* tensor) +{ + qmckl_vector vector = qmckl_vector_of_tensor(*tensor); + qmckl_exit_code rc = + qmckl_vector_of_double(context, target, size_max, &vector); + *tensor = qmckl_tensor_of_vector(vector, tensor->order, tensor->size); + return rc; +} + #+end_src + + + + ** Tests #+begin_src c :comments link :tangle (eval c_test) @@ -534,7 +706,7 @@ qmckl_matrix_of_tensor(const qmckl_tensor tensor, for (int64_t i=0 ; isize[0] = A.size[0]; + C->size[1] = B.size[1]; + rc = qmckl_dgemm (context, 'N', 'N', + C->size[0], C->size[1], A.size[1], + alpha, + A.data, A.size[0], + B.data, B.size[0], + beta, + C->data, C->size[0]); + break; + case 1: + if (A.size[0] != B.size[0]) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_matmul", + "A and B have incompatible dimensions"); + } + C->size[0] = A.size[1]; + C->size[1] = B.size[1]; + rc = qmckl_dgemm (context, 'T', 'N', + C->size[0], C->size[1], A.size[0], + alpha, + A.data, A.size[0], + B.data, B.size[0], + beta, + C->data, C->size[0]); + break; + case 2: + if (A.size[1] != B.size[1]) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_matmul", + "A and B have incompatible dimensions"); + } + C->size[0] = A.size[0]; + C->size[1] = B.size[0]; + rc = qmckl_dgemm (context, 'N', 'T', + C->size[0], C->size[1], A.size[1], + alpha, + A.data, A.size[0], + B.data, B.size[0], + beta, + C->data, C->size[0]); + break; + case 3: + if (A.size[0] != B.size[1]) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_matmul", + "A and B have incompatible dimensions"); + } + C->size[0] = A.size[1]; + C->size[1] = B.size[0]; + rc = qmckl_dgemm (context, 'T', 'T', + C->size[0], C->size[1], A.size[0], + alpha, + A.data, A.size[0], + B.data, B.size[0], + beta, + C->data, C->size[0]); + break; + } + return rc; +} + #+end_src + + +*** TODO Test :noexport: + ** ~qmckl_adjugate~ Given a matrix $\mathbf{A}$, the adjugate matrix @@ -1693,6 +2043,90 @@ qmckl_exit_code test_qmckl_adjugate(qmckl_context context); assert(QMCKL_SUCCESS == test_qmckl_adjugate(context)); #+end_src +** ~qmckl_transpose~ + + Transposes a matrix: $A^\dagger_{ji} = A_{ij}$. + + | Variable | Type | In/Out | Description | + |----------+---------------+--------+-------------------| + | context | qmckl_context | in | Global state | + | A | qmckl_matrix | in | Input matrix | + | At | qmckl_matrix | out | Transposed matrix | + + #+begin_src c :tangle (eval h_private_func) :comments org +qmckl_exit_code +qmckl_transpose (qmckl_context context, + const qmckl_matrix A, + qmckl_matrix At ); + #+end_src + + + #+begin_src c :tangle (eval c) :comments org +qmckl_exit_code +qmckl_transpose (qmckl_context context, + const qmckl_matrix A, + qmckl_matrix At ) +{ + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_INVALID_CONTEXT; + } + + if (A.size[0] < 1) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_transpose", + "Invalid size for A"); + } + + if (At.data == NULL) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_transpose", + "Output matrix not allocated"); + } + + if (At.size[0] != A.size[1] || At.size[1] != A.size[0]) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_transpose", + "Invalid size for At"); + } + + for (int64_t j=0 ; j +#include #include #ifdef HAVE_CONFIG_H #include "config.h" @@ -68,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[walk_num][3][num]~ | New set of electron coordinates | - | ~coord_old~ | ~double[walk_num][3][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 | @@ -107,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 @@ -145,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) { @@ -168,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 @@ -203,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 @@ -394,27 +395,34 @@ 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 | - |---------+---------------------------+---------------------------| - | C | ~[walk_num][elec_num][3]~ | ~[walk_num][3][elec_num]~ | - | Fortran | ~(3,elec_num,walk_num)~ | ~(elec_num,3,walk_num)~ | + | | Normal | Transposed | + |---------+--------------------------+--------------------------| + | C | ~[walk_num*elec_num][3]~ | ~[3][walk_num*elec_num]~ | + | Fortran | ~(3,walk_num*elec_num)~ | ~(walk_num*elec_num, 3)~ | #+begin_src c :comments org :tangle (eval h_func) :exports none -qmckl_exit_code qmckl_get_electron_coord (const qmckl_context context, const char transp, double* const coord); +qmckl_exit_code +qmckl_get_electron_coord (const qmckl_context context, + const char transp, + double* const coord, + 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, const char transp, double* const coord) { - - if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { - return QMCKL_INVALID_CONTEXT; - } - +qmckl_get_electron_coord (const qmckl_context context, + const char transp, + double* const coord, + const int64_t size_max) +{ if (transp != 'N' && transp != 'T') { return qmckl_failwith( context, QMCKL_INVALID_ARG_2, @@ -429,42 +437,32 @@ qmckl_get_electron_coord (const qmckl_context context, const char transp, double "coord is a null pointer"); } + if (size_max <= 0) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_4, + "qmckl_get_electron_coord", + "size_max should be > 0"); + } + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_INVALID_CONTEXT; + } + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - if ( !(ctx->electron.provided) ) { + if (!ctx->electron.provided) { return qmckl_failwith( context, QMCKL_NOT_PROVIDED, "qmckl_get_electron_coord", - "electron data is not provided"); + NULL); } - int64_t elec_num = ctx->electron.num; - int64_t walk_num = ctx->electron.walk_num; + 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]); - assert (ctx->electron.coord_new != NULL); - - double* ptr1 = ctx->electron.coord_new; - double* ptr2 = coord; - - if (transp == 'N') { - - for (int64_t i=0 ; ielectron.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", @@ -521,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", @@ -530,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; @@ -636,18 +640,19 @@ 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 + interface integer(c_int32_t) function qmckl_set_electron_walk_num(context, walk_num) bind(C) use, intrinsic :: iso_c_binding 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 @@ -659,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. @@ -686,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, @@ -698,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, @@ -709,41 +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", - "destination array is 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; @@ -758,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 @@ -865,9 +847,10 @@ assert(rc == QMCKL_SUCCESS); double elec_coord2[walk_num*3*elec_num]; -rc = qmckl_get_electron_coord (context, 'N', elec_coord2); +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 ; ++i) { +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] ); } @@ -882,7 +865,7 @@ for (int64_t i=0 ; i<3*elec_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 @@ -969,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; @@ -995,7 +978,7 @@ qmckl_exit_code qmckl_provide_ee_distance(qmckl_context context) | ~context~ | ~qmckl_context~ | in | Global state | | ~elec_num~ | ~int64_t~ | in | Number of electrons | | ~walk_num~ | ~int64_t~ | in | Number of walkers | - | ~coord~ | ~double[walk_num][3][elec_num]~ | in | Electron coordinates | + | ~coord~ | ~double[3][walk_num][elec_num]~ | in | Electron coordinates | | ~ee_distance~ | ~double[walk_num][elec_num][elec_num]~ | out | Electron-electron distances | #+begin_src f90 :comments org :tangle (eval f) :noweb yes @@ -1006,10 +989,11 @@ integer function qmckl_compute_ee_distance_f(context, elec_num, walk_num, coord, integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: elec_num integer*8 , intent(in) :: walk_num - double precision , intent(in) :: coord(elec_num,3,walk_num) + double precision , intent(in) :: coord(elec_num,walk_num,3) double precision , intent(out) :: ee_distance(elec_num,elec_num,walk_num) - integer*8 :: k + integer*8 :: k, i, j + double precision :: x, y, z info = QMCKL_SUCCESS @@ -1030,8 +1014,8 @@ integer function qmckl_compute_ee_distance_f(context, elec_num, walk_num, coord, do k=1,walk_num info = qmckl_distance(context, 'T', 'T', elec_num, elec_num, & - coord(1,1,k), elec_num, & - coord(1,1,k), elec_num, & + coord(1,k,1), elec_num * walk_num, & + coord(1,k,1), elec_num * walk_num, & ee_distance(1,1,k), elec_num) if (info /= QMCKL_SUCCESS) then exit @@ -1075,7 +1059,7 @@ qmckl_exit_code qmckl_compute_ee_distance ( #+end_src *** Test - + #+begin_src python :results output :exports none import numpy as np @@ -1210,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; @@ -1274,8 +1258,8 @@ integer function qmckl_compute_ee_distance_rescaled_f(context, elec_num, rescale do k=1,walk_num info = qmckl_distance_rescaled(context, 'T', 'T', elec_num, elec_num, & - coord(1,1,k), elec_num, & - coord(1,1,k), elec_num, & + coord(1,k,1), elec_num * walk_num, & + coord(1,k,1), elec_num * walk_num, & ee_distance_rescaled(1,1,k), elec_num, rescale_factor_kappa_ee) if (info /= QMCKL_SUCCESS) then exit @@ -1376,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$. @@ -1414,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 @@ -1456,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; @@ -1604,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. \[ @@ -1641,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 @@ -1770,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")) @@ -1897,8 +1881,8 @@ 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->nucleus.coord, + ctx->electron.coord_new.data, + ctx->nucleus.coord.data, ctx->electron.en_distance); if (rc != QMCKL_SUCCESS) { return rc; @@ -1938,7 +1922,7 @@ integer function qmckl_compute_en_distance_f(context, elec_num, nucl_num, walk_n 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,walk_num,3) double precision , intent(in) :: nucl_coord(nucl_num,3) double precision , intent(out) :: en_distance(elec_num,nucl_num,walk_num) @@ -1968,7 +1952,7 @@ integer function qmckl_compute_en_distance_f(context, elec_num, nucl_num, walk_n do k=1,walk_num info = qmckl_distance(context, 'T', 'T', elec_num, nucl_num, & - elec_coord(1,1,k), elec_num, & + elec_coord(1,k,1), elec_num * walk_num, & nucl_coord, nucl_num, & en_distance(1,1,k), elec_num) if (info /= QMCKL_SUCCESS) then @@ -2005,7 +1989,7 @@ qmckl_exit_code qmckl_compute_en_distance ( 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,walk_num,3) real (c_double ) , intent(in) :: nucl_coord(elec_num,3) real (c_double ) , intent(out) :: en_distance(elec_num,nucl_num,walk_num) @@ -2052,10 +2036,10 @@ assert(qmckl_electron_provided(context)); rc = qmckl_set_nucleus_num (context, nucl_num); assert(rc == QMCKL_SUCCESS); -rc = qmckl_set_nucleus_charge (context, charge); +rc = qmckl_set_nucleus_charge (context, charge, nucl_num); assert (rc == QMCKL_SUCCESS); -rc = qmckl_set_nucleus_coord (context, 'T', nucl_coord); +rc = qmckl_set_nucleus_coord (context, 'T', nucl_coord, 3*nucl_num); assert (rc == QMCKL_SUCCESS); assert(qmckl_nucleus_provided(context)); @@ -2088,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. \[ @@ -2103,6 +2087,7 @@ assert(fabs(en_distance[1][0][1] - 3.1804527583077356) < 1.e-12); qmckl_exit_code qmckl_get_electron_en_distance_rescaled(qmckl_context context, double* distance_rescaled); #+end_src + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_electron_en_distance_rescaled(qmckl_context context, double* distance_rescaled) { @@ -2173,8 +2158,8 @@ 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->nucleus.coord, + ctx->electron.coord_new.data, + ctx->nucleus.coord.data, ctx->electron.en_distance_rescaled); if (rc != QMCKL_SUCCESS) { return rc; @@ -2217,7 +2202,7 @@ integer function qmckl_compute_en_distance_rescaled_f(context, elec_num, nucl_nu integer*8 , intent(in) :: nucl_num double precision , intent(in) :: rescale_factor_kappa_en 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,walk_num,3) double precision , intent(in) :: nucl_coord(nucl_num,3) double precision , intent(out) :: en_distance_rescaled(elec_num,nucl_num,walk_num) @@ -2253,7 +2238,7 @@ integer function qmckl_compute_en_distance_rescaled_f(context, elec_num, nucl_nu do k=1,walk_num info = qmckl_distance_rescaled(context, 'T', 'T', elec_num, nucl_num, & - elec_coord(1,1,k), elec_num, & + elec_coord(1,k,1), elec_num*walk_num, & nucl_coord, nucl_num, & en_distance_rescaled(1,1,k), elec_num, rescale_factor_kappa_en) if (info /= QMCKL_SUCCESS) then @@ -2292,7 +2277,7 @@ qmckl_exit_code qmckl_compute_en_distance_rescaled ( integer (c_int64_t) , intent(in) , value :: nucl_num real (c_double ) , intent(in) , value :: rescale_factor_kappa_en 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,walk_num,3) real (c_double ) , intent(in) :: nucl_coord(elec_num,3) real (c_double ) , intent(out) :: en_distance_rescaled(elec_num,nucl_num,walk_num) @@ -2341,10 +2326,10 @@ assert(qmckl_electron_provided(context)); rc = qmckl_set_nucleus_num (context, nucl_num); assert(rc == QMCKL_SUCCESS); -rc = qmckl_set_nucleus_charge (context, charge); +rc = qmckl_set_nucleus_charge (context, charge, nucl_num); assert (rc == QMCKL_SUCCESS); -rc = qmckl_set_nucleus_coord (context, 'T', nucl_coord); +rc = qmckl_set_nucleus_coord (context, 'T', nucl_coord, 3*nucl_num); assert (rc == QMCKL_SUCCESS); assert(qmckl_nucleus_provided(context)); @@ -2378,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$. @@ -2461,8 +2446,8 @@ 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->nucleus.coord, + ctx->electron.coord_new.data, + ctx->nucleus.coord.data, ctx->electron.en_distance_rescaled_deriv_e); if (rc != QMCKL_SUCCESS) { return rc; @@ -2506,7 +2491,7 @@ integer function qmckl_compute_en_distance_rescaled_deriv_e_f(context, elec_num, integer*8 , intent(in) :: nucl_num double precision , intent(in) :: rescale_factor_kappa_en 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,walk_num,3) double precision , intent(in) :: nucl_coord(nucl_num,3) double precision , intent(out) :: en_distance_rescaled_deriv_e(elec_num,nucl_num,walk_num) @@ -2542,7 +2527,7 @@ integer function qmckl_compute_en_distance_rescaled_deriv_e_f(context, elec_num, do k=1,walk_num info = qmckl_distance_rescaled_deriv_e(context, 'T', 'T', elec_num, nucl_num, & - elec_coord(1,1,k), elec_num, & + elec_coord(1,k,1), elec_num*walk_num, & nucl_coord, nucl_num, & en_distance_rescaled_deriv_e(1,1,k), elec_num, rescale_factor_kappa_en) if (info /= QMCKL_SUCCESS) then @@ -2581,7 +2566,7 @@ qmckl_exit_code qmckl_compute_en_distance_rescaled_deriv_e ( integer (c_int64_t) , intent(in) , value :: nucl_num real (c_double ) , intent(in) , value :: rescale_factor_kappa_en 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,walk_num,3) real (c_double ) , intent(in) :: nucl_coord(elec_num,3) real (c_double ) , intent(out) :: en_distance_rescaled_deriv_e(elec_num,nucl_num,walk_num) @@ -2593,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 @@ -2611,10 +2596,10 @@ assert(rc == QMCKL_SUCCESS); rc = qmckl_set_nucleus_rescale_factor (context, nucl_rescale_factor_kappa); assert(rc == QMCKL_SUCCESS); -rc = qmckl_set_nucleus_charge (context, charge); +rc = qmckl_set_nucleus_charge (context, charge, nucl_num); assert (rc == QMCKL_SUCCESS); -rc = qmckl_set_nucleus_coord (context, 'T', nucl_coord); +rc = qmckl_set_nucleus_coord (context, 'T', nucl_coord, 3*nucl_num); assert (rc == QMCKL_SUCCESS); assert(qmckl_nucleus_provided(context)); @@ -2656,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 @@ -2684,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 @@ -2732,7 +2717,7 @@ qmckl_exit_code qmckl_provide_en_potential(qmckl_context context) ctx->electron.num, ctx->nucleus.num, ctx->electron.walk_num, - ctx->nucleus.charge, + ctx->nucleus.charge.data, ctx->electron.en_distance, ctx->electron.en_pot); if (rc != QMCKL_SUCCESS) { @@ -2822,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")) @@ -2858,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: @@ -2872,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.) @@ -2890,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 @@ -2907,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 @@ -2924,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 @@ -2949,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 @@ -2958,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 @@ -2967,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) @@ -2982,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 bc48594..ad6826f 100644 --- a/org/qmckl_jastrow.org +++ b/org/qmckl_jastrow.org @@ -1,4 +1,5 @@ #+TITLE: Jastrow Factor + #+SETUPFILE: ../tools/theme.setup #+INCLUDE: ../tools/lib.org @@ -103,26 +104,25 @@ int main() { computed data: - | Variable | Type | In/Out | Description | - |-------------------------------+---------------------------------------------------------------------+---------------------------------------------------------------------------------------------------------+---------------------------------| - | ~dim_cord_vect~ | ~int64_t~ | Number of unique C coefficients | | - | ~dim_cord_vect_date~ | ~uint64_t~ | Number of unique C coefficients | | - | ~asymp_jasb~ | ~double[2]~ | Asymptotic component | | - | ~asymp_jasb_date~ | ~uint64_t~ | Asymptotic component | | - | ~cord_vect_full~ | ~double[dim_cord_vect][nucl_num]~ | vector of non-zero coefficients | | - | ~cord_vect_full_date~ | ~uint64_t~ | Keep track of changes here | | - | ~lkpm_combined_index~ | ~int64_t[4][dim_cord_vect]~ | Transform l,k,p, and m into consecutive indices | | - | ~lkpm_combined_index_date~ | ~uint64_t~ | Transform l,k,p, and m into consecutive indices | | - | ~tmp_c~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]~ | out | vector of non-zero coefficients | - | ~dtmp_c~ | ~double[walk_num][elec_num][4][nucl_num][0:cord_num][0:cord_num-1]~ | vector of non-zero coefficients | | - | ~een_rescaled_e~ | ~double[walk_num][walk_num][elec_num][elec_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord | | - | ~een_rescaled_e_date~ | ~uint64_t~ | Keep track of the date of creation | | - | ~een_rescaled_n~ | ~double[walk_num][elec_num][nucl_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord | | - | ~een_rescaled_n_date~ | ~uint64_t~ | Keep track of the date of creation | | - | ~een_rescaled_e_deriv_e~ | ~double[walk_num][elec_num][4][elec_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons | | - | ~een_rescaled_e_deriv_e_date~ | ~uint64_t~ | Keep track of the date of creation | | - | ~een_rescaled_n_deriv_e~ | ~double[walk_num][elec_num][4][nucl_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons | | - | ~een_rescaled_n_deriv_e_date~ | ~uint64_t~ | Keep track of the date of creation | | + | Variable | Type | In/Out | Description | + |----------------------------+---------------------------------------------------------------------+-------------------------------------------------+---------------------------------| + | ~dim_cord_vect~ | ~int64_t~ | Number of unique C coefficients | | + | ~dim_cord_vect_date~ | ~uint64_t~ | Number of unique C coefficients | | + | ~asymp_jasb~ | ~double[2]~ | Asymptotic component | | + | ~asymp_jasb_date~ | ~uint64_t~ | Asymptotic component | | + | ~cord_vect_full~ | ~double[dim_cord_vect][nucl_num]~ | vector of non-zero coefficients | | + | ~cord_vect_full_date~ | ~uint64_t~ | Keep track of changes here | | + | ~lkpm_combined_index~ | ~int64_t[4][dim_cord_vect]~ | Transform l,k,p, and m into consecutive indices | | + | ~lkpm_combined_index_date~ | ~uint64_t~ | Transform l,k,p, and m into consecutive indices | | + | ~tmp_c~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]~ | out | vector of non-zero coefficients | + | ~dtmp_c~ | ~double[walk_num][elec_num][4][nucl_num][0:cord_num][0:cord_num-1]~ | vector of non-zero coefficients | | + + | ~een_rescaled_n~ | ~double[walk_num][elec_num][nucl_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord | | + | ~een_rescaled_n_date~ | ~uint64_t~ | Keep track of the date of creation | | + | ~een_rescaled_e_deriv_e~ | ~double[walk_num][elec_num][4][elec_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons | | + | ~een_rescaled_e_deriv_e_date~ | ~uint64_t~ | Keep track of the date of creation | | + | ~een_rescaled_n_deriv_e~ | ~double[walk_num][elec_num][4][nucl_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons | | + | ~een_rescaled_n_deriv_e_date~ | ~uint64_t~ | Keep track of the date of creation | | For H2O we have the following data: @@ -1093,7 +1093,7 @@ assert(rc == QMCKL_SUCCESS); double elec_coord2[walk_num*3*elec_num]; -rc = qmckl_get_electron_coord (context, 'N', elec_coord2); +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 ; ++i) { assert( elec_coord[i] == elec_coord2[i] ); @@ -1131,15 +1131,15 @@ assert(k == nucl_rescale_factor_kappa); double nucl_coord2[3*nucl_num]; -rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2); +rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2, 3*nucl_num); assert(rc == QMCKL_NOT_PROVIDED); -rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0])); +rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), 3*nucl_num); assert(rc == QMCKL_SUCCESS); assert(!qmckl_nucleus_provided(context)); -rc = qmckl_get_nucleus_coord (context, 'N', nucl_coord2); +rc = qmckl_get_nucleus_coord (context, 'N', nucl_coord2, nucl_num*3); assert(rc == QMCKL_SUCCESS); for (int64_t k=0 ; k<3 ; ++k) { for (int64_t i=0 ; ielectron.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,8 +3924,8 @@ 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->nucleus.coord, + ctx->electron.coord_new.data, + ctx->nucleus.coord.data, ctx->electron.en_distance, ctx->jastrow.een_rescaled_n, ctx->jastrow.een_rescaled_n_deriv_e); diff --git a/org/qmckl_local_energy.org b/org/qmckl_local_energy.org index 68c1525..1a1cacc 100644 --- a/org/qmckl_local_energy.org +++ b/org/qmckl_local_energy.org @@ -576,10 +576,10 @@ assert(rc == QMCKL_SUCCESS); rc = qmckl_set_nucleus_num (context, nucl_num); assert(rc == QMCKL_SUCCESS); -rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0])); +rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), nucl_num*3); assert(rc == QMCKL_SUCCESS); -rc = qmckl_set_nucleus_charge(context, nucl_charge); +rc = qmckl_set_nucleus_charge(context, nucl_charge, nucl_num); assert(rc == QMCKL_SUCCESS); assert(qmckl_nucleus_provided(context)); diff --git a/org/qmckl_mo.org b/org/qmckl_mo.org index a07bba4..462828d 100644 --- a/org/qmckl_mo.org +++ b/org/qmckl_mo.org @@ -720,10 +720,10 @@ assert(rc == QMCKL_SUCCESS); rc = qmckl_set_nucleus_num (context, nucl_num); assert(rc == QMCKL_SUCCESS); -rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0])); +rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), nucl_num*3); assert(rc == QMCKL_SUCCESS); -rc = qmckl_set_nucleus_charge(context, nucl_charge); +rc = qmckl_set_nucleus_charge(context, nucl_charge, nucl_num); assert(rc == QMCKL_SUCCESS); assert(qmckl_nucleus_provided(context)); diff --git a/org/qmckl_nucleus.org b/org/qmckl_nucleus.org index 75c91b9..a562618 100644 --- a/org/qmckl_nucleus.org +++ b/org/qmckl_nucleus.org @@ -67,42 +67,42 @@ int main() { The following data stored in the context: - |------------------------+----------------+-------------------------------------------| - | ~uninitialized~ | int32_t | Keeps bit set for uninitialized data | - | ~num~ | int64_t | Total number of nuclei | - | ~provided~ | bool | If true, ~nucleus~ is valid | - | ~charge~ | double[num] | Nuclear charges | - | ~coord~ | double[3][num] | Nuclear coordinates, in transposed format | - | ~coord_date~ | int64_t | Nuclear coordinates, date if modified | - | ~rescale_factor_kappa~ | double | The distance scaling factor | + |------------------------+--------------+-------------------------------------------| + | ~uninitialized~ | int32_t | Keeps bit set for uninitialized data | + | ~num~ | int64_t | Total number of nuclei | + | ~provided~ | bool | If true, ~nucleus~ is valid | + | ~charge~ | qmckl_vector | Nuclear charges | + | ~coord~ | qmckl_matrix | Nuclear coordinates, in transposed format | + | ~coord_date~ | int64_t | Nuclear coordinates, date if modified | + | ~rescale_factor_kappa~ | double | The distance scaling factor | Computed data: - |-----------------------------+------------------+------------------------------------------------------------| - | ~nn_distance~ | double[num][num] | Nucleus-nucleus distances | - | ~nn_distance_date~ | int64_t | Date when Nucleus-nucleus distances were computed | - | ~nn_distance_rescaled~ | double[num][num] | Nucleus-nucleus rescaled distances | - | ~nn_distance_rescaled_date~ | int64_t | Date when Nucleus-nucleus rescaled distances were computed | - | ~repulsion~ | double | Nuclear repulsion energy | - | ~repulsion_date~ | int64_t | Date when the nuclear repulsion energy was computed | + |-----------------------------+--------------+------------------------------------------------------------| + | ~nn_distance~ | qmckl_matrix | Nucleus-nucleus distances | + | ~nn_distance_date~ | int64_t | Date when Nucleus-nucleus distances were computed | + | ~nn_distance_rescaled~ | qmckl_matrix | Nucleus-nucleus rescaled distances | + | ~nn_distance_rescaled_date~ | int64_t | Date when Nucleus-nucleus rescaled distances were computed | + | ~repulsion~ | double | Nuclear repulsion energy | + | ~repulsion_date~ | int64_t | Date when the nuclear repulsion energy was computed | ** Data structure #+begin_src c :comments org :tangle (eval h_private_type) typedef struct qmckl_nucleus_struct { - int64_t num; - int64_t repulsion_date; - int64_t nn_distance_date; - int64_t nn_distance_rescaled_date; - int64_t coord_date; - double* coord; - double* charge; - double* nn_distance; - double* nn_distance_rescaled; - double repulsion; - double rescale_factor_kappa; - int32_t uninitialized; - bool provided; + int64_t num; + int64_t repulsion_date; + int64_t nn_distance_date; + int64_t nn_distance_rescaled_date; + int64_t coord_date; + qmckl_vector charge; + qmckl_matrix coord; + qmckl_matrix nn_distance; + qmckl_matrix nn_distance_rescaled; + double repulsion; + double rescale_factor_kappa; + int32_t uninitialized; + bool provided; } qmckl_nucleus_struct; #+end_src @@ -138,15 +138,6 @@ qmckl_exit_code qmckl_init_nucleus(qmckl_context context) { #+end_src ** Access functions - - #+begin_src c :comments org :tangle (eval h_func) :exports none -qmckl_exit_code qmckl_get_nucleus_num (const qmckl_context context, int64_t* const num); -qmckl_exit_code qmckl_get_nucleus_charge (const qmckl_context context, double* const charge); -qmckl_exit_code qmckl_get_nucleus_coord (const qmckl_context context, const char transp, double* const coord); - -qmckl_exit_code qmckl_get_nucleus_rescale_factor (const qmckl_context context, double* const rescale_factor_kappa); - #+end_src - #+NAME:post #+begin_src c :exports none if ( (ctx->nucleus.uninitialized & mask) != 0) { @@ -154,6 +145,13 @@ if ( (ctx->nucleus.uninitialized & mask) != 0) { } #+end_src + + #+begin_src c :comments org :tangle (eval h_func) :exports none +qmckl_exit_code +qmckl_get_nucleus_num(const qmckl_context context, + int64_t* const num); + #+end_src + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_nucleus_num (const qmckl_context context, int64_t* const num) { @@ -187,10 +185,22 @@ qmckl_get_nucleus_num (const qmckl_context context, int64_t* const num) { return QMCKL_SUCCESS; } + #+end_src + #+begin_src c :comments org :tangle (eval h_func) :exports none qmckl_exit_code -qmckl_get_nucleus_charge (const qmckl_context context, double* const charge) { +qmckl_get_nucleus_charge(const qmckl_context context, + double* const charge, + const int64_t size_max); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code +qmckl_get_nucleus_charge (const qmckl_context context, + double* const charge, + const int64_t size_max) +{ if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_INVALID_CONTEXT; @@ -215,16 +225,31 @@ qmckl_get_nucleus_charge (const qmckl_context context, double* const charge) { "nucleus data is not provided"); } - assert (ctx->nucleus.charge != NULL); + assert (ctx->nucleus.charge.data != NULL); - int64_t nucl_num = ctx->nucleus.num; + if (ctx->nucleus.num > size_max) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_get_nucleus_charge", + "Array too small"); + } - memcpy(charge, ctx->nucleus.charge, nucl_num*sizeof(double)); + qmckl_exit_code rc; + rc = qmckl_double_of_vector(context, ctx->nucleus.charge, charge, size_max); - return QMCKL_SUCCESS; + return rc; } + #+end_src + + #+begin_src c :comments org :tangle (eval h_func) :exports none +qmckl_exit_code +qmckl_get_nucleus_rescale_factor(const qmckl_context context, + double* const rescale_factor_kappa); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_nucleus_rescale_factor (const qmckl_context context, double* const rescale_factor_kappa) @@ -250,10 +275,24 @@ qmckl_get_nucleus_rescale_factor (const qmckl_context context, return QMCKL_SUCCESS; } + #+end_src + #+begin_src c :comments org :tangle (eval h_func) :exports none qmckl_exit_code -qmckl_get_nucleus_coord (const qmckl_context context, const char transp, double* const coord) { +qmckl_get_nucleus_coord(const qmckl_context context, + const char transp, + double* const coord, + const int64_t size_max); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code +qmckl_get_nucleus_coord (const qmckl_context context, + const char transp, + double* const coord, + const int64_t size_max) +{ if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_INVALID_CONTEXT; @@ -285,25 +324,21 @@ qmckl_get_nucleus_coord (const qmckl_context context, const char transp, double* "nucleus data is not provided"); } - int64_t nucl_num = ctx->nucleus.num; + assert (ctx->nucleus.coord.data != NULL); - assert (ctx->nucleus.coord != NULL); + qmckl_exit_code rc; if (transp == 'N') { - qmckl_exit_code rc; - - rc = qmckl_transpose(context, nucl_num, 3, - ctx->nucleus.coord, nucl_num, - coord, 3); + qmckl_matrix At = qmckl_matrix_alloc(context, 3, ctx->nucleus.coord.size[0]); + rc = qmckl_transpose(context, ctx->nucleus.coord, At); if (rc != QMCKL_SUCCESS) return rc; - + rc = qmckl_double_of_matrix(context, At, coord, size_max); + qmckl_matrix_free(context, At); } else { - - memcpy(coord, ctx->nucleus.coord, 3*nucl_num*sizeof(double)); - + rc = qmckl_double_of_matrix(context, ctx->nucleus.coord, coord, size_max); } - return QMCKL_SUCCESS; + return rc; } #+end_src @@ -330,51 +365,6 @@ bool qmckl_nucleus_provided(const qmckl_context context) { ** Initialization functions - To set the data relative to the nuclei in the context, the - following functions need to be called. - - #+begin_src c :comments org :tangle (eval h_func) -qmckl_exit_code qmckl_set_nucleus_num (qmckl_context context, const int64_t num); -qmckl_exit_code qmckl_set_nucleus_charge (qmckl_context context, const double* charge); -qmckl_exit_code qmckl_set_nucleus_coord (qmckl_context context, const char transp, const double* coord); - -qmckl_exit_code qmckl_set_nucleus_rescale_factor (qmckl_context context, const double rescale_factor_kappa); - #+end_src - - #+begin_src f90 :tangle (eval fh_func) :comments org :exports none -interface - integer(c_int32_t) function qmckl_set_nucleus_num(context, num) & - bind(C) - use, intrinsic :: iso_c_binding - import - implicit none - integer (c_int64_t) , intent(in) , value :: context - integer (c_int64_t) , intent(in) , value :: num - end function -end interface -interface - integer(c_int32_t) function qmckl_set_nucleus_charge(context, charge) & - bind(C) - use, intrinsic :: iso_c_binding - import - implicit none - integer (c_int64_t) , intent(in) , value :: context - real (c_double) , intent(in) :: charge(*) - end function -end interface -interface - integer(c_int32_t) function qmckl_set_nucleus_coord(context, transp, coord) & - bind(C) - use, intrinsic :: iso_c_binding - import - implicit none - integer (c_int64_t) , intent(in) , value :: context - character(c_char) , intent(in) , value :: transp - real (c_double) , intent(in) :: coord(*) - end function -end interface - #+end_src - #+NAME:pre2 #+begin_src c :exports none if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { @@ -392,12 +382,22 @@ ctx->nucleus.provided = (ctx->nucleus.uninitialized == 0); return QMCKL_SUCCESS; #+end_src + To set the data relative to the nuclei in the context, the + following functions need to be called. - To set the number of nuclei, use + #+begin_src c :comments org :tangle (eval h_func) +qmckl_exit_code +qmckl_set_nucleus_num(qmckl_context context, + const int64_t num); + #+end_src + + Sets the number of nuclei. #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code -qmckl_set_nucleus_num(qmckl_context context, const int64_t num) { +qmckl_set_nucleus_num(qmckl_context context, + const int64_t num) +{ <> if (num <= 0) { @@ -415,11 +415,35 @@ qmckl_set_nucleus_num(qmckl_context context, const int64_t num) { } #+end_src - The following function sets the nuclear charges of all the atoms. + #+begin_src f90 :tangle (eval fh_func) :comments org :exports none +interface + integer(c_int32_t) function qmckl_set_nucleus_num(context, num) & + bind(C) + use, intrinsic :: iso_c_binding + import + implicit none + integer (c_int64_t) , intent(in) , value :: context + integer (c_int64_t) , intent(in) , value :: num + end function +end interface + #+end_src + + + #+begin_src c :comments org :tangle (eval h_func) +qmckl_exit_code +qmckl_set_nucleus_charge(qmckl_context context, + const double* charge, + const int64_t size_max); + #+end_src + + Sets the nuclear charges of all the atoms. #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code -qmckl_set_nucleus_charge(qmckl_context context, const double* charge) { +qmckl_set_nucleus_charge(qmckl_context context, + const double* charge, + const int64_t size_max) +{ <> if (charge == NULL) { @@ -437,34 +461,126 @@ qmckl_set_nucleus_charge(qmckl_context context, const double* charge) { rc = qmckl_get_nucleus_num(context, &num); if (rc != QMCKL_SUCCESS) return rc; - if (ctx->nucleus.charge != NULL) { - qmckl_free(context, ctx->nucleus.charge); - ctx->nucleus.charge= NULL; - } - - qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; - mem_info.size = num*sizeof(double); - - assert (ctx->nucleus.charge == NULL); - ctx->nucleus.charge = (double*) qmckl_malloc(context, mem_info); - if (ctx->nucleus.charge == NULL) { + if (num > size_max) { return qmckl_failwith( context, - QMCKL_ALLOCATION_FAILED, + QMCKL_INVALID_ARG_3, "qmckl_set_nucleus_charge", - NULL); + "Array too small"); } - ctx->nucleus.charge= memcpy(ctx->nucleus.charge, charge, num*sizeof(double)); - assert (ctx->nucleus.charge != NULL); + + ctx->nucleus.charge = qmckl_vector_alloc(context, num); + rc = qmckl_vector_of_double(context, charge, num, &(ctx->nucleus.charge)); <> } #+end_src - The following function sets the rescale parameter for the nuclear distances. + #+begin_src f90 :tangle (eval fh_func) :comments org :exports none +interface + integer(c_int32_t) function qmckl_set_nucleus_charge(context, charge, size_max) & + bind(C) + use, intrinsic :: iso_c_binding + import + implicit none + integer (c_int64_t) , intent(in) , value :: context + real (c_double) , intent(in) :: charge(*) + integer (c_int64_t) , intent(in) , value :: size_max + end function +end interface + #+end_src + + + + #+begin_src c :comments org :tangle (eval h_func) +qmckl_exit_code +qmckl_set_nucleus_coord(qmckl_context context, + const char transp, + const double* coord, + const int64_t size_max); + #+end_src + + Sets the nuclear coordinates of all the atoms. The coordinates + are be given in atomic units. #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code -qmckl_set_nucleus_rescale_factor(qmckl_context context, const double rescale_factor_kappa) { +qmckl_set_nucleus_coord(qmckl_context context, + const char transp, + const double* coord, + const int64_t size_max) +{ + <> + + qmckl_exit_code rc; + + int32_t mask = 1 << 2; + + const int64_t nucl_num = (int64_t) ctx->nucleus.num; + + if (ctx->nucleus.coord.data != NULL) { + rc = qmckl_matrix_free(context, ctx->nucleus.coord); + if (rc != QMCKL_SUCCESS) return rc; + } + + ctx->nucleus.coord = qmckl_matrix_alloc(context, nucl_num, 3); + + if (ctx->nucleus.coord.data == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_set_nucleus_coord", + NULL); + } + + if (size_max < 3*nucl_num) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_4, + "qmckl_set_nucleus_coord", + "Array too small"); + } + + if (transp == 'N') { + qmckl_matrix At; + At = qmckl_matrix_alloc(context, 3, nucl_num); + rc = qmckl_matrix_of_double(context, coord, 3*nucl_num, &At); + if (rc != QMCKL_SUCCESS) return rc; + rc = qmckl_transpose(context, At, ctx->nucleus.coord); + } else { + rc = qmckl_matrix_of_double(context, coord, nucl_num*3, &(ctx->nucleus.coord)); + } + if (rc != QMCKL_SUCCESS) return rc; + + <> +} + #+end_src + + #+begin_src f90 :tangle (eval fh_func) :comments org :exports none +interface + integer(c_int32_t) function qmckl_set_nucleus_coord(context, transp, coord, size_max) & + bind(C) + use, intrinsic :: iso_c_binding + import + implicit none + integer (c_int64_t) , intent(in) , value :: context + character(c_char) , intent(in) , value :: transp + real (c_double) , intent(in) :: coord(*) + integer (c_int64_t) , intent(in) , value :: size_max + end function +end interface + #+end_src + + #+begin_src c :comments org :tangle (eval h_func) +qmckl_exit_code +qmckl_set_nucleus_rescale_factor(qmckl_context context, + const double kappa); + #+end_src + + Sets the rescale parameter for the nuclear distances. + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code +qmckl_set_nucleus_rescale_factor(qmckl_context context, + const double rescale_factor_kappa) +{ <> if (rescale_factor_kappa <= 0.0) { @@ -480,50 +596,17 @@ qmckl_set_nucleus_rescale_factor(qmckl_context context, const double rescale_fac } #+end_src - The following function sets the nuclear coordinates of all the - atoms. The coordinates should be given in atomic units. - - #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -qmckl_exit_code -qmckl_set_nucleus_coord(qmckl_context context, const char transp, const double* coord) { - <> - - int64_t nucl_num = (int64_t) 0; - qmckl_exit_code rc; - - int32_t mask = 1 << 2; - - rc = qmckl_get_nucleus_num(context, &nucl_num); - if (rc != QMCKL_SUCCESS) return rc; - - if (ctx->nucleus.coord != NULL) { - qmckl_free(context, ctx->nucleus.coord); - ctx->nucleus.coord = NULL; - } - - qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; - mem_info.size = 3*nucl_num*sizeof(double); - - assert(ctx->nucleus.coord == NULL); - - ctx->nucleus.coord = (double*) qmckl_malloc(context, mem_info); - if (coord == NULL) { - return qmckl_failwith( context, - QMCKL_ALLOCATION_FAILED, - "qmckl_set_nucleus_coord", - NULL); - } - if (transp == 'N') { - rc = qmckl_transpose(context, 3, nucl_num, - coord, 3, - ctx->nucleus.coord, nucl_num); - if (rc != QMCKL_SUCCESS) return rc; - } else { - memcpy(ctx->nucleus.coord, coord, 3*nucl_num*sizeof(double)); - } - - <> -} + #+begin_src f90 :tangle (eval fh_func) :comments org :exports none +interface + integer(c_int32_t) function qmckl_set_rescale_factor(context, kappa) & + bind(C) + use, intrinsic :: iso_c_binding + import + implicit none + integer (c_int64_t) , intent(in) , value :: context + real (c_double) , intent(in) , value :: kappa + end function +end interface #+end_src ** Test @@ -568,15 +651,15 @@ assert(k == nucl_rescale_factor_kappa); double nucl_coord2[3*nucl_num]; -rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2); +rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2, 3*nucl_num); assert(rc == QMCKL_NOT_PROVIDED); -rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0])); +rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), 3*nucl_num); assert(rc == QMCKL_SUCCESS); assert(!qmckl_nucleus_provided(context)); -rc = qmckl_get_nucleus_coord (context, 'N', nucl_coord2); +rc = qmckl_get_nucleus_coord (context, 'N', nucl_coord2, 3*nucl_num); assert(rc == QMCKL_SUCCESS); for (size_t k=0 ; k<3 ; ++k) { for (size_t i=0 ; inucleus.num * ctx->nucleus.num; - memcpy(distance, ctx->nucleus.nn_distance, sze * sizeof(double)); + const int64_t sze = ctx->nucleus.num * ctx->nucleus.num; + if (sze > size_max) { + return qmckl_failwith(context, + QMCKL_INVALID_ARG_3, + "qmckl_get_nucleus_nn_distance", + "Array too small"); + } + rc = qmckl_double_of_matrix(context, ctx->nucleus.nn_distance, distance, size_max); - return QMCKL_SUCCESS; + return rc; } #+end_src #+begin_src f90 :tangle (eval fh_func) :comments org :exports none interface - integer(c_int32_t) function qmckl_get_nucleus_nn_distance(context, distance) & + integer(c_int32_t) function qmckl_get_nucleus_nn_distance(context, distance, size_max) & bind(C) use, intrinsic :: iso_c_binding import implicit none integer (c_int64_t) , intent(in) , value :: context real (c_double ) , intent(out) :: distance(*) + integer (c_int64_t) , intent(in) , value :: size_max end function end interface #+end_src @@ -678,26 +774,23 @@ qmckl_exit_code qmckl_provide_nn_distance(qmckl_context context) if (!ctx->nucleus.provided) return QMCKL_NOT_PROVIDED; /* Allocate array */ - if (ctx->nucleus.nn_distance == NULL) { + if (ctx->nucleus.nn_distance.data == NULL) { + ctx->nucleus.nn_distance = + qmckl_matrix_alloc(context, ctx->nucleus.num, ctx->nucleus.num); - qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; - mem_info.size = ctx->nucleus.num * ctx->nucleus.num * sizeof(double); - double* nn_distance = (double*) qmckl_malloc(context, mem_info); - - if (nn_distance == NULL) { + if (ctx->nucleus.nn_distance.data == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_nn_distance", NULL); } - ctx->nucleus.nn_distance = nn_distance; } qmckl_exit_code rc = qmckl_compute_nn_distance(context, ctx->nucleus.num, - ctx->nucleus.coord, - ctx->nucleus.nn_distance); + ctx->nucleus.coord.data, + ctx->nucleus.nn_distance.data); if (rc != QMCKL_SUCCESS) { return rc; } @@ -788,7 +881,7 @@ qmckl_exit_code qmckl_compute_nn_distance ( assert(qmckl_nucleus_provided(context)); double distance[nucl_num*nucl_num]; -rc = qmckl_get_nucleus_nn_distance(context, distance); +rc = qmckl_get_nucleus_nn_distance(context, distance, nucl_num*nucl_num); assert(distance[0] == 0.); assert(distance[1] == distance[nucl_num]); assert(fabs(distance[1]-2.070304721365169) < 1.e-12); @@ -800,11 +893,17 @@ assert(fabs(distance[1]-2.070304721365169) < 1.e-12); *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes -qmckl_exit_code qmckl_get_nucleus_nn_distance_rescaled(qmckl_context context, double* distance_rescaled); +qmckl_exit_code +qmckl_get_nucleus_nn_distance_rescaled(qmckl_context context, + double* distance_rescaled, + const int64_t size_max); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -qmckl_exit_code qmckl_get_nucleus_nn_distance_rescaled(qmckl_context context, double* distance_rescaled) +qmckl_exit_code +qmckl_get_nucleus_nn_distance_rescaled(qmckl_context context, + double* distance_rescaled, + const int64_t size_max) { /* Check input parameters */ if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { @@ -817,13 +916,35 @@ qmckl_exit_code qmckl_get_nucleus_nn_distance_rescaled(qmckl_context context, do qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - size_t sze = ctx->nucleus.num * ctx->nucleus.num; - memcpy(distance_rescaled, ctx->nucleus.nn_distance_rescaled, sze * sizeof(double)); - - return QMCKL_SUCCESS; + const int64_t sze = ctx->nucleus.num * ctx->nucleus.num; + if (sze > size_max) { + return qmckl_failwith(context, + QMCKL_INVALID_ARG_3, + "qmckl_get_nucleus_nn_distance_rescaled", + "Array too small"); + } + rc = qmckl_double_of_matrix(context, + ctx->nucleus.nn_distance_rescaled, + distance_rescaled, + size_max); + return rc; } #+end_src + #+begin_src f90 :tangle (eval fh_func) :comments org :exports none +interface + integer(c_int32_t) function qmckl_get_nucleus_nn_distance_rescaled(context, distance_rescaled, size_max) & + bind(C) + use, intrinsic :: iso_c_binding + import + implicit none + integer (c_int64_t) , intent(in) , value :: context + real (c_double ) , intent(out) :: distance_rescaled(*) + integer (c_int64_t) , intent(in) , value :: size_max + end function +end interface + #+end_src + *** Provide :noexport: #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none @@ -844,27 +965,24 @@ qmckl_exit_code qmckl_provide_nn_distance_rescaled(qmckl_context context) if (!ctx->nucleus.provided) return QMCKL_NOT_PROVIDED; /* Allocate array */ - if (ctx->nucleus.nn_distance_rescaled == NULL) { + if (ctx->nucleus.nn_distance_rescaled.data == NULL) { + ctx->nucleus.nn_distance_rescaled = + qmckl_matrix_alloc(context, ctx->nucleus.num, ctx->nucleus.num); - qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; - mem_info.size = ctx->nucleus.num * ctx->nucleus.num * sizeof(double); - double* nn_distance_rescaled = (double*) qmckl_malloc(context, mem_info); - - if (nn_distance_rescaled == NULL) { + if (ctx->nucleus.nn_distance_rescaled.data == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_nn_distance_rescaled", NULL); } - ctx->nucleus.nn_distance_rescaled = nn_distance_rescaled; } qmckl_exit_code rc = qmckl_compute_nn_distance_rescaled(context, ctx->nucleus.num, ctx->nucleus.rescale_factor_kappa, - ctx->nucleus.coord, - ctx->nucleus.nn_distance_rescaled); + ctx->nucleus.coord.data, + ctx->nucleus.nn_distance_rescaled.data); if (rc != QMCKL_SUCCESS) { return rc; } @@ -959,7 +1077,7 @@ qmckl_exit_code qmckl_compute_nn_distance_rescaled ( //assert(qmckl_nucleus_provided(context)); // //double distance[nucl_num*nucl_num]; -//rc = qmckl_get_nucleus_nn_distance(context, distance); +//rc = qmckl_get_nucleus_nn_distance(context, distance, nucl_num*nucl_num); //assert(distance[0] == 0.); //assert(distance[1] == distance[nucl_num]); //assert(fabs(distance[1]-2.070304721365169) < 1.e-12); @@ -975,11 +1093,11 @@ qmckl_exit_code qmckl_compute_nn_distance_rescaled ( *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes -qmckl_exit_code qmckl_get_nucleus_repulsion(qmckl_context context, double* energy); +qmckl_exit_code qmckl_get_nucleus_repulsion(qmckl_context context, double* const energy); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -qmckl_exit_code qmckl_get_nucleus_repulsion(qmckl_context context, double* energy) +qmckl_exit_code qmckl_get_nucleus_repulsion(qmckl_context context, double* const energy) { /* Check input parameters */ if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { @@ -1037,8 +1155,8 @@ qmckl_exit_code qmckl_provide_nucleus_repulsion(qmckl_context context) rc = qmckl_compute_nucleus_repulsion(context, ctx->nucleus.num, - ctx->nucleus.charge, - ctx->nucleus.nn_distance, + ctx->nucleus.charge.data, + ctx->nucleus.nn_distance.data, &(ctx->nucleus.repulsion)); if (rc != QMCKL_SUCCESS) { return rc; diff --git a/org/qmckl_point.org b/org/qmckl_point.org index e4683de..d0f4614 100644 --- a/org/qmckl_point.org +++ b/org/qmckl_point.org @@ -17,11 +17,14 @@ walkers. #ifndef QMCKL_POINT_HPT #define QMCKL_POINT_HPT #include +#include "qmckl_blas_private_type.h" #+end_src #+begin_src c :tangle (eval h_private_func) #ifndef QMCKL_POINT_HPF #define QMCKL_POINT_HPF +#include "qmckl_blas_private_type.h" +#include "qmckl_blas_private_func.h" #+end_src #+begin_src c :tangle (eval c_test) :noweb yes @@ -33,6 +36,8 @@ walkers. #endif #include "chbrclf.h" +#include "qmckl_blas_private_type.h" +#include "qmckl_blas_private_func.h" int main() { qmckl_context context; @@ -61,44 +66,41 @@ int main() { #include "qmckl.h" #include "qmckl_context_private_type.h" #include "qmckl_memory_private_type.h" -#include "qmckl_memory_private_func.h" +#include "qmckl_blas_private_type.h" + #include "qmckl_point_private_func.h" +#include "qmckl_memory_private_func.h" +#include "qmckl_blas_private_func.h" #+end_src * Context The following data stored in the context: - | Variable | Type | Description | - |-----------+---------------+------------------------| - | ~num~ | ~int64_t~ | Total number of points | - | ~coord_x~ | ~double[num]~ | X coordinates | - | ~coord_y~ | ~double[num]~ | Y coordinates | - | ~coord_z~ | ~double[num]~ | Z coordinates | + | 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 'transposed' and 'normal' storage follows the convention: - - | | Normal | Transposed | - |---------+------------------+------------------| - | C | ~[point_num][3]~ | ~[3][point_num]~ | - | Fortran | ~(3,point_num)~ | ~(point_num,3)~ | + We consider that the matrix is stored 'transposed' and 'normal' + corresponds to the 3 \times ~num~ matrix. ** Data structure #+begin_src c :comments org :tangle (eval h_private_type) typedef struct qmckl_point_struct { - double* coord_x; - double* coord_y; - double* coord_z; - int64_t num; + 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,16 +164,17 @@ 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, double* const coord, const int64_t size_max); #+end_src @@ -193,6 +186,7 @@ qmckl_exit_code qmckl_get_point(const qmckl_context context, #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_point(const qmckl_context context, + const char transp, double* const coord, const int64_t size_max) { @@ -210,13 +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_x != NULL); - assert (ctx->point->coord_y != NULL); - assert (ctx->point->coord_z != NULL); + assert (ctx->point.coord.data != NULL); if (size_max < 3*point_num) { return qmckl_failwith( context, @@ -225,212 +217,120 @@ qmckl_get_point(const qmckl_context context, "size_max too small"); } - - double * ptr = coord; - for (int64_t i=0 ; ipoint->coord_x[i]; ++ptr; - ,*ptr = ctx->point->coord_y[i]; ++ptr; - ,*ptr = ctx->point->coord_z[i]; ++ptr; + 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 ); + 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); } - - return QMCKL_SUCCESS; + if (rc != QMCKL_SUCCESS) return rc; + + return rc; } #+end_src #+begin_src f90 :comments org :tangle (eval fh_func) :noweb yes interface - integer(c_int32_t) function qmckl_get_point(context, coord, size_max) bind(C) + integer(c_int32_t) function qmckl_get_point(context, transp, coord, size_max) bind(C) use, intrinsic :: iso_c_binding import implicit none integer (c_int64_t) , intent(in) , value :: context + character(c_char) , intent(in) , value :: transp real (c_double ) , intent(out) :: coord(*) integer (c_int64_t) , intent(in) :: size_max end function end interface #+end_src - - #+begin_src c :comments org :tangle (eval h_func) :exports none -qmckl_exit_code qmckl_get_point_xyz (const qmckl_context context, - double* const coord_x, - double* const coord_y, - double* const coord_z, - const int64_t size_max); - #+end_src - - Returns the point coordinates in three different arrays, one for - each component x,y,z. - The pointers are assumed to point on a memory block of size - ~size_max~ \ge ~point_num~. - - #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -qmckl_exit_code -qmckl_get_point_xyz (const qmckl_context context, - double* const coord_x, - double* const coord_y, - double* const coord_z, - const int64_t size_max) -{ - - if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { - return QMCKL_INVALID_CONTEXT; - } - - if (coord_x == NULL) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_2, - "qmckl_get_point_coord_xyz", - "coord_x is a null pointer"); - } - - if (coord_y == NULL) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_3, - "qmckl_get_point_coord_xyz", - "coord_y is a null pointer"); - } - - if (coord_z == NULL) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_4, - "qmckl_get_point_coord_xyz", - "coord_z is a null pointer"); - } - - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; - assert (ctx != NULL); - assert (ctx->point != NULL); - - int64_t point_num = ctx->point->num; - - assert (ctx->point->coord_x != NULL); - assert (ctx->point->coord_y != NULL); - assert (ctx->point->coord_z != NULL); - - if (size_max < point_num) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_5, - "qmckl_get_point_coord_xyz", - "size_max too small"); - } - - memcpy(coord_x, ctx->point->coord_x, point_num*sizeof(double)); - memcpy(coord_y, ctx->point->coord_y, point_num*sizeof(double)); - memcpy(coord_z, ctx->point->coord_z, point_num*sizeof(double)); - - return QMCKL_SUCCESS; -} - #+end_src - - #+begin_src f90 :comments org :tangle (eval fh_func) :noweb yes -interface - integer(c_int32_t) function qmckl_get_point_xyz(context, & - coord_x, coord_y, coord_z, size_max) bind(C) - use, intrinsic :: iso_c_binding - import - implicit none - - integer (c_int64_t) , intent(in) , value :: context - real (c_double ) , intent(out) :: coord_x(*) - real (c_double ) , intent(out) :: coord_y(*) - real (c_double ) , intent(out) :: coord_z(*) - integer (c_int64_t) , intent(in) :: size_max - end function -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. - - #+NAME: check_alloc - #+begin_src c :exports none - if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { - return QMCKL_NULL_CONTEXT; - } - - qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; - assert (ctx != NULL); - assert (ctx->point != NULL); - - if (ctx->point->num < num) { - - if (ctx->point->coord_x != NULL) { - qmckl_free(context, ctx->point->coord_x); - ctx->point->coord_x = NULL; - } - - if (ctx->point->coord_y != NULL) { - qmckl_free(context, ctx->point->coord_y); - ctx->point->coord_y = NULL; - } - - if (ctx->point->coord_z != NULL) { - qmckl_free(context, ctx->point->coord_z); - ctx->point->coord_z = NULL; - } - - qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; - mem_info.size = num*sizeof(double); - - ctx->point->coord_x = (double*) qmckl_malloc(context, mem_info); - if (ctx->point->coord_x == NULL) { - return qmckl_failwith( context, - QMCKL_ALLOCATION_FAILED, - "qmckl_set_point", - NULL); - } - - ctx->point->coord_y = (double*) qmckl_malloc(context, mem_info); - if (ctx->point->coord_y == NULL) { - return qmckl_failwith( context, - QMCKL_ALLOCATION_FAILED, - "qmckl_set_point", - NULL); - } - - ctx->point->coord_z = (double*) qmckl_malloc(context, mem_info); - if (ctx->point->coord_z == NULL) { - return qmckl_failwith( context, - QMCKL_ALLOCATION_FAILED, - "qmckl_set_point", - NULL); - } - }; - - ctx->point->num = num; - #+end_src 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 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 qmckl_set_point (qmckl_context context, + const char transp, const double* coord, const int64_t num) { - - <> - for (int64_t i=0 ; ipoint->coord_x[i] = coord[3*i ]; - ctx->point->coord_y[i] = coord[3*i+1]; - ctx->point->coord_z[i] = coord[3*i+2]; + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; } + if (transp != 'N' && transp != 'T') { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "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); + + qmckl_exit_code rc; + if (ctx->point.num < num) { + + 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) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_set_point", + NULL); + } + + }; + + ctx->point.num = num; + + if (transp == 'T') { + 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]; + } + } + + /* Increment the date of the context */ + ctx->date += 1UL; + ctx->point.date = ctx->date; + return QMCKL_SUCCESS; } @@ -439,76 +339,38 @@ qmckl_set_point (qmckl_context context, #+begin_src f90 :comments org :tangle (eval fh_func) :noweb yes interface integer(c_int32_t) function qmckl_set_point(context, & - coord_x, coord_y, coord_z, size_max) bind(C) + transp, coord, size_max) bind(C) use, intrinsic :: iso_c_binding import implicit none integer (c_int64_t) , intent(in) , value :: context - real (c_double ) , intent(in) :: coord_x(*) - real (c_double ) , intent(in) :: coord_y(*) - real (c_double ) , intent(in) :: coord_z(*) + character(c_char) , intent(in) , value :: transp + real (c_double ) , intent(in) :: coord(*) integer (c_int64_t) , intent(in) , value :: size_max end function end interface #+end_src - #+begin_src c :comments org :tangle (eval h_func) -qmckl_exit_code qmckl_set_point_xyz (qmckl_context context, - const double* coord_x, - const double* coord_y, - const double* coord_z, - const int64_t num); - #+end_src - - #+begin_src c :comments org :tangle (eval c) :noweb yes -qmckl_exit_code -qmckl_set_point_xyz (qmckl_context context, - const double* coord_x, - const double* coord_y, - const double* coord_z, - const int64_t num) -{ - - <> - - memcpy(ctx->point->coord_x, coord_x, num*sizeof(double)); - memcpy(ctx->point->coord_y, coord_y, num*sizeof(double)); - memcpy(ctx->point->coord_z, coord_z, num*sizeof(double)); - - return QMCKL_SUCCESS; -} - #+end_src - - #+begin_src f90 :comments org :tangle (eval fh_func) :noweb yes -interface - integer(c_int32_t) function qmckl_set_point_xyz(context, & - coord_x, coord_y, coord_z, size_max) bind(C) - use, intrinsic :: iso_c_binding - import - implicit none - - integer (c_int64_t) , intent(in) , value :: context - real (c_double ) , intent(in) :: coord_x(*) - real (c_double ) , intent(in) :: coord_y(*) - real (c_double ) , intent(in) :: coord_z(*) - integer (c_int64_t) , intent(in) , value :: size_max - end function -end interface - #+end_src ** Test #+begin_src c :tangle (eval c_test) /* Reference input data */ int64_t point_num = chbrclf_elec_num; -double* coord = &(chbrclf_elec_coord[0][0][0]); +const double* coord = &(chbrclf_elec_coord[0][0][0]); /* --- */ qmckl_exit_code rc; +double coord2[point_num*3]; +double coord3[point_num*3]; -rc = qmckl_set_point (context, coord, point_num); + +rc = qmckl_get_point (context, 'N', coord2, (point_num*3)); +assert(rc == QMCKL_NOT_PROVIDED); + +rc = qmckl_set_point (context, 'N', coord, point_num); assert(rc == QMCKL_SUCCESS); int64_t n; @@ -516,25 +378,30 @@ rc = qmckl_get_point_num (context, &n); assert(rc == QMCKL_SUCCESS); assert(n == point_num); -double coord2[point_num*3]; -double coord_x[point_num]; -double coord_y[point_num]; -double coord_z[point_num]; - -rc = qmckl_get_point_xyz (context, coord_x, coord_y, coord_z, point_num); -assert(rc == QMCKL_SUCCESS); - -rc = qmckl_get_point (context, coord2, (point_num*3)); +rc = qmckl_get_point (context, 'N', coord2, (point_num*3)); assert(rc == QMCKL_SUCCESS); for (int64_t i=0 ; i<3*point_num ; ++i) { assert( coord[i] == coord2[i] ); } +rc = qmckl_get_point (context, 'T', coord2, (point_num*3)); +assert(rc == QMCKL_SUCCESS); + for (int64_t i=0 ; i 0~ - - ~n > 0~ - - ~lda >= m~ - - ~ldb >= n~ - - ~A~ is allocated with at least $m \times n \times 8$ bytes - - ~B~ is allocated with at least $n \times m \times 8$ bytes - -*** C header - - #+CALL: generate_c_header(table=qmckl_transpose_args,rettyp="qmckl_exit_code",fname="qmckl_transpose") - - #+RESULTS: - #+begin_src c :tangle (eval h_func) :comments org - qmckl_exit_code qmckl_transpose ( - const qmckl_context context, - const int64_t m, - const int64_t n, - const double* A, - const int64_t lda, - double* const B, - const int64_t ldb ); - #+end_src - -*** Source - #+begin_src f90 :tangle (eval f) -integer function qmckl_transpose_f(context, m, n, A, LDA, B, LDB) & - result(info) - use qmckl - implicit none - integer(qmckl_context) , intent(in) :: context - integer*8 , intent(in) :: m, n - integer*8 , intent(in) :: lda - real*8 , intent(in) :: A(lda,*) - integer*8 , intent(in) :: ldb - real*8 , intent(out) :: B(ldb,*) - - integer*8 :: i,j - - info = QMCKL_SUCCESS - - if (context == QMCKL_NULL_CONTEXT) then - info = QMCKL_INVALID_CONTEXT - return - endif - - if (m <= 0_8) then - info = QMCKL_INVALID_ARG_2 - return - endif - - if (n <= 0_8) then - info = QMCKL_INVALID_ARG_3 - return - endif - - if (LDA < m) then - info = QMCKL_INVALID_ARG_5 - return - endif - - if (LDB < n) then - info = QMCKL_INVALID_ARG_7 - return - endif - - do j=1,m - do i=1,n - B(i,j) = A(j,i) - end do - end do - -end function qmckl_transpose_f - #+end_src - -*** C interface :noexport: - - #+CALL: generate_c_interface(table=qmckl_transpose_args,rettyp="qmckl_exit_code",fname="qmckl_transpose") - - #+RESULTS: - #+begin_src f90 :tangle (eval f) :comments org :exports none - integer(c_int32_t) function qmckl_transpose & - (context, m, n, A, lda, B, ldb) & - bind(C) result(info) - - use, intrinsic :: iso_c_binding - implicit none - - integer (c_int64_t) , intent(in) , value :: context - integer (c_int64_t) , intent(in) , value :: m - integer (c_int64_t) , intent(in) , value :: n - real (c_double ) , intent(in) :: A(lda,*) - integer (c_int64_t) , intent(in) , value :: lda - real (c_double ) , intent(out) :: B(ldb,*) - integer (c_int64_t) , intent(in) , value :: ldb - - integer(c_int32_t), external :: qmckl_transpose_f - info = qmckl_transpose_f & - (context, m, n, A, lda, B, ldb) - - end function qmckl_transpose - #+end_src - - #+CALL: generate_f_interface(table=qmckl_transpose_args,rettyp="qmckl_exit_code",fname="qmckl_transpose") - - #+RESULTS: - #+begin_src f90 :tangle (eval fh_func) :comments org :exports none - interface - integer(c_int32_t) function qmckl_transpose & - (context, m, n, A, lda, B, ldb) & - bind(C) - use, intrinsic :: iso_c_binding - import - implicit none - - integer (c_int64_t) , intent(in) , value :: context - integer (c_int64_t) , intent(in) , value :: m - integer (c_int64_t) , intent(in) , value :: n - real (c_double ) , intent(in) :: A(lda,*) - integer (c_int64_t) , intent(in) , value :: lda - real (c_double ) , intent(out) :: B(ldb,*) - integer (c_int64_t) , intent(in) , value :: ldb - - end function qmckl_transpose - end interface - #+end_src - -*** Test :noexport: - #+begin_src f90 :tangle (eval f_test) -integer(qmckl_exit_code) function test_qmckl_transpose(context) bind(C) - use qmckl - implicit none - integer(qmckl_context), intent(in), value :: context - - double precision, allocatable :: A(:,:), B(:,:) - integer*8 :: m, n, LDA, LDB - integer*8 :: i,j - double precision :: x - - m = 5 - n = 6 - LDA = m+3 - LDB = n+1 - - allocate( A(LDA,n), B(LDB,m) ) - - A = 0.d0 - B = 0.d0 - do j=1,n - do i=1,m - A(i,j) = -10.d0 + dble(i+j) - end do - end do - - test_qmckl_transpose = qmckl_transpose(context, m, n, A, LDA, B, LDB) - - if (test_qmckl_transpose /= QMCKL_SUCCESS) return - - test_qmckl_transpose = QMCKL_FAILURE - - x = 0.d0 - do j=1,n - do i=1,m - x = x + (A(i,j)-B(j,i))**2 - end do - end do - if (dabs(x) <= 1.d-15) then - test_qmckl_transpose = QMCKL_SUCCESS - endif - - deallocate(A,B) -end function test_qmckl_transpose - #+end_src - - #+begin_src c :comments link :tangle (eval c_test) -qmckl_exit_code test_qmckl_transpose(qmckl_context context); -assert(QMCKL_SUCCESS == test_qmckl_transpose(context)); - #+end_src - -* End of files :noexport: - - #+begin_src c :comments link :tangle (eval c_test) - assert (qmckl_context_destroy(context) == QMCKL_SUCCESS); - return 0; -} - - #+end_src - - -# -*- mode: org -*- -# vim: syntax=c