TREX/qmckl
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2025-01-06 19:33:14 +01: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

374
org/qmckl_ao.org
View File

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

2
org/qmckl_context.org
View File

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

322
org/qmckl_electron.org
View File

@ -69,22 +69,22 @@ int main() {
The following data stored in the context:
| Variable | Type | Description |
|---------------------------+----------------------------+-------------------------------------------|
| ~uninitialized~ | ~int32_t~ | Keeps bit set for uninitialized data |
| ~num~ | ~int64_t~ | Total number of electrons |
| ~up_num~ | ~int64_t~ | Number of up-spin electrons |
| ~down_num~ | ~int64_t~ | Number of down-spin electrons |
| ~walk_num~ | ~int64_t~ | Number of walkers |
| ~rescale_factor_kappa_ee~ | ~double~ | The distance scaling factor |
| ~rescale_factor_kappa_en~ | ~double~ | The distance scaling factor |
| ~provided~ | ~bool~ | If true, ~electron~ is valid |
| ~coord_new~ | ~double[3][walk_num][num]~ | New set of electron coordinates |
| ~coord_old~ | ~double[3][walk_num][num]~ | Old set of electron coordinates |
| ~coord_new_date~ | ~uint64_t~ | Last modification date of the coordinates |
| Variable | Type | Description |
|---------------------------+----------------+---------------------------------------------------------------|
| ~uninitialized~ | ~int32_t~ | Keeps bit set for uninitialized data |
| ~num~ | ~int64_t~ | Total number of electrons |
| ~up_num~ | ~int64_t~ | Number of up-spin electrons |
| ~down_num~ | ~int64_t~ | Number of down-spin electrons |
| ~walk_num~ | ~int64_t~ | Number of walkers |
| ~rescale_factor_kappa_ee~ | ~double~ | The distance scaling factor |
| ~rescale_factor_kappa_en~ | ~double~ | The distance scaling factor |
| ~provided~ | ~bool~ | If true, ~electron~ is valid |
| ~coord_new~ | ~qmckl_matrix~ | Current set of electron coordinates. Pointer to ~ctx->points~ |
| ~coord_old~ | ~qmckl_matrix~ | Old set of electron coordinates |
| ~coord_new_date~ | ~uint64_t~ | Last modification date of the coordinates |
Computed data:
| Variable | Type | Description |
|-------------------------------------+--------------------------------------+----------------------------------------------------------------------|
| ~ee_distance~ | ~double[walk_num][num][num]~ | Electron-electron distances |
@ -108,33 +108,33 @@ int main() {
#+begin_src c :comments org :tangle (eval h_private_type)
typedef struct qmckl_electron_struct {
int64_t num;
int64_t up_num;
int64_t down_num;
int64_t walk_num;
double rescale_factor_kappa_ee;
double rescale_factor_kappa_en;
int64_t coord_new_date;
int64_t ee_distance_date;
int64_t en_distance_date;
int64_t ee_pot_date;
int64_t en_pot_date;
int64_t ee_distance_rescaled_date;
int64_t ee_distance_rescaled_deriv_e_date;
int64_t en_distance_rescaled_date;
int64_t en_distance_rescaled_deriv_e_date;
double* coord_new;
double* coord_old;
double* ee_distance;
double* en_distance;
double* ee_pot;
double* en_pot;
double* ee_distance_rescaled;
double* ee_distance_rescaled_deriv_e;
double* en_distance_rescaled;
double* en_distance_rescaled_deriv_e;
int32_t uninitialized;
bool provided;
int64_t num;
int64_t up_num;
int64_t down_num;
int64_t walk_num;
double rescale_factor_kappa_ee;
double rescale_factor_kappa_en;
int64_t coord_new_date;
int64_t ee_distance_date;
int64_t en_distance_date;
int64_t ee_pot_date;
int64_t en_pot_date;
int64_t ee_distance_rescaled_date;
int64_t ee_distance_rescaled_deriv_e_date;
int64_t en_distance_rescaled_date;
int64_t en_distance_rescaled_deriv_e_date;
qmckl_matrix coord_new;
qmckl_matrix coord_old;
double* ee_distance;
double* en_distance;
double* ee_pot;
double* en_pot;
double* ee_distance_rescaled;
double* ee_distance_rescaled_deriv_e;
double* en_distance_rescaled;
double* en_distance_rescaled_deriv_e;
int32_t uninitialized;
bool provided;
} qmckl_electron_struct;
#+end_src
@ -146,10 +146,10 @@ typedef struct qmckl_electron_struct {
Some values are initialized by default, and are not concerned by
this mechanism.
#+begin_src c :comments org :tangle (eval h_private_func)
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_exit_code qmckl_init_electron(qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c)
qmckl_exit_code qmckl_init_electron(qmckl_context context) {
@ -169,8 +169,8 @@ qmckl_exit_code qmckl_init_electron(qmckl_context context) {
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src c :comments org :tangle (eval h_func)
bool qmckl_electron_provided (const qmckl_context context);
#+end_src
@ -204,7 +204,7 @@ if ( (ctx->electron.uninitialized & mask) != 0) {
return NULL;
}
#+end_src
*** Number of electrons
#+begin_src c :comments org :tangle (eval h_func) :exports none
@ -395,7 +395,7 @@ qmckl_get_electron_rescale_factor_en (const qmckl_context context, double* const
*** Electron coordinates
Returns the current electron coordinates. The pointer is assumed
to point on a memory block of size ~size_max~ \ge ~3 * elec_num * walk_num~.
to point on a memory block of size ~size_max~ \ge ~3 * elec_num * walk_num~.
The order of the indices is:
| | Normal | Transposed |
@ -412,6 +412,10 @@ qmckl_get_electron_coord (const qmckl_context context,
const int64_t size_max);
#+end_src
As the ~coord_new~ attribute is a pointer equal to ~points~,
returning the current electron coordinates is equivalent to
returning the current points.
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_electron_coord (const qmckl_context context,
@ -419,10 +423,6 @@ qmckl_get_electron_coord (const qmckl_context context,
double* const coord,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_INVALID_CONTEXT;
}
if (transp != 'N' && transp != 'T') {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
@ -437,46 +437,32 @@ qmckl_get_electron_coord (const qmckl_context context,
"coord is a null pointer");
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
if ( !(ctx->electron.provided) ) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_get_electron_coord",
"electron data is not provided");
}
int64_t elec_num = ctx->electron.num;
int64_t walk_num = ctx->electron.walk_num;
if (size_max < elec_num*walk_num*3) {
if (size_max <= 0) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_4,
"qmckl_get_electron_coord",
"Array too small");
"size_max should be > 0");
}
assert (ctx->electron.coord_new != NULL);
double* ptr1 = ctx->electron.coord_new;
double* ptr2 = coord;
if (transp == 'N') {
for (int64_t i=0 ; 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 {
memcpy(ptr2, ptr1, 3*elec_num*walk_num*sizeof(double));
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_INVALID_CONTEXT;
}
return QMCKL_SUCCESS;
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
if (!ctx->electron.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_get_electron_coord",
NULL);
}
assert (ctx->point.coord.data == ctx->electron.coord_new.data);
assert (ctx->point.coord.size[0] == ctx->electron.coord_new.size[0]);
assert (ctx->point.coord.size[1] == ctx->electron.coord_new.size[1]);
return qmckl_get_point(context, transp, coord, size_max);
}
#+end_src
@ -512,20 +498,23 @@ ctx->electron.uninitialized &= ~mask;
ctx->electron.provided = (ctx->electron.uninitialized == 0);
if (ctx->electron.provided) {
if (ctx->electron.coord_new != NULL) {
qmckl_free(context, ctx->electron.coord_new);
ctx->electron.coord_new = NULL;
if (ctx->electron.coord_new.data != NULL) {
const qmckl_exit_code rc = qmckl_matrix_free(context, ctx->electron.coord_new);
assert (rc == QMCKL_SUCCESS);
}
if (ctx->electron.coord_old != NULL) {
qmckl_free(context, ctx->electron.coord_old);
ctx->electron.coord_old = NULL;
if (ctx->electron.coord_old.data != NULL) {
const qmckl_exit_code rc = qmckl_matrix_free(context, ctx->electron.coord_old);
assert (rc == QMCKL_SUCCESS);
}
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->electron.num * ctx->electron.walk_num * 3 * sizeof(double);
const int64_t walk_num = ctx->electron.walk_num;
const int64_t elec_num = ctx->electron.num;
double* coord_new = (double*) qmckl_malloc(context, mem_info);
if (coord_new == NULL) {
const int64_t point_num = walk_num * elec_num;
qmckl_matrix coord_new = qmckl_matrix_alloc(context, point_num, 3);
if (coord_new.data == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_set_electron_num",
@ -533,8 +522,9 @@ if (ctx->electron.provided) {
}
ctx->electron.coord_new = coord_new;
double* coord_old = (double*) qmckl_malloc(context, mem_info);
if (coord_old == NULL) {
qmckl_matrix coord_old = qmckl_matrix_alloc(context, point_num, 3);
if (coord_old.data == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_set_electron_num",
@ -542,6 +532,8 @@ if (ctx->electron.provided) {
}
ctx->electron.coord_old = coord_old;
ctx->point.num = point_num;
ctx->point.coord = coord_new;
}
return QMCKL_SUCCESS;
@ -648,9 +640,9 @@ interface
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(in) , value :: alpha
integer (c_int64_t) , intent(in) , value :: beta
integer (c_int64_t) , intent(in) , value :: beta
end function
end interface
@ -660,7 +652,7 @@ interface
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(in) , value :: walk_num
end function
end interface
@ -672,6 +664,9 @@ end interface
overwritten. This can be done only when the data relative to
electrons have been set.
~size_max~ should be equal to ~elec_num * walk_num * 3~, to be symmetric
with ~qmckl_get_electron_coord~.
Important: changing the electron coordinates increments the date
in the context.
@ -699,10 +694,7 @@ qmckl_set_electron_coord(qmckl_context context,
"coord is a null pointer");
}
int64_t elec_num;
qmckl_exit_code rc;
rc = qmckl_get_electron_num(context, &elec_num);
if (rc != QMCKL_SUCCESS) return rc;
const int64_t elec_num = ctx->electron.num;
if (elec_num == 0L) {
return qmckl_failwith( context,
@ -711,9 +703,9 @@ qmckl_set_electron_coord(qmckl_context context,
"elec_num is not set");
}
int64_t walk_num;
rc = qmckl_get_electron_walk_num(context, &walk_num);
if (rc != QMCKL_SUCCESS) return rc;
const int64_t walk_num = ctx->electron.walk_num;
const int64_t point_num = ctx->point.num ;
if (walk_num == 0L) {
return qmckl_failwith( context,
@ -722,40 +714,18 @@ qmckl_set_electron_coord(qmckl_context context,
"walk_num is not set");
}
if (size_max < walk_num*3*elec_num) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_4,
"qmckl_set_electron_coord",
"Array too small");
}
/* If num and walk_num are set, the arrays should be allocated */
assert (ctx->electron.coord_old != NULL);
assert (ctx->electron.coord_new != NULL);
/* Increment the date of the context */
ctx->date += 1UL;
assert(point_num == elec_num * walk_num);
/* Swap pointers */
double * swap;
swap = ctx->electron.coord_old;
ctx->electron.coord_old = ctx->electron.coord_new;
ctx->electron.coord_new = swap;
ctx->point.coord = ctx->electron.coord_old;
ctx->electron.coord_old = ctx->electron.coord_new ;
double* ptr1 = ctx->electron.coord_new;
if (transp == 'N') {
qmckl_exit_code rc;
rc = qmckl_set_point(context, transp, coord, size_max/3);
assert (rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i<walk_num*elec_num ; ++i) {
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 {
ctx->electron.coord_new = ctx->point.coord ;
memcpy(ptr1, coord, 3*elec_num*walk_num*sizeof(double));
}
ctx->electron.coord_new_date = ctx->date;
return QMCKL_SUCCESS;
@ -770,8 +740,8 @@ interface
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
character , intent(in) , value :: transp
integer (c_int64_t) , intent(in) , value :: context
character , intent(in) , value :: transp
double precision , intent(in) :: coord(*)
integer (c_int64_t) , intent(in) , value :: size_max
end function
@ -880,6 +850,7 @@ double elec_coord2[walk_num*3*elec_num];
rc = qmckl_get_electron_coord (context, 'N', elec_coord2, walk_num*3*elec_num);
assert(rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i<3*elec_num*walk_num ; ++i) {
printf("%f %f\n", elec_coord[i], elec_coord2[i]);
assert( elec_coord[i] == elec_coord2[i] );
}
@ -894,7 +865,7 @@ for (int64_t i=0 ; i<3*elec_num*walk_num ; ++i) {
the dependencies are more recent than the date of the data to
compute. If it is the case, then the data is recomputed and the
current date is stored.
** Electron-electron distances
*** Get
@ -981,7 +952,7 @@ qmckl_exit_code qmckl_provide_ee_distance(qmckl_context context)
qmckl_compute_ee_distance(context,
ctx->electron.num,
ctx->electron.walk_num,
ctx->electron.coord_new,
ctx->electron.coord_new.data,
ctx->electron.ee_distance);
if (rc != QMCKL_SUCCESS) {
return rc;
@ -1088,7 +1059,7 @@ qmckl_exit_code qmckl_compute_ee_distance (
#+end_src
*** Test
#+begin_src python :results output :exports none
import numpy as np
@ -1223,7 +1194,7 @@ qmckl_exit_code qmckl_provide_ee_distance_rescaled(qmckl_context context)
ctx->electron.num,
ctx->electron.rescale_factor_kappa_en,
ctx->electron.walk_num,
ctx->electron.coord_new,
ctx->electron.coord_new.data,
ctx->electron.ee_distance_rescaled);
if (rc != QMCKL_SUCCESS) {
return rc;
@ -1389,12 +1360,12 @@ assert(fabs(ee_distance_rescaled[elec_num*elec_num+1]-0.9985724058042633) < 1.e-
#+end_src
** Electron-electron rescaled distance gradients and laplacian with respect to electron coords
** Electron-electron rescaled distance gradients and laplacian with respect to electron coords
The rescaled distances which is given as $R = (1 - \exp{-\kappa r})/\kappa$
The rescaled distances which is given as $R = (1 - \exp{-\kappa r})/\kappa$
needs to be perturbed with respect to the electorn coordinates.
This data is stored in the ~ee_distance_rescaled_deriv_e~ tensor. The
The first three elements of this three index tensor ~[4][num][num]~ gives the
This data is stored in the ~ee_distance_rescaled_deriv_e~ tensor. The
The first three elements of this three index tensor ~[4][num][num]~ gives the
derivatives in the x, y, and z directions $dx, dy, dz$ and the last index
gives the Laplacian $\partial x^2 + \partial y^2 + \partial z^2$.
@ -1427,7 +1398,7 @@ qmckl_exit_code qmckl_get_electron_ee_distance_rescaled_deriv_e(qmckl_context co
#+end_src
*** Provide :noexport:
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_provide_ee_distance_rescaled_deriv_e(qmckl_context context);
#+end_src
@ -1469,7 +1440,7 @@ qmckl_exit_code qmckl_provide_ee_distance_rescaled_deriv_e(qmckl_context context
ctx->electron.num,
ctx->electron.rescale_factor_kappa_en,
ctx->electron.walk_num,
ctx->electron.coord_new,
ctx->electron.coord_new.data,
ctx->electron.ee_distance_rescaled_deriv_e);
if (rc != QMCKL_SUCCESS) {
return rc;
@ -1617,7 +1588,7 @@ rc = qmckl_get_electron_ee_distance_rescaled_deriv_e(context, ee_distance_rescal
#+end_src
** Electron-electron potential
~ee_pot~ calculates the ~ee~ potential energy.
\[
@ -1654,7 +1625,7 @@ qmckl_exit_code qmckl_get_electron_ee_potential(qmckl_context context, double* c
return QMCKL_SUCCESS;
}
#+end_src
*** Provide :noexport:
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
@ -1783,7 +1754,7 @@ end function qmckl_compute_ee_potential_f
const int64_t elec_num,
const int64_t walk_num,
const double* ee_distance,
double* const ee_pot );
double* const ee_pot );
#+end_src
#+CALL: generate_c_interface(table=qmckl_ee_potential_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
@ -1910,7 +1881,7 @@ qmckl_exit_code qmckl_provide_en_distance(qmckl_context context)
ctx->electron.num,
ctx->nucleus.num,
ctx->electron.walk_num,
ctx->electron.coord_new,
ctx->electron.coord_new.data,
ctx->nucleus.coord.data,
ctx->electron.en_distance);
if (rc != QMCKL_SUCCESS) {
@ -2101,7 +2072,7 @@ assert(fabs(en_distance[1][0][1] - 3.1804527583077356) < 1.e-12);
** Electron-nucleus rescaled distances
~en_distance_rescaled~ stores the matrix of the rescaled distances between
~en_distance_rescaled~ stores the matrix of the rescaled distances between
electrons and nuclei.
\[
@ -2187,7 +2158,7 @@ qmckl_exit_code qmckl_provide_en_distance_rescaled(qmckl_context context)
ctx->nucleus.num,
ctx->electron.rescale_factor_kappa_en,
ctx->electron.walk_num,
ctx->electron.coord_new,
ctx->electron.coord_new.data,
ctx->nucleus.coord.data,
ctx->electron.en_distance_rescaled);
if (rc != QMCKL_SUCCESS) {
@ -2374,7 +2345,6 @@ assert (rc == QMCKL_SUCCESS);
assert(fabs(en_distance_rescaled[0][0][0] - 0.9994721712909764) < 1.e-12);
// (1,2,1)
printf("%f\n%f\n", en_distance_rescaled[0][1][0] , 0.9998448354439821);
assert(fabs(en_distance_rescaled[0][1][0] - 0.9998448354439821) < 1.e-12);
// (2,1,1)
@ -2393,10 +2363,10 @@ assert(fabs(en_distance_rescaled[1][0][1] - 0.9584331688679852) < 1.e-12);
** Electron-nucleus rescaled distance gradients and laplacian with respect to electron coords
The rescaled distances which is given as $R = (1 - \exp{-\kappa r})/\kappa$
The rescaled distances which is given as $R = (1 - \exp{-\kappa r})/\kappa$
needs to be perturbed with respect to the nuclear coordinates.
This data is stored in the ~en_distance_rescaled_deriv_e~ tensor. The
The first three elements of this three index tensor ~[4][nucl_num][elec_num]~ gives the
This data is stored in the ~en_distance_rescaled_deriv_e~ tensor. The
The first three elements of this three index tensor ~[4][nucl_num][elec_num]~ gives the
derivatives in the x, y, and z directions $dx, dy, dz$ and the last index
gives the Laplacian $\partial x^2 + \partial y^2 + \partial z^2$.
@ -2476,7 +2446,7 @@ qmckl_exit_code qmckl_provide_en_distance_rescaled_deriv_e(qmckl_context context
ctx->nucleus.num,
ctx->electron.rescale_factor_kappa_en,
ctx->electron.walk_num,
ctx->electron.coord_new,
ctx->electron.coord_new.data,
ctx->nucleus.coord.data,
ctx->electron.en_distance_rescaled_deriv_e);
if (rc != QMCKL_SUCCESS) {
@ -2608,7 +2578,7 @@ qmckl_exit_code qmckl_compute_en_distance_rescaled_deriv_e (
#+end_src
*** Test
#+begin_src python :results output :exports none
import numpy as np
@ -2671,7 +2641,7 @@ assert (rc == QMCKL_SUCCESS);
where \(\mathcal{V}_{en}\) is the ~en~ potential, \[r_{iA}\] the ~en~
distance and \[Z_A\] is the nuclear charge.
*** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
@ -2699,7 +2669,7 @@ qmckl_exit_code qmckl_get_electron_en_potential(qmckl_context context, double* c
return QMCKL_SUCCESS;
}
#+end_src
*** Provide :noexport:
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
@ -2837,7 +2807,7 @@ end function qmckl_compute_en_potential_f
const int64_t walk_num,
const double* charge,
const double* en_distance,
double* const en_pot );
double* const en_pot );
#+end_src
#+CALL: generate_c_interface(table=qmckl_en_potential_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
@ -2873,7 +2843,7 @@ double en_pot[walk_num];
rc = qmckl_get_electron_en_potential(context, &(en_pot[0]));
assert (rc == QMCKL_SUCCESS);
#+end_src
** Generate initial coordinates
*** Compute :noexport:
@ -2887,11 +2857,11 @@ subroutine draw_init_points
integer :: iwalk
logical, allocatable :: do_elec(:)
integer :: acc_num
real, allocatable :: xmin(:,:)
integer :: i, j, k, l, kk
real :: norm
allocate (do_elec(elec_num), xmin(3,elec_num))
xmin = -huge(1.)
@ -2905,7 +2875,7 @@ subroutine draw_init_points
call rinfo( irp_here, 'Norm : ', norm )
call rinfo( irp_here, 'mo_scale: ' , mo_scale )
mo_coef_transp = mo_coef_transp/norm
double precision :: qmc_ranf
real :: mo_max
do i=1,elec_alpha_num
@ -2922,7 +2892,7 @@ subroutine draw_init_points
xmin(2,i) = nucl_coord(l,2)
xmin(3,i) = nucl_coord(l,3)
enddo
call iinfo(irp_here, 'Det num = ', det_num )
do k=1,elec_beta_num
i = k+elec_alpha_num
@ -2939,10 +2909,10 @@ subroutine draw_init_points
xmin(2,i) = nucl_coord(l,2)
xmin(3,i) = nucl_coord(l,3)
enddo
call rinfo( irp_here, 'time step =', time_step )
do iwalk=1,walk_num
print *, 'Generating initial positions for walker', iwalk
print *, 'Generating initial positions for walker', iwalk
acc_num = 0
do_elec = .True.
integer :: iter
@ -2964,7 +2934,7 @@ subroutine draw_init_points
TOUCH elec_coord
re_compute = minval(nucl_elec_dist(1:nucl_num,1:elec_num))
enddo
do i=1,elec_alpha_num
if (do_elec(i)) then
if ( mo_value_transp(i,i)**2 >= qmc_ranf()) then
@ -2973,7 +2943,7 @@ subroutine draw_init_points
endif
endif
enddo
do i=1,elec_beta_num
if (do_elec(i+elec_alpha_num)) then
if ( mo_value_transp(i,i+elec_alpha_num)**2 >= qmc_ranf()) then
@ -2982,9 +2952,9 @@ subroutine draw_init_points
endif
endif
enddo
enddo
do l=1,3
do i=1,elec_num+1
elec_coord_full(i,l,iwalk) = elec_coord(i,l)
@ -2997,7 +2967,7 @@ subroutine draw_init_points
endif
SOFT_TOUCH elec_coord elec_coord_full
deallocate (do_elec, xmin)
end
# end_src

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

94
org/qmckl_point.org
View File

@ -77,10 +77,11 @@ int main() {
The following data stored in the context:
| Variable | Type | Description |
|----------+----------------+------------------------|
| ~num~ | ~int64_t~ | Total number of points |
| ~coord~ | ~qmckl_matrix~ | ~num~ \times 3 matrix3 |
| Variable | Type | Description |
|----------+----------------+-------------------------------------------|
| ~num~ | ~int64_t~ | Total number of points |
| ~date~ | ~uint64_t~ | Last modification date of the coordinates |
| ~coord~ | ~qmckl_matrix~ | ~num~ \times 3 matrix3 |
We consider that the matrix is stored 'transposed' and 'normal'
corresponds to the 3 \times ~num~ matrix.
@ -90,15 +91,16 @@ int main() {
#+begin_src c :comments org :tangle (eval h_private_type)
typedef struct qmckl_point_struct {
int64_t num;
uint64_t date;
qmckl_matrix coord;
} qmckl_point_struct;
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func)
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_exit_code qmckl_init_point(qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c)
qmckl_exit_code qmckl_init_point(qmckl_context context) {
@ -109,16 +111,7 @@ qmckl_exit_code qmckl_init_point(qmckl_context context) {
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = sizeof(qmckl_point_struct);
ctx->point = (qmckl_point_struct*) qmckl_malloc(context, mem_info);
if (ctx->point == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_init_point",
NULL);
}
memset(ctx->point, 0, sizeof(qmckl_point_struct));
memset(&(ctx->point), 0, sizeof(qmckl_point_struct));
return QMCKL_SUCCESS;
}
@ -139,7 +132,7 @@ qmckl_exit_code qmckl_get_point_num (const qmckl_context context, int64_t* const
#+end_src
Returns the number of points stored in the context.
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_point_num (const qmckl_context context, int64_t* const num) {
@ -157,10 +150,9 @@ qmckl_get_point_num (const qmckl_context context, int64_t* const num) {
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
assert (ctx->point != NULL);
assert (ctx->point->num > (int64_t) 0);
,*num = ctx->point->num;
assert (ctx->point.num >= (int64_t) 0);
,*num = ctx->point.num;
return QMCKL_SUCCESS;
}
#+end_src
@ -172,14 +164,14 @@ interface
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(out) :: num
integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(out) :: num
end function
end interface
#+end_src
*** Point coordinates
#+begin_src c :comments org :tangle (eval h_func) :exports none
qmckl_exit_code qmckl_get_point(const qmckl_context context,
const char transp,
@ -212,12 +204,11 @@ qmckl_get_point(const qmckl_context context,
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
assert (ctx->point != NULL);
int64_t point_num = ctx->point->num;
int64_t point_num = ctx->point.num;
if (point_num == 0) return QMCKL_NOT_PROVIDED;
assert (ctx->point->coord.data != NULL);
assert (ctx->point.coord.data != NULL);
if (size_max < 3*point_num) {
return qmckl_failwith( context,
@ -229,16 +220,16 @@ qmckl_get_point(const qmckl_context context,
qmckl_exit_code rc;
if (transp == 'N') {
qmckl_matrix At = qmckl_matrix_alloc( context, 3, point_num);
rc = qmckl_transpose( context, ctx->point->coord, At );
rc = qmckl_transpose( context, ctx->point.coord, At );
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_double_of_matrix( context, At, coord, size_max);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_matrix_free(context, At);
} else {
rc = qmckl_double_of_matrix( context, ctx->point->coord, coord, size_max);
rc = qmckl_double_of_matrix( context, ctx->point.coord, coord, size_max);
}
if (rc != QMCKL_SUCCESS) return rc;
return rc;
}
@ -259,23 +250,23 @@ interface
end interface
#+end_src
** Initialization functions
When the data is set in the context, if the arrays are large
enough, we overwrite the data contained in them.
To set the data relative to the points in the context, one of the
following functions need to be called.
following functions need to be called.
#+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code qmckl_set_point (qmckl_context context,
const char transp,
const double* coord,
const int64_t size_max);
const int64_t num);
#+end_src
Copy a sequence of (x,y,z) into the context.
Copy a sequence of ~num~ points $(x,y,z)$ into the context.
#+begin_src c :comments org :tangle (eval c) :noweb yes
qmckl_exit_code
@ -284,7 +275,7 @@ qmckl_set_point (qmckl_context context,
const double* coord,
const int64_t num)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
@ -295,28 +286,27 @@ qmckl_set_point (qmckl_context context,
"qmckl_set_point",
"transp should be 'N' or 'T'");
}
if (coord == NULL) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_set_point",
"coord is a NULL pointer");
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
assert (ctx->point != NULL);
qmckl_exit_code rc;
if (ctx->point->num < num) {
if (ctx->point.num < num) {
if (ctx->point->coord.data != NULL) {
rc = qmckl_matrix_free(context, ctx->point->coord);
return rc;
if (ctx->point.coord.data != NULL) {
rc = qmckl_matrix_free(context, ctx->point.coord);
assert (rc == QMCKL_SUCCESS);
}
ctx->point->coord = qmckl_matrix_alloc(context, num, 3);
if (ctx->point->coord.data == NULL) {
ctx->point.coord = qmckl_matrix_alloc(context, num, 3);
if (ctx->point.coord.data == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_set_point",
@ -325,18 +315,22 @@ qmckl_set_point (qmckl_context context,
};
ctx->point->num = num;
ctx->point.num = num;
if (transp == 'T') {
memcpy(ctx->point->coord.data, coord, 3*num*sizeof(double));
memcpy(ctx->point.coord.data, coord, 3*num*sizeof(double));
} else {
for (int64_t i=0 ; i<num ; ++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, 2) = coord[3*i+2];
qmckl_mat(ctx->point.coord, i, 0) = coord[3*i ];
qmckl_mat(ctx->point.coord, i, 1) = coord[3*i+1];
qmckl_mat(ctx->point.coord, i, 2) = coord[3*i+2];
}
}
/* Increment the date of the context */
ctx->date += 1UL;
ctx->point.date = ctx->date;
return QMCKL_SUCCESS;
}