TREX/qmckl
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-08-14 16:58:38 +02:00
Code Issues Projects Releases Wiki Activity

Introduced points in electron

Browse Source
This commit is contained in:
Anthony Scemama 2022-01-26 17:06:51 +01:00
parent 4b36005ca0
commit 0c9a50a681
5 changed files with 373 additions and 423 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

310
org/qmckl_ao.org
View File

@ -279,11 +279,11 @@ typedef struct qmckl_ao_basis_struct {
int32_t * nucleus_max_ang_mom; int32_t * nucleus_max_ang_mom;
double * nucleus_range; double * nucleus_range;
double * primitive_vgl; double * primitive_vgl;
int64_t primitive_vgl_date; uint64_t primitive_vgl_date;
double * shell_vgl; double * shell_vgl;
int64_t shell_vgl_date; uint64_t shell_vgl_date;
double * ao_vgl; double * ao_vgl;
int64_t ao_vgl_date; uint64_t ao_vgl_date;
int32_t uninitialized; int32_t uninitialized;
bool provided; bool provided;
@ -2457,15 +2457,15 @@ for (int64_t i=0 ; i < ao_num ; ++i) {
The following data is computed as described in the next sections: The following data is computed as described in the next sections:
|----------------------+----------------------------------+-----------------------------------------------------------------------------------------------| |----------------------+-----------------------------------+----------------------------------------------------------------------------------------------|
| Variable | Type | Description | | Variable | Type | Description |
|----------------------+----------------------------------+-----------------------------------------------------------------------------------------------| |----------------------+-----------------------------------+----------------------------------------------------------------------------------------------|
| ~primitive_vgl~ | ~double[5][elec_num][prim_num]~ | Value, gradients, Laplacian of the primitives at electron positions | | ~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 electron positions | | ~primitive_vgl_date~ | ~uint64_t~ | Last modification date of Value, gradients, Laplacian of the primitives at current positions |
| ~shell_vgl~ | ~double[5][elec_num][shell_num]~ | Value, gradients, Laplacian of the primitives at electron 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 electron positions | | ~shell_vgl_date~ | ~uint64_t~ | Last modification date of Value, gradients, Laplacian of the AOs at current positions |
| ~ao_vgl~ | ~double[5][elec_num][ao_num]~ | Value, gradients, Laplacian of the primitives at electron 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 electron positions | | ~ao_vgl_date~ | ~uint64_t~ | Last modification date of Value, gradients, Laplacian of the AOs at current positions |
*** After initialization *** 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; qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL); 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) { if (size_max < sze) {
return qmckl_failwith( context, return qmckl_failwith( context,
QMCKL_INVALID_ARG_3, 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; qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL); 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) { if (size_max < sze) {
return qmckl_failwith( context, return qmckl_failwith( context,
QMCKL_INVALID_ARG_3, 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; qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL); 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) { if (size_max < sze) {
return qmckl_failwith( context, return qmckl_failwith( context,
QMCKL_INVALID_ARG_3, QMCKL_INVALID_ARG_3,
@ -3046,16 +3046,16 @@ assert(0 == test_qmckl_ao_gaussian_vgl(context));
#+NAME: qmckl_ao_basis_primitive_gaussian_vgl_args #+NAME: qmckl_ao_basis_primitive_gaussian_vgl_args
| Variable | Type | In/Out | Description | | Variable | Type | In/Out | Description |
|----------------------+---------------------------------+--------+--------------------------------------------------| |----------------------+----------------------------------+--------+--------------------------------------------------|
| ~context~ | ~qmckl_context~ | in | Global state | | ~context~ | ~qmckl_context~ | in | Global state |
| ~prim_num~ | ~int64_t~ | in | Number of primitives | | ~prim_num~ | ~int64_t~ | in | Number of primitives |
| ~elec_num~ | ~int64_t~ | in | Number of electrons | | ~point_num~ | ~int64_t~ | in | Number of points considered |
| ~nucl_num~ | ~int64_t~ | in | Number of nuclei | | ~nucl_num~ | ~int64_t~ | in | Number of nuclei |
| ~nucleus_prim_index~ | ~int64_t[nucl_num]~ | in | Index of the 1st primitive of each nucleus | | ~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 | | ~coord~ | ~double[3][point_num]~ | in | Coordinates |
| ~nucl_coord~ | ~double[3][elec_num]~ | in | Nuclear coordinates | | ~nucl_coord~ | ~double[3][nucl_num]~ | in | Nuclear coordinates |
| ~expo~ | ~double[prim_num]~ | in | Exponents of the primitives | | ~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 | | ~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")) #+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,10 +3064,10 @@ assert(0 == test_qmckl_ao_gaussian_vgl(context));
qmckl_exit_code qmckl_compute_ao_basis_primitive_gaussian_vgl ( qmckl_exit_code qmckl_compute_ao_basis_primitive_gaussian_vgl (
const qmckl_context context, const qmckl_context context,
const int64_t prim_num, const int64_t prim_num,
const int64_t elec_num, const int64_t point_num,
const int64_t nucl_num, const int64_t nucl_num,
const int64_t* nucleus_prim_index, const int64_t* nucleus_prim_index,
const double* elec_coord, const double* coord,
const double* nucl_coord, const double* nucl_coord,
const double* expo, const double* expo,
double* const primitive_vgl ); double* const primitive_vgl );
@ -3076,8 +3076,8 @@ assert(0 == test_qmckl_ao_gaussian_vgl(context));
#+begin_src f90 :comments org :tangle (eval f) :noweb yes #+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_ao_basis_primitive_gaussian_vgl_f( & integer function qmckl_compute_ao_basis_primitive_gaussian_vgl_f( &
context, prim_num, elec_num, nucl_num, & context, prim_num, point_num, nucl_num, &
nucleus_prim_index, elec_coord, nucl_coord, & nucleus_prim_index, coord, nucl_coord, &
expo, primitive_vgl) & expo, primitive_vgl) &
result(info) result(info)
@ -3086,14 +3086,14 @@ integer function qmckl_compute_ao_basis_primitive_gaussian_vgl_f( &
integer(qmckl_context), intent(in) :: context integer(qmckl_context), intent(in) :: context
integer*8 , intent(in) :: prim_num integer*8 , intent(in) :: prim_num
integer*8 , intent(in) :: nucl_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) 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) :: nucl_coord(nucl_num,3)
double precision , intent(in) :: expo(prim_num) 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 double precision :: x, y, z, two_a, ar2, r2, v, cutoff
info = QMCKL_SUCCESS info = QMCKL_SUCCESS
@ -3104,10 +3104,10 @@ integer function qmckl_compute_ao_basis_primitive_gaussian_vgl_f( &
do inucl=1,nucl_num do inucl=1,nucl_num
! C is zero-based, so shift bounds by one ! C is zero-based, so shift bounds by one
do iprim = nucleus_prim_index(inucl)+1, nucleus_prim_index(inucl+1) do iprim = nucleus_prim_index(inucl)+1, nucleus_prim_index(inucl+1)
do ielec = 1, elec_num do ipoint = 1, point_num
x = elec_coord(ielec,1) - nucl_coord(inucl,1) x = coord(ipoint,1) - nucl_coord(inucl,1)
y = elec_coord(ielec,2) - nucl_coord(inucl,2) y = coord(ipoint,2) - nucl_coord(inucl,2)
z = elec_coord(ielec,3) - nucl_coord(inucl,3) z = coord(ipoint,3) - nucl_coord(inucl,3)
r2 = x*x + y*y + z*z r2 = x*x + y*y + z*z
ar2 = expo(iprim)*r2 ar2 = expo(iprim)*r2
@ -3116,11 +3116,11 @@ integer function qmckl_compute_ao_basis_primitive_gaussian_vgl_f( &
v = dexp(-ar2) v = dexp(-ar2)
two_a = -2.d0 * expo(iprim) * v two_a = -2.d0 * expo(iprim) * v
primitive_vgl(iprim, ielec, 1) = v primitive_vgl(iprim, ipoint, 1) = v
primitive_vgl(iprim, ielec, 2) = two_a * x primitive_vgl(iprim, ipoint, 2) = two_a * x
primitive_vgl(iprim, ielec, 3) = two_a * y primitive_vgl(iprim, ipoint, 3) = two_a * y
primitive_vgl(iprim, ielec, 4) = two_a * z primitive_vgl(iprim, ipoint, 4) = two_a * z
primitive_vgl(iprim, ielec, 5) = two_a * (3.d0 - 2.d0*ar2) primitive_vgl(iprim, ipoint, 5) = two_a * (3.d0 - 2.d0*ar2)
end do end do
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 & integer(c_int32_t) function qmckl_compute_ao_basis_primitive_gaussian_vgl &
(context, & (context, &
prim_num, & prim_num, &
elec_num, & point_num, &
nucl_num, & nucl_num, &
nucleus_prim_index, & nucleus_prim_index, &
elec_coord, & coord, &
nucl_coord, & nucl_coord, &
expo, & expo, &
primitive_vgl) & 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 :: context
integer (c_int64_t) , intent(in) , value :: prim_num 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) , value :: nucl_num
integer (c_int64_t) , intent(in) :: nucleus_prim_index(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) :: coord(point_num,3)
real (c_double ) , intent(in) :: nucl_coord(elec_num,3) real (c_double ) , intent(in) :: nucl_coord(nucl_num,3)
real (c_double ) , intent(in) :: expo(prim_num) 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 integer(c_int32_t), external :: qmckl_compute_ao_basis_primitive_gaussian_vgl_f
info = qmckl_compute_ao_basis_primitive_gaussian_vgl_f & info = qmckl_compute_ao_basis_primitive_gaussian_vgl_f &
(context, & (context, &
prim_num, & prim_num, &
elec_num, & point_num, &
nucl_num, & nucl_num, &
nucleus_prim_index, & nucleus_prim_index, &
elec_coord, & coord, &
nucl_coord, & nucl_coord, &
expo, & expo, &
primitive_vgl) primitive_vgl)
@ -3201,13 +3201,13 @@ qmckl_exit_code qmckl_provide_ao_basis_primitive_vgl(qmckl_context context)
} }
/* Compute if necessary */ /* 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 */ /* Allocate array */
if (ctx->ao_basis.primitive_vgl == NULL) { if (ctx->ao_basis.primitive_vgl == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; 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); sizeof(double);
double* primitive_vgl = (double*) qmckl_malloc(context, mem_info); double* primitive_vgl = (double*) qmckl_malloc(context, mem_info);
@ -3224,10 +3224,10 @@ qmckl_exit_code qmckl_provide_ao_basis_primitive_vgl(qmckl_context context)
if (ctx->ao_basis.type == 'G') { if (ctx->ao_basis.type == 'G') {
rc = qmckl_compute_ao_basis_primitive_gaussian_vgl(context, rc = qmckl_compute_ao_basis_primitive_gaussian_vgl(context,
ctx->ao_basis.prim_num, ctx->ao_basis.prim_num,
ctx->electron.num, ctx->point.num,
ctx->nucleus.num, ctx->nucleus.num,
ctx->ao_basis.nucleus_prim_index, ctx->ao_basis.nucleus_prim_index,
ctx->electron.coord_new, ctx->point.coord.data,
ctx->nucleus.coord.data, ctx->nucleus.coord.data,
ctx->ao_basis.exponent, ctx->ao_basis.exponent,
ctx->ao_basis.primitive_vgl); ctx->ao_basis.primitive_vgl);
@ -3312,15 +3312,15 @@ print ( "[7][4][26] : %e"% lf(a,x,y))
assert(qmckl_electron_provided(context)); 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); 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]), 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 (rc == QMCKL_SUCCESS);
assert( fabs(prim_vgl[0][26][7] - ( 1.0501570432064878E-003)) < 1.e-14 ); 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 // l : primitives
k=0; k=0;
for (j=0 ; j<elec_num ; ++j) { for (j=0 ; j<point_num ; ++j) {
for (i=0 ; i<nucl_num ; ++i) { for (i=0 ; i<nucl_num ; ++i) {
r2 = nucl_elec_dist[i][j]; r2 = nucl_point_dist[i][j];
if (r2 < nucl_radius2[i]) { if (r2 < nucl_radius2[i]) {
@ -3373,21 +3373,21 @@ for (j=0 ; j<elec_num ; ++j) {
#+NAME: qmckl_ao_basis_shell_gaussian_vgl_args #+NAME: qmckl_ao_basis_shell_gaussian_vgl_args
| Variable | Type | In/Out | Description | | Variable | Type | In/Out | Description |
|---------------------+----------------------------------+--------+----------------------------------------------| |---------------------+-----------------------------------+--------+----------------------------------------------|
| ~context~ | ~qmckl_context~ | in | Global state | | ~context~ | ~qmckl_context~ | in | Global state |
| ~prim_num~ | ~int64_t~ | in | Number of primitives | | ~prim_num~ | ~int64_t~ | in | Number of primitives |
| ~shell_num~ | ~int64_t~ | in | Number of shells | | ~shell_num~ | ~int64_t~ | in | Number of shells |
| ~elec_num~ | ~int64_t~ | in | Number of electrons | | ~point_num~ | ~int64_t~ | in | Number of points |
| ~nucl_num~ | ~int64_t~ | in | Number of nuclei | | ~nucl_num~ | ~int64_t~ | in | Number of nuclei |
| ~nucleus_shell_num~ | ~int64_t[nucl_num]~ | in | Number of shells for each nucleus | | ~nucleus_shell_num~ | ~int64_t[nucl_num]~ | in | Number of shells for each nucleus |
| ~nucleus_index~ | ~int64_t[nucl_num]~ | in | Index of the 1st shell of each nucleus | | ~nucleus_index~ | ~int64_t[nucl_num]~ | in | Index of the 1st shell of each nucleus |
| ~shell_prim_index~ | ~int64_t[shell_num]~ | in | Index of the 1st primitive of each shell | | ~shell_prim_index~ | ~int64_t[shell_num]~ | in | Index of the 1st primitive of each shell |
| ~shell_prim_num~ | ~int64_t[shell_num]~ | in | Number of primitives per shell | | ~shell_prim_num~ | ~int64_t[shell_num]~ | in | Number of primitives per shell |
| ~elec_coord~ | ~double[3][elec_num]~ | in | Electron coordinates | | ~coord~ | ~double[3][point_num]~ | in | Coordinates |
| ~nucl_coord~ | ~double[3][elec_num]~ | in | Nuclear coordinates | | ~nucl_coord~ | ~double[3][nucl_num]~ | in | Nuclear coordinates |
| ~expo~ | ~double[prim_num]~ | in | Exponents of the primitives | | ~expo~ | ~double[prim_num]~ | in | Exponents of the primitives |
| ~coef_normalized~ | ~double[prim_num]~ | in | Coefficients of the primitives | | ~coef_normalized~ | ~double[prim_num]~ | in | Coefficients of the primitives |
| ~shell_vgl~ | ~double[5][elec_num][shell_num]~ | out | Value, gradients and Laplacian of the shells | | ~shell_vgl~ | ~double[5][point_num][shell_num]~ | out | Value, gradients and Laplacian of the shells |
#+CALL: generate_c_header(table=qmckl_ao_basis_shell_gaussian_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_ao_basis_shell_gaussian_vgl")) #+CALL: generate_c_header(table=qmckl_ao_basis_shell_gaussian_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_ao_basis_shell_gaussian_vgl"))
@ -3397,13 +3397,13 @@ for (j=0 ; j<elec_num ; ++j) {
const qmckl_context context, const qmckl_context context,
const int64_t prim_num, const int64_t prim_num,
const int64_t shell_num, const int64_t shell_num,
const int64_t elec_num, const int64_t point_num,
const int64_t nucl_num, const int64_t nucl_num,
const int64_t* nucleus_shell_num, const int64_t* nucleus_shell_num,
const int64_t* nucleus_index, const int64_t* nucleus_index,
const int64_t* shell_prim_index, const int64_t* shell_prim_index,
const int64_t* shell_prim_num, const int64_t* shell_prim_num,
const double* elec_coord, const double* coord,
const double* nucl_coord, const double* nucl_coord,
const double* expo, const double* expo,
const double* coef_normalized, const double* coef_normalized,
@ -3412,9 +3412,9 @@ for (j=0 ; j<elec_num ; ++j) {
#+begin_src f90 :comments org :tangle (eval f) :noweb yes #+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_ao_basis_shell_gaussian_vgl_f( & integer function qmckl_compute_ao_basis_shell_gaussian_vgl_f( &
context, prim_num, shell_num, elec_num, nucl_num, & context, prim_num, shell_num, point_num, nucl_num, &
nucleus_shell_num, nucleus_index, shell_prim_index, & nucleus_shell_num, nucleus_index, shell_prim_index, &
shell_prim_num, elec_coord, nucl_coord, expo, & shell_prim_num, coord, nucl_coord, expo, &
coef_normalized, shell_vgl) & coef_normalized, shell_vgl) &
result(info) result(info)
use qmckl use qmckl
@ -3423,18 +3423,18 @@ integer function qmckl_compute_ao_basis_shell_gaussian_vgl_f( &
integer*8 , intent(in) :: prim_num integer*8 , intent(in) :: prim_num
integer*8 , intent(in) :: shell_num integer*8 , intent(in) :: shell_num
integer*8 , intent(in) :: nucl_num integer*8 , intent(in) :: nucl_num
integer*8 , intent(in) :: elec_num integer*8 , intent(in) :: point_num
integer*8 , intent(in) :: nucleus_shell_num(nucl_num) integer*8 , intent(in) :: nucleus_shell_num(nucl_num)
integer*8 , intent(in) :: nucleus_index(nucl_num) integer*8 , intent(in) :: nucleus_index(nucl_num)
integer*8 , intent(in) :: shell_prim_index(shell_num) integer*8 , intent(in) :: shell_prim_index(shell_num)
integer*8 , intent(in) :: shell_prim_num(shell_num) integer*8 , intent(in) :: shell_prim_num(shell_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) double precision , intent(in) :: nucl_coord(nucl_num,3)
double precision , intent(in) :: expo(prim_num) double precision , intent(in) :: expo(prim_num)
double precision , intent(in) :: coef_normalized(prim_num) double precision , intent(in) :: coef_normalized(prim_num)
double precision , intent(out) :: shell_vgl(shell_num,elec_num,5) double precision , intent(out) :: shell_vgl(shell_num,point_num,5)
integer*8 :: inucl, iprim, ielec, ishell integer*8 :: inucl, iprim, ipoint, ishell
integer*8 :: ishell_start, ishell_end integer*8 :: ishell_start, ishell_end
integer*8 :: iprim_start , iprim_end integer*8 :: iprim_start , iprim_end
double precision :: x, y, z, two_a, ar2, r2, v, cutoff double precision :: x, y, z, two_a, ar2, r2, v, cutoff
@ -3451,21 +3451,21 @@ integer function qmckl_compute_ao_basis_shell_gaussian_vgl_f( &
ishell_start = nucleus_index(inucl) + 1 ishell_start = nucleus_index(inucl) + 1
ishell_end = nucleus_index(inucl) + nucleus_shell_num(inucl) ishell_end = nucleus_index(inucl) + nucleus_shell_num(inucl)
do ielec = 1, elec_num do ipoint = 1, point_num
x = elec_coord(ielec,1) - nucl_coord(inucl,1) x = coord(ipoint,1) - nucl_coord(inucl,1)
y = elec_coord(ielec,2) - nucl_coord(inucl,2) y = coord(ipoint,2) - nucl_coord(inucl,2)
z = elec_coord(ielec,3) - nucl_coord(inucl,3) z = coord(ipoint,3) - nucl_coord(inucl,3)
r2 = x*x + y*y + z*z r2 = x*x + y*y + z*z
do ishell=ishell_start, ishell_end do ishell=ishell_start, ishell_end
shell_vgl(ishell, ielec, 1) = 0.d0 shell_vgl(ishell, ipoint, 1) = 0.d0
shell_vgl(ishell, ielec, 2) = 0.d0 shell_vgl(ishell, ipoint, 2) = 0.d0
shell_vgl(ishell, ielec, 3) = 0.d0 shell_vgl(ishell, ipoint, 3) = 0.d0
shell_vgl(ishell, ielec, 4) = 0.d0 shell_vgl(ishell, ipoint, 4) = 0.d0
shell_vgl(ishell, ielec, 5) = 0.d0 shell_vgl(ishell, ipoint, 5) = 0.d0
iprim_start = shell_prim_index(ishell) + 1 iprim_start = shell_prim_index(ishell) + 1
iprim_end = shell_prim_index(ishell) + shell_prim_num(ishell) iprim_end = shell_prim_index(ishell) + shell_prim_num(ishell)
@ -3480,20 +3480,20 @@ integer function qmckl_compute_ao_basis_shell_gaussian_vgl_f( &
v = coef_normalized(iprim) * dexp(-ar2) v = coef_normalized(iprim) * dexp(-ar2)
two_a = -2.d0 * expo(iprim) * v two_a = -2.d0 * expo(iprim) * v
shell_vgl(ishell, ielec, 1) = & shell_vgl(ishell, ipoint, 1) = &
shell_vgl(ishell, ielec, 1) + v shell_vgl(ishell, ipoint, 1) + v
shell_vgl(ishell, ielec, 2) = & shell_vgl(ishell, ipoint, 2) = &
shell_vgl(ishell, ielec, 2) + two_a * x shell_vgl(ishell, ipoint, 2) + two_a * x
shell_vgl(ishell, ielec, 3) = & shell_vgl(ishell, ipoint, 3) = &
shell_vgl(ishell, ielec, 3) + two_a * y shell_vgl(ishell, ipoint, 3) + two_a * y
shell_vgl(ishell, ielec, 4) = & shell_vgl(ishell, ipoint, 4) = &
shell_vgl(ishell, ielec, 4) + two_a * z shell_vgl(ishell, ipoint, 4) + two_a * z
shell_vgl(ishell, ielec, 5) = & shell_vgl(ishell, ipoint, 5) = &
shell_vgl(ishell, ielec, 5) + two_a * (3.d0 - 2.d0*ar2) shell_vgl(ishell, ipoint, 5) + two_a * (3.d0 - 2.d0*ar2)
end do end do
@ -3513,13 +3513,13 @@ end function qmckl_compute_ao_basis_shell_gaussian_vgl_f
(context, & (context, &
prim_num, & prim_num, &
shell_num, & shell_num, &
elec_num, & point_num, &
nucl_num, & nucl_num, &
nucleus_shell_num, & nucleus_shell_num, &
nucleus_index, & nucleus_index, &
shell_prim_index, & shell_prim_index, &
shell_prim_num, & shell_prim_num, &
elec_coord, & coord, &
nucl_coord, & nucl_coord, &
expo, & expo, &
coef_normalized, & coef_normalized, &
@ -3532,30 +3532,30 @@ end function qmckl_compute_ao_basis_shell_gaussian_vgl_f
integer (c_int64_t) , intent(in) , value :: context integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(in) , value :: prim_num integer (c_int64_t) , intent(in) , value :: prim_num
integer (c_int64_t) , intent(in) , value :: shell_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 integer (c_int64_t) , intent(in) , value :: nucl_num
integer (c_int64_t) , intent(in) :: nucleus_shell_num(nucl_num) integer (c_int64_t) , intent(in) :: nucleus_shell_num(nucl_num)
integer (c_int64_t) , intent(in) :: nucleus_index(nucl_num) integer (c_int64_t) , intent(in) :: nucleus_index(nucl_num)
integer (c_int64_t) , intent(in) :: shell_prim_index(shell_num) integer (c_int64_t) , intent(in) :: shell_prim_index(shell_num)
integer (c_int64_t) , intent(in) :: shell_prim_num(shell_num) integer (c_int64_t) , intent(in) :: shell_prim_num(shell_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(elec_num,3) real (c_double ) , intent(in) :: nucl_coord(nucl_num,3)
real (c_double ) , intent(in) :: expo(prim_num) real (c_double ) , intent(in) :: expo(prim_num)
real (c_double ) , intent(in) :: coef_normalized(prim_num) real (c_double ) , intent(in) :: coef_normalized(prim_num)
real (c_double ) , intent(out) :: shell_vgl(shell_num,elec_num,5) real (c_double ) , intent(out) :: shell_vgl(shell_num,point_num,5)
integer(c_int32_t), external :: qmckl_compute_ao_basis_shell_gaussian_vgl_f integer(c_int32_t), external :: qmckl_compute_ao_basis_shell_gaussian_vgl_f
info = qmckl_compute_ao_basis_shell_gaussian_vgl_f & info = qmckl_compute_ao_basis_shell_gaussian_vgl_f &
(context, & (context, &
prim_num, & prim_num, &
shell_num, & shell_num, &
elec_num, & point_num, &
nucl_num, & nucl_num, &
nucleus_shell_num, & nucleus_shell_num, &
nucleus_index, & nucleus_index, &
shell_prim_index, & shell_prim_index, &
shell_prim_num, & shell_prim_num, &
elec_coord, & coord, &
nucl_coord, & nucl_coord, &
expo, & expo, &
coef_normalized, & coef_normalized, &
@ -3591,21 +3591,14 @@ qmckl_exit_code qmckl_provide_ao_basis_shell_vgl(qmckl_context context)
NULL); NULL);
} }
if(!(ctx->electron.provided)) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_electron",
NULL);
}
/* Compute if necessary */ /* 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 */ /* Allocate array */
if (ctx->ao_basis.shell_vgl == NULL) { if (ctx->ao_basis.shell_vgl == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; 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); double* shell_vgl = (double*) qmckl_malloc(context, mem_info);
if (shell_vgl == NULL) { if (shell_vgl == NULL) {
@ -3622,13 +3615,13 @@ qmckl_exit_code qmckl_provide_ao_basis_shell_vgl(qmckl_context context)
rc = qmckl_compute_ao_basis_shell_gaussian_vgl(context, rc = qmckl_compute_ao_basis_shell_gaussian_vgl(context,
ctx->ao_basis.prim_num, ctx->ao_basis.prim_num,
ctx->ao_basis.shell_num, ctx->ao_basis.shell_num,
ctx->electron.num, ctx->point.num,
ctx->nucleus.num, ctx->nucleus.num,
ctx->ao_basis.nucleus_shell_num, ctx->ao_basis.nucleus_shell_num,
ctx->ao_basis.nucleus_index, ctx->ao_basis.nucleus_index,
ctx->ao_basis.shell_prim_index, ctx->ao_basis.shell_prim_index,
ctx->ao_basis.shell_prim_num, ctx->ao_basis.shell_prim_num,
ctx->electron.coord_new, ctx->point.coord.data,
ctx->nucleus.coord.data, ctx->nucleus.coord.data,
ctx->ao_basis.exponent, ctx->ao_basis.exponent,
ctx->ao_basis.coefficient_normalized, ctx->ao_basis.coefficient_normalized,
@ -3680,7 +3673,7 @@ elec_15_w2 = np.array( [ -2.20180344582,-1.9113150239, 2.2193744778600002 ] )
nucl_1 = np.array( [ 1.096243353458458e+00, 8.907054016973815e-01, 7.777092280258892e-01 ] ) nucl_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 ] ) 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 x = elec_26_w1 ; y = nucl_1
a = [( 8.236000E+03, -1.130000E-04 * 6.1616545431994848e+02 ), a = [( 8.236000E+03, -1.130000E-04 * 6.1616545431994848e+02 ),
( 1.235000E+03, -8.780000E-04 * 1.4847738511079908e+02 ), ( 1.235000E+03, -8.780000E-04 * 1.4847738511079908e+02 ),
@ -3726,7 +3719,7 @@ print ( "[1][4][26] : %25.15e"% lf(a,x,y))
assert(qmckl_electron_provided(context)); 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); assert(rc == QMCKL_SUCCESS);
@ -4341,13 +4334,13 @@ end function test_qmckl_ao_polynomial_vgl
#+NAME: qmckl_ao_vgl_args #+NAME: qmckl_ao_vgl_args
| Variable | Type | In/Out | Description | | Variable | Type | In/Out | Description |
|-----------------------+----------------------------------+--------+----------------------------------------------| |-----------------------+-----------------------------------+--------+----------------------------------------------|
| ~context~ | ~qmckl_context~ | in | Global state | | ~context~ | ~qmckl_context~ | in | Global state |
| ~ao_num~ | ~int64_t~ | in | Number of AOs | | ~ao_num~ | ~int64_t~ | in | Number of AOs |
| ~shell_num~ | ~int64_t~ | in | Number of shells | | ~shell_num~ | ~int64_t~ | in | Number of shells |
| ~elec_num~ | ~int64_t~ | in | Number of electrons | | ~point_num~ | ~int64_t~ | in | Number of points |
| ~nucl_num~ | ~int64_t~ | in | Number of nuclei | | ~nucl_num~ | ~int64_t~ | in | Number of nuclei |
| ~elec_coord~ | ~double[3][elec_num]~ | in | Electron coordinates | | ~coord~ | ~double[3][point_num]~ | in | Coordinates |
| ~nucl_coord~ | ~double[3][nucl_num]~ | in | Nuclear 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_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_shell_num~ | ~int64_t[nucl_num]~ | in | Number of shells per nucleus |
@ -4355,13 +4348,13 @@ end function test_qmckl_ao_polynomial_vgl
| ~nucleus_max_ang_mom~ | ~int32_t[nucl_num]~ | in | Maximum angular momentum per nucleus | | ~nucleus_max_ang_mom~ | ~int32_t[nucl_num]~ | in | Maximum angular momentum per nucleus |
| ~shell_ang_mom~ | ~int32_t[shell_num]~ | in | Angular momentum of each shell | | ~shell_ang_mom~ | ~int32_t[shell_num]~ | in | Angular momentum of each shell |
| ~ao_factor~ | ~double[ao_num]~ | in | Normalization factor of the AOs | | ~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 | | ~shell_vgl~ | ~double[5][point_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 | | ~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 #+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_ao_vgl_f(context, & integer function qmckl_compute_ao_vgl_f(context, &
ao_num, shell_num, elec_num, nucl_num, & ao_num, shell_num, point_num, nucl_num, &
elec_coord, nucl_coord, nucleus_index, nucleus_shell_num, & coord, nucl_coord, nucleus_index, nucleus_shell_num, &
nucleus_range, nucleus_max_ang_mom, shell_ang_mom, & nucleus_range, nucleus_max_ang_mom, shell_ang_mom, &
ao_factor, shell_vgl, ao_vgl) & ao_factor, shell_vgl, ao_vgl) &
result(info) result(info)
@ -4370,9 +4363,9 @@ integer function qmckl_compute_ao_vgl_f(context, &
integer(qmckl_context), intent(in) :: context integer(qmckl_context), intent(in) :: context
integer*8 , intent(in) :: ao_num integer*8 , intent(in) :: ao_num
integer*8 , intent(in) :: shell_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 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) double precision , intent(in) :: nucl_coord(nucl_num,3)
integer*8 , intent(in) :: nucleus_index(nucl_num) integer*8 , intent(in) :: nucleus_index(nucl_num)
integer*8 , intent(in) :: nucleus_shell_num(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) :: nucleus_max_ang_mom(nucl_num)
integer , intent(in) :: shell_ang_mom(shell_num) integer , intent(in) :: shell_ang_mom(shell_num)
double precision , intent(in) :: ao_factor(ao_num) double precision , intent(in) :: ao_factor(ao_num)
double precision , intent(in) :: shell_vgl(shell_num,elec_num,5) double precision , intent(in) :: shell_vgl(shell_num,point_num,5)
double precision , intent(out) :: ao_vgl(ao_num,elec_num,5) double precision , intent(out) :: ao_vgl(ao_num,point_num,5)
double precision :: e_coord(3), n_coord(3) double precision :: e_coord(3), n_coord(3)
integer*8 :: n_poly integer*8 :: n_poly
integer :: l, il, k integer :: l, il, k
integer*8 :: ielec, inucl, ishell integer*8 :: ipoint, inucl, ishell
integer*8 :: ishell_start, ishell_end integer*8 :: ishell_start, ishell_end
integer :: lstart(0:20) integer :: lstart(0:20)
double precision :: x, y, z, r2 double precision :: x, y, z, r2
@ -4409,17 +4402,17 @@ integer function qmckl_compute_ao_vgl_f(context, &
! TODO : Use numerical precision here ! TODO : Use numerical precision here
cutoff = -dlog(1.d-15) cutoff = -dlog(1.d-15)
do ielec = 1, elec_num do ipoint = 1, point_num
e_coord(1) = elec_coord(ielec,1) e_coord(1) = coord(ipoint,1)
e_coord(2) = elec_coord(ielec,2) e_coord(2) = coord(ipoint,2)
e_coord(3) = elec_coord(ielec,3) e_coord(3) = coord(ipoint,3)
k=1 k=1
do inucl=1,nucl_num do inucl=1,nucl_num
n_coord(1) = nucl_coord(inucl,1) n_coord(1) = nucl_coord(inucl,1)
n_coord(2) = nucl_coord(inucl,2) n_coord(2) = nucl_coord(inucl,2)
n_coord(3) = nucl_coord(inucl,3) 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) x = e_coord(1) - n_coord(1)
y = e_coord(2) - n_coord(2) y = e_coord(2) - n_coord(2)
z = e_coord(3) - n_coord(3) z = e_coord(3) - n_coord(3)
@ -4442,35 +4435,35 @@ integer function qmckl_compute_ao_vgl_f(context, &
l = shell_ang_mom(ishell) l = shell_ang_mom(ishell)
do il = lstart(l), lstart(l+1)-1 do il = lstart(l), lstart(l+1)-1
! Value ! Value
ao_vgl(k,ielec,1) = & ao_vgl(k,ipoint,1) = &
poly_vgl(1,il) * shell_vgl(ishell,ielec,1) * ao_factor(k) poly_vgl(1,il) * shell_vgl(ishell,ipoint,1) * ao_factor(k)
! Grad_x ! Grad_x
ao_vgl(k,ielec,2) = ( & ao_vgl(k,ipoint,2) = ( &
poly_vgl(2,il) * shell_vgl(ishell,ielec,1) + & poly_vgl(2,il) * shell_vgl(ishell,ipoint,1) + &
poly_vgl(1,il) * shell_vgl(ishell,ielec,2) & poly_vgl(1,il) * shell_vgl(ishell,ipoint,2) &
) * ao_factor(k) ) * ao_factor(k)
! Grad_y ! Grad_y
ao_vgl(k,ielec,3) = ( & ao_vgl(k,ipoint,3) = ( &
poly_vgl(3,il) * shell_vgl(ishell,ielec,1) + & poly_vgl(3,il) * shell_vgl(ishell,ipoint,1) + &
poly_vgl(1,il) * shell_vgl(ishell,ielec,3) & poly_vgl(1,il) * shell_vgl(ishell,ipoint,3) &
) * ao_factor(k) ) * ao_factor(k)
! Grad_z ! Grad_z
ao_vgl(k,ielec,4) = ( & ao_vgl(k,ipoint,4) = ( &
poly_vgl(4,il) * shell_vgl(ishell,ielec,1) + & poly_vgl(4,il) * shell_vgl(ishell,ipoint,1) + &
poly_vgl(1,il) * shell_vgl(ishell,ielec,4) & poly_vgl(1,il) * shell_vgl(ishell,ipoint,4) &
) * ao_factor(k) ) * ao_factor(k)
! Lapl_z ! Lapl_z
ao_vgl(k,ielec,5) = ( & ao_vgl(k,ipoint,5) = ( &
poly_vgl(5,il) * shell_vgl(ishell,ielec,1) + & poly_vgl(5,il) * shell_vgl(ishell,ipoint,1) + &
poly_vgl(1,il) * shell_vgl(ishell,ielec,5) + & poly_vgl(1,il) * shell_vgl(ishell,ipoint,5) + &
2.d0 * ( & 2.d0 * ( &
poly_vgl(2,il) * shell_vgl(ishell,ielec,2) + & poly_vgl(2,il) * shell_vgl(ishell,ipoint,2) + &
poly_vgl(3,il) * shell_vgl(ishell,ielec,3) + & poly_vgl(3,il) * shell_vgl(ishell,ipoint,3) + &
poly_vgl(4,il) * shell_vgl(ishell,ielec,4) ) & poly_vgl(4,il) * shell_vgl(ishell,ipoint,4) ) &
) * ao_factor(k) ) * ao_factor(k)
k = k+1 k = k+1
@ -4492,9 +4485,9 @@ end function qmckl_compute_ao_vgl_f
const qmckl_context context, const qmckl_context context,
const int64_t ao_num, const int64_t ao_num,
const int64_t shell_num, const int64_t shell_num,
const int64_t elec_num, const int64_t point_num,
const int64_t nucl_num, const int64_t nucl_num,
const double* elec_coord, const double* coord,
const double* nucl_coord, const double* nucl_coord,
const int64_t* nucleus_index, const int64_t* nucleus_index,
const int64_t* nucleus_shell_num, const int64_t* nucleus_shell_num,
@ -4514,9 +4507,9 @@ end function qmckl_compute_ao_vgl_f
(context, & (context, &
ao_num, & ao_num, &
shell_num, & shell_num, &
elec_num, & point_num, &
nucl_num, & nucl_num, &
elec_coord, & coord, &
nucl_coord, & nucl_coord, &
nucleus_index, & nucleus_index, &
nucleus_shell_num, & 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 :: context
integer (c_int64_t) , intent(in) , value :: ao_num 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 :: 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 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) 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_index(nucl_num)
integer (c_int64_t) , intent(in) :: nucleus_shell_num(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) :: nucleus_max_ang_mom(nucl_num)
integer (c_int32_t) , intent(in) :: shell_ang_mom(shell_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) :: ao_factor(ao_num)
real (c_double ) , intent(in) :: shell_vgl(shell_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,elec_num,5) real (c_double ) , intent(out) :: ao_vgl(ao_num,point_num,5)
integer(c_int32_t), external :: qmckl_compute_ao_vgl_f integer(c_int32_t), external :: qmckl_compute_ao_vgl_f
info = qmckl_compute_ao_vgl_f & info = qmckl_compute_ao_vgl_f &
(context, & (context, &
ao_num, & ao_num, &
shell_num, & shell_num, &
elec_num, & point_num, &
nucl_num, & nucl_num, &
elec_coord, & coord, &
nucl_coord, & nucl_coord, &
nucleus_index, & nucleus_index, &
nucleus_shell_num, & nucleus_shell_num, &
@ -4595,15 +4588,8 @@ qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context)
NULL); NULL);
} }
if(!(ctx->electron.provided)) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_electron",
NULL);
}
/* Compute if necessary */ /* 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; 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) { if (ctx->ao_basis.ao_vgl == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; 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); double* ao_vgl = (double*) qmckl_malloc(context, mem_info);
if (ao_vgl == NULL) { if (ao_vgl == NULL) {
@ -4632,9 +4618,9 @@ qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context)
rc = qmckl_compute_ao_vgl(context, rc = qmckl_compute_ao_vgl(context,
ctx->ao_basis.ao_num, ctx->ao_basis.ao_num,
ctx->ao_basis.shell_num, ctx->ao_basis.shell_num,
ctx->electron.num, ctx->point.num,
ctx->nucleus.num, ctx->nucleus.num,
ctx->electron.coord_new, ctx->point.coord.data,
ctx->nucleus.coord.data, ctx->nucleus.coord.data,
ctx->ao_basis.nucleus_index, ctx->ao_basis.nucleus_index,
ctx->ao_basis.nucleus_shell_num, ctx->ao_basis.nucleus_shell_num,
@ -4752,7 +4738,7 @@ assert (rc == QMCKL_SUCCESS);
assert(qmckl_electron_provided(context)); 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); assert(rc == QMCKL_SUCCESS);

2
org/qmckl_context.org
View File

@ -124,7 +124,7 @@ typedef struct qmckl_context_struct {
uint64_t date; uint64_t date;
/* Points */ /* Points */
qmckl_point_struct *point; qmckl_point_struct point;
/* -- Molecular system -- */ /* -- Molecular system -- */
qmckl_nucleus_struct nucleus; qmckl_nucleus_struct nucleus;

158
org/qmckl_electron.org
View File

@ -70,7 +70,7 @@ int main() {
The following data stored in the context: The following data stored in the context:
| Variable | Type | Description | | Variable | Type | Description |
|---------------------------+----------------------------+-------------------------------------------| |---------------------------+----------------+---------------------------------------------------------------|
| ~uninitialized~ | ~int32_t~ | Keeps bit set for uninitialized data | | ~uninitialized~ | ~int32_t~ | Keeps bit set for uninitialized data |
| ~num~ | ~int64_t~ | Total number of electrons | | ~num~ | ~int64_t~ | Total number of electrons |
| ~up_num~ | ~int64_t~ | Number of up-spin electrons | | ~up_num~ | ~int64_t~ | Number of up-spin electrons |
@ -79,8 +79,8 @@ int main() {
| ~rescale_factor_kappa_ee~ | ~double~ | The distance scaling factor | | ~rescale_factor_kappa_ee~ | ~double~ | The distance scaling factor |
| ~rescale_factor_kappa_en~ | ~double~ | The distance scaling factor | | ~rescale_factor_kappa_en~ | ~double~ | The distance scaling factor |
| ~provided~ | ~bool~ | If true, ~electron~ is valid | | ~provided~ | ~bool~ | If true, ~electron~ is valid |
| ~coord_new~ | ~double[3][walk_num][num]~ | New set of electron coordinates | | ~coord_new~ | ~qmckl_matrix~ | Current set of electron coordinates. Pointer to ~ctx->points~ |
| ~coord_old~ | ~double[3][walk_num][num]~ | Old set of electron coordinates | | ~coord_old~ | ~qmckl_matrix~ | Old set of electron coordinates |
| ~coord_new_date~ | ~uint64_t~ | Last modification date of the coordinates | | ~coord_new_date~ | ~uint64_t~ | Last modification date of the coordinates |
Computed data: Computed data:
@ -123,8 +123,8 @@ typedef struct qmckl_electron_struct {
int64_t ee_distance_rescaled_deriv_e_date; int64_t ee_distance_rescaled_deriv_e_date;
int64_t en_distance_rescaled_date; int64_t en_distance_rescaled_date;
int64_t en_distance_rescaled_deriv_e_date; int64_t en_distance_rescaled_deriv_e_date;
double* coord_new; qmckl_matrix coord_new;
double* coord_old; qmckl_matrix coord_old;
double* ee_distance; double* ee_distance;
double* en_distance; double* en_distance;
double* ee_pot; double* ee_pot;
@ -412,6 +412,10 @@ qmckl_get_electron_coord (const qmckl_context context,
const int64_t size_max); const int64_t size_max);
#+end_src #+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 #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_exit_code
qmckl_get_electron_coord (const qmckl_context context, qmckl_get_electron_coord (const qmckl_context context,
@ -419,10 +423,6 @@ qmckl_get_electron_coord (const qmckl_context context,
double* const coord, double* const coord,
const int64_t size_max) const int64_t size_max)
{ {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_INVALID_CONTEXT;
}
if (transp != 'N' && transp != 'T') { if (transp != 'N' && transp != 'T') {
return qmckl_failwith( context, return qmckl_failwith( context,
QMCKL_INVALID_ARG_2, QMCKL_INVALID_ARG_2,
@ -437,46 +437,32 @@ qmckl_get_electron_coord (const qmckl_context context,
"coord is a null pointer"); "coord is a null pointer");
} }
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; if (size_max <= 0) {
assert (ctx != NULL);
if ( !(ctx->electron.provided) ) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_get_electron_coord",
"electron data is not provided");
}
int64_t elec_num = ctx->electron.num;
int64_t walk_num = ctx->electron.walk_num;
if (size_max < elec_num*walk_num*3) {
return qmckl_failwith( context, return qmckl_failwith( context,
QMCKL_INVALID_ARG_4, QMCKL_INVALID_ARG_4,
"qmckl_get_electron_coord", "qmckl_get_electron_coord",
"Array too small"); "size_max should be > 0");
} }
assert (ctx->electron.coord_new != NULL); if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_INVALID_CONTEXT;
double* ptr1 = ctx->electron.coord_new;
double* ptr2 = coord;
if (transp == 'N') {
for (int64_t i=0 ; i<elec_num*walk_num ; ++i) {
coord[3*i] = ptr1[i] ;
coord[3*i+1] = ptr1[i+elec_num*walk_num] ;
coord[3*i+2] = ptr1[i+2*elec_num*walk_num] ;
} }
} else { qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
memcpy(ptr2, ptr1, 3*elec_num*walk_num*sizeof(double));
if (!ctx->electron.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_get_electron_coord",
NULL);
} }
return QMCKL_SUCCESS; assert (ctx->point.coord.data == ctx->electron.coord_new.data);
assert (ctx->point.coord.size[0] == ctx->electron.coord_new.size[0]);
assert (ctx->point.coord.size[1] == ctx->electron.coord_new.size[1]);
return qmckl_get_point(context, transp, coord, size_max);
} }
#+end_src #+end_src
@ -512,20 +498,23 @@ ctx->electron.uninitialized &= ~mask;
ctx->electron.provided = (ctx->electron.uninitialized == 0); ctx->electron.provided = (ctx->electron.uninitialized == 0);
if (ctx->electron.provided) { if (ctx->electron.provided) {
if (ctx->electron.coord_new != NULL) { if (ctx->electron.coord_new.data != NULL) {
qmckl_free(context, ctx->electron.coord_new); const qmckl_exit_code rc = qmckl_matrix_free(context, ctx->electron.coord_new);
ctx->electron.coord_new = NULL; assert (rc == QMCKL_SUCCESS);
} }
if (ctx->electron.coord_old != NULL) { if (ctx->electron.coord_old.data != NULL) {
qmckl_free(context, ctx->electron.coord_old); const qmckl_exit_code rc = qmckl_matrix_free(context, ctx->electron.coord_old);
ctx->electron.coord_old = NULL; assert (rc == QMCKL_SUCCESS);
} }
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; const int64_t walk_num = ctx->electron.walk_num;
mem_info.size = ctx->electron.num * ctx->electron.walk_num * 3 * sizeof(double); const int64_t elec_num = ctx->electron.num;
double* coord_new = (double*) qmckl_malloc(context, mem_info); const int64_t point_num = walk_num * elec_num;
if (coord_new == NULL) {
qmckl_matrix coord_new = qmckl_matrix_alloc(context, point_num, 3);
if (coord_new.data == NULL) {
return qmckl_failwith( context, return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED, QMCKL_ALLOCATION_FAILED,
"qmckl_set_electron_num", "qmckl_set_electron_num",
@ -533,8 +522,9 @@ if (ctx->electron.provided) {
} }
ctx->electron.coord_new = coord_new; ctx->electron.coord_new = coord_new;
double* coord_old = (double*) qmckl_malloc(context, mem_info); qmckl_matrix coord_old = qmckl_matrix_alloc(context, point_num, 3);
if (coord_old == NULL) {
if (coord_old.data == NULL) {
return qmckl_failwith( context, return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED, QMCKL_ALLOCATION_FAILED,
"qmckl_set_electron_num", "qmckl_set_electron_num",
@ -542,6 +532,8 @@ if (ctx->electron.provided) {
} }
ctx->electron.coord_old = coord_old; ctx->electron.coord_old = coord_old;
ctx->point.num = point_num;
ctx->point.coord = coord_new;
} }
return QMCKL_SUCCESS; return QMCKL_SUCCESS;
@ -672,6 +664,9 @@ end interface
overwritten. This can be done only when the data relative to overwritten. This can be done only when the data relative to
electrons have been set. 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 Important: changing the electron coordinates increments the date
in the context. in the context.
@ -699,10 +694,7 @@ qmckl_set_electron_coord(qmckl_context context,
"coord is a null pointer"); "coord is a null pointer");
} }
int64_t elec_num; const int64_t elec_num = ctx->electron.num;
qmckl_exit_code rc;
rc = qmckl_get_electron_num(context, &elec_num);
if (rc != QMCKL_SUCCESS) return rc;
if (elec_num == 0L) { if (elec_num == 0L) {
return qmckl_failwith( context, return qmckl_failwith( context,
@ -711,9 +703,9 @@ qmckl_set_electron_coord(qmckl_context context,
"elec_num is not set"); "elec_num is not set");
} }
int64_t walk_num; const int64_t walk_num = ctx->electron.walk_num;
rc = qmckl_get_electron_walk_num(context, &walk_num);
if (rc != QMCKL_SUCCESS) return rc; const int64_t point_num = ctx->point.num ;
if (walk_num == 0L) { if (walk_num == 0L) {
return qmckl_failwith( context, return qmckl_failwith( context,
@ -722,40 +714,18 @@ qmckl_set_electron_coord(qmckl_context context,
"walk_num is not set"); "walk_num is not set");
} }
if (size_max < walk_num*3*elec_num) { assert(point_num == elec_num * walk_num);
return qmckl_failwith( context,
QMCKL_INVALID_ARG_4,
"qmckl_set_electron_coord",
"Array too small");
}
/* If num and walk_num are set, the arrays should be allocated */
assert (ctx->electron.coord_old != NULL);
assert (ctx->electron.coord_new != NULL);
/* Increment the date of the context */
ctx->date += 1UL;
/* Swap pointers */ /* Swap pointers */
double * swap; ctx->point.coord = ctx->electron.coord_old;
swap = ctx->electron.coord_old; ctx->electron.coord_old = ctx->electron.coord_new ;
ctx->electron.coord_old = ctx->electron.coord_new;
ctx->electron.coord_new = swap;
double* ptr1 = ctx->electron.coord_new; qmckl_exit_code rc;
if (transp == 'N') { rc = qmckl_set_point(context, transp, coord, size_max/3);
assert (rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i<walk_num*elec_num ; ++i) { ctx->electron.coord_new = ctx->point.coord ;
ptr1[i] = coord[3*i];
ptr1[i+walk_num*elec_num] = coord[3*i+1];
ptr1[i+2*walk_num*elec_num] = coord[3*i+2];
}
} else {
memcpy(ptr1, coord, 3*elec_num*walk_num*sizeof(double));
}
ctx->electron.coord_new_date = ctx->date; ctx->electron.coord_new_date = ctx->date;
return QMCKL_SUCCESS; return QMCKL_SUCCESS;
@ -880,6 +850,7 @@ double elec_coord2[walk_num*3*elec_num];
rc = qmckl_get_electron_coord (context, 'N', elec_coord2, walk_num*3*elec_num); rc = qmckl_get_electron_coord (context, 'N', elec_coord2, walk_num*3*elec_num);
assert(rc == QMCKL_SUCCESS); assert(rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i<3*elec_num*walk_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] ); assert( elec_coord[i] == elec_coord2[i] );
} }
@ -981,7 +952,7 @@ qmckl_exit_code qmckl_provide_ee_distance(qmckl_context context)
qmckl_compute_ee_distance(context, qmckl_compute_ee_distance(context,
ctx->electron.num, ctx->electron.num,
ctx->electron.walk_num, ctx->electron.walk_num,
ctx->electron.coord_new, ctx->electron.coord_new.data,
ctx->electron.ee_distance); ctx->electron.ee_distance);
if (rc != QMCKL_SUCCESS) { if (rc != QMCKL_SUCCESS) {
return rc; return rc;
@ -1223,7 +1194,7 @@ qmckl_exit_code qmckl_provide_ee_distance_rescaled(qmckl_context context)
ctx->electron.num, ctx->electron.num,
ctx->electron.rescale_factor_kappa_en, ctx->electron.rescale_factor_kappa_en,
ctx->electron.walk_num, ctx->electron.walk_num,
ctx->electron.coord_new, ctx->electron.coord_new.data,
ctx->electron.ee_distance_rescaled); ctx->electron.ee_distance_rescaled);
if (rc != QMCKL_SUCCESS) { if (rc != QMCKL_SUCCESS) {
return rc; return rc;
@ -1469,7 +1440,7 @@ qmckl_exit_code qmckl_provide_ee_distance_rescaled_deriv_e(qmckl_context context
ctx->electron.num, ctx->electron.num,
ctx->electron.rescale_factor_kappa_en, ctx->electron.rescale_factor_kappa_en,
ctx->electron.walk_num, ctx->electron.walk_num,
ctx->electron.coord_new, ctx->electron.coord_new.data,
ctx->electron.ee_distance_rescaled_deriv_e); ctx->electron.ee_distance_rescaled_deriv_e);
if (rc != QMCKL_SUCCESS) { if (rc != QMCKL_SUCCESS) {
return rc; return rc;
@ -1910,7 +1881,7 @@ qmckl_exit_code qmckl_provide_en_distance(qmckl_context context)
ctx->electron.num, ctx->electron.num,
ctx->nucleus.num, ctx->nucleus.num,
ctx->electron.walk_num, ctx->electron.walk_num,
ctx->electron.coord_new, ctx->electron.coord_new.data,
ctx->nucleus.coord.data, ctx->nucleus.coord.data,
ctx->electron.en_distance); ctx->electron.en_distance);
if (rc != QMCKL_SUCCESS) { if (rc != QMCKL_SUCCESS) {
@ -2187,7 +2158,7 @@ qmckl_exit_code qmckl_provide_en_distance_rescaled(qmckl_context context)
ctx->nucleus.num, ctx->nucleus.num,
ctx->electron.rescale_factor_kappa_en, ctx->electron.rescale_factor_kappa_en,
ctx->electron.walk_num, ctx->electron.walk_num,
ctx->electron.coord_new, ctx->electron.coord_new.data,
ctx->nucleus.coord.data, ctx->nucleus.coord.data,
ctx->electron.en_distance_rescaled); ctx->electron.en_distance_rescaled);
if (rc != QMCKL_SUCCESS) { if (rc != QMCKL_SUCCESS) {
@ -2374,7 +2345,6 @@ assert (rc == QMCKL_SUCCESS);
assert(fabs(en_distance_rescaled[0][0][0] - 0.9994721712909764) < 1.e-12); assert(fabs(en_distance_rescaled[0][0][0] - 0.9994721712909764) < 1.e-12);
// (1,2,1) // (1,2,1)
printf("%f\n%f\n", en_distance_rescaled[0][1][0] , 0.9998448354439821);
assert(fabs(en_distance_rescaled[0][1][0] - 0.9998448354439821) < 1.e-12); assert(fabs(en_distance_rescaled[0][1][0] - 0.9998448354439821) < 1.e-12);
// (2,1,1) // (2,1,1)
@ -2476,7 +2446,7 @@ qmckl_exit_code qmckl_provide_en_distance_rescaled_deriv_e(qmckl_context context
ctx->nucleus.num, ctx->nucleus.num,
ctx->electron.rescale_factor_kappa_en, ctx->electron.rescale_factor_kappa_en,
ctx->electron.walk_num, ctx->electron.walk_num,
ctx->electron.coord_new, ctx->electron.coord_new.data,
ctx->nucleus.coord.data, ctx->nucleus.coord.data,
ctx->electron.en_distance_rescaled_deriv_e); ctx->electron.en_distance_rescaled_deriv_e);
if (rc != QMCKL_SUCCESS) { if (rc != QMCKL_SUCCESS) {

4
org/qmckl_jastrow.org
View File

@ -3260,7 +3260,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_e_deriv_e(qmckl_context context)
ctx->electron.num, ctx->electron.num,
ctx->jastrow.cord_num, ctx->jastrow.cord_num,
ctx->electron.rescale_factor_kappa_ee, ctx->electron.rescale_factor_kappa_ee,
ctx->electron.coord_new, ctx->electron.coord_new.data,
ctx->electron.ee_distance, ctx->electron.ee_distance,
ctx->jastrow.een_rescaled_e, ctx->jastrow.een_rescaled_e,
ctx->jastrow.een_rescaled_e_deriv_e); ctx->jastrow.een_rescaled_e_deriv_e);
@ -3924,7 +3924,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_n_deriv_e(qmckl_context context)
ctx->nucleus.num, ctx->nucleus.num,
ctx->jastrow.cord_num, ctx->jastrow.cord_num,
ctx->electron.rescale_factor_kappa_en, ctx->electron.rescale_factor_kappa_en,
ctx->electron.coord_new, ctx->electron.coord_new.data,
ctx->nucleus.coord.data, ctx->nucleus.coord.data,
ctx->electron.en_distance, ctx->electron.en_distance,
ctx->jastrow.een_rescaled_n, ctx->jastrow.een_rescaled_n,

60
org/qmckl_point.org
View File

@ -78,8 +78,9 @@ int main() {
The following data stored in the context: The following data stored in the context:
| Variable | Type | Description | | Variable | Type | Description |
|----------+----------------+------------------------| |----------+----------------+-------------------------------------------|
| ~num~ | ~int64_t~ | Total number of points | | ~num~ | ~int64_t~ | Total number of points |
| ~date~ | ~uint64_t~ | Last modification date of the coordinates |
| ~coord~ | ~qmckl_matrix~ | ~num~ \times 3 matrix3 | | ~coord~ | ~qmckl_matrix~ | ~num~ \times 3 matrix3 |
We consider that the matrix is stored 'transposed' and 'normal' We consider that the matrix is stored 'transposed' and 'normal'
@ -90,6 +91,7 @@ int main() {
#+begin_src c :comments org :tangle (eval h_private_type) #+begin_src c :comments org :tangle (eval h_private_type)
typedef struct qmckl_point_struct { typedef struct qmckl_point_struct {
int64_t num; int64_t num;
uint64_t date;
qmckl_matrix coord; qmckl_matrix coord;
} qmckl_point_struct; } qmckl_point_struct;
@ -109,16 +111,7 @@ qmckl_exit_code qmckl_init_point(qmckl_context context) {
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL); assert (ctx != NULL);
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; memset(&(ctx->point), 0, sizeof(qmckl_point_struct));
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));
return QMCKL_SUCCESS; return QMCKL_SUCCESS;
} }
@ -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; qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL); assert (ctx != NULL);
assert (ctx->point != NULL);
assert (ctx->point->num > (int64_t) 0); assert (ctx->point.num >= (int64_t) 0);
,*num = ctx->point->num; ,*num = ctx->point.num;
return QMCKL_SUCCESS; return QMCKL_SUCCESS;
} }
#+end_src #+end_src
@ -212,12 +204,11 @@ qmckl_get_point(const qmckl_context context,
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL); 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; if (point_num == 0) return QMCKL_NOT_PROVIDED;
assert (ctx->point->coord.data != NULL); assert (ctx->point.coord.data != NULL);
if (size_max < 3*point_num) { if (size_max < 3*point_num) {
return qmckl_failwith( context, return qmckl_failwith( context,
@ -229,13 +220,13 @@ qmckl_get_point(const qmckl_context context,
qmckl_exit_code rc; qmckl_exit_code rc;
if (transp == 'N') { if (transp == 'N') {
qmckl_matrix At = qmckl_matrix_alloc( context, 3, point_num); qmckl_matrix At = qmckl_matrix_alloc( context, 3, point_num);
rc = qmckl_transpose( context, ctx->point->coord, At ); rc = qmckl_transpose( context, ctx->point.coord, At );
if (rc != QMCKL_SUCCESS) return rc; if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_double_of_matrix( context, At, coord, size_max); rc = qmckl_double_of_matrix( context, At, coord, size_max);
if (rc != QMCKL_SUCCESS) return rc; if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_matrix_free(context, At); rc = qmckl_matrix_free(context, At);
} else { } else {
rc = qmckl_double_of_matrix( context, ctx->point->coord, coord, size_max); rc = qmckl_double_of_matrix( context, ctx->point.coord, coord, size_max);
} }
if (rc != QMCKL_SUCCESS) return rc; if (rc != QMCKL_SUCCESS) return rc;
@ -272,10 +263,10 @@ end interface
qmckl_exit_code qmckl_set_point (qmckl_context context, qmckl_exit_code qmckl_set_point (qmckl_context context,
const char transp, const char transp,
const double* coord, const double* coord,
const int64_t size_max); const int64_t num);
#+end_src #+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 #+begin_src c :comments org :tangle (eval c) :noweb yes
qmckl_exit_code qmckl_exit_code
@ -305,18 +296,17 @@ qmckl_set_point (qmckl_context context,
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL); assert (ctx != NULL);
assert (ctx->point != NULL);
qmckl_exit_code rc; qmckl_exit_code rc;
if (ctx->point->num < num) { if (ctx->point.num < num) {
if (ctx->point->coord.data != NULL) { if (ctx->point.coord.data != NULL) {
rc = qmckl_matrix_free(context, ctx->point->coord); rc = qmckl_matrix_free(context, ctx->point.coord);
return rc; assert (rc == QMCKL_SUCCESS);
} }
ctx->point->coord = qmckl_matrix_alloc(context, num, 3); ctx->point.coord = qmckl_matrix_alloc(context, num, 3);
if (ctx->point->coord.data == NULL) { if (ctx->point.coord.data == NULL) {
return qmckl_failwith( context, return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED, QMCKL_ALLOCATION_FAILED,
"qmckl_set_point", "qmckl_set_point",
@ -325,18 +315,22 @@ qmckl_set_point (qmckl_context context,
}; };
ctx->point->num = num; ctx->point.num = num;
if (transp == 'T') { if (transp == 'T') {
memcpy(ctx->point->coord.data, coord, 3*num*sizeof(double)); memcpy(ctx->point.coord.data, coord, 3*num*sizeof(double));
} else { } else {
for (int64_t i=0 ; i<num ; ++i) { for (int64_t i=0 ; i<num ; ++i) {
qmckl_mat(ctx->point->coord, i, 0) = coord[3*i ]; qmckl_mat(ctx->point.coord, i, 0) = coord[3*i ];
qmckl_mat(ctx->point->coord, i, 1) = coord[3*i+1]; qmckl_mat(ctx->point.coord, i, 1) = coord[3*i+1];
qmckl_mat(ctx->point->coord, i, 2) = coord[3*i+2]; 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; return QMCKL_SUCCESS;
} }