TREX/qmckl
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-06-30 00:44:52 +02:00
Code Issues Projects Releases Wiki Activity

Merge branch 'TREX-CoE:master' into master

Browse Source
This commit is contained in:
vijay 2022-01-26 18:34:48 +01:00 committed by GitHub
commit 402fe7c9c8
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
12 changed files with 1337 additions and 1178 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

480
org/qmckl_blas.org
View File

@ -38,8 +38,9 @@
#include "qmckl.h"
#include "qmckl_context_private_type.h"
#include "qmckl_memory_private_type.h"
#include "qmckl_memory_private_func.h"
#include "qmckl_blas_private_type.h"
#include "qmckl_memory_private_func.h"
#include "qmckl_blas_private_func.h"
#+end_src
@ -379,23 +380,18 @@ qmckl_tensor_of_vector(const qmckl_vector vector,
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_vector
qmckl_vector_of_matrix(const qmckl_matrix matrix,
const int64_t size);
qmckl_vector_of_matrix(const qmckl_matrix matrix);
#+end_src
Reshapes a matrix into a vector.
#+begin_src c :comments org :tangle (eval c)
qmckl_vector
qmckl_vector_of_matrix(const qmckl_matrix matrix,
const int64_t size)
qmckl_vector_of_matrix(const qmckl_matrix matrix)
{
/* Always true */
assert (matrix.size[0] * matrix.size[1] == size);
qmckl_vector result;
result.size = size;
result.size = matrix.size[0] * matrix.size[1];
result.data = matrix.data;
return result;
@ -438,27 +434,23 @@ qmckl_tensor_of_matrix(const qmckl_matrix matrix,
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_vector
qmckl_vector_of_tensor(const qmckl_tensor tensor,
const int64_t size);
qmckl_vector_of_tensor(const qmckl_tensor tensor);
#+end_src
Reshapes a tensor into a vector.
#+begin_src c :comments org :tangle (eval c)
qmckl_vector
qmckl_vector_of_tensor(const qmckl_tensor tensor,
const int64_t size)
qmckl_vector_of_tensor(const qmckl_tensor tensor)
{
/* Always true */
int64_t prod_size = (int64_t) 1;
for (int64_t i=0 ; i<tensor.order ; i++) {
int64_t prod_size = (int64_t) tensor.size[0];
for (int64_t i=1 ; i<tensor.order ; i++) {
prod_size *= tensor.size[i];
}
assert (prod_size == size);
qmckl_vector result;
result.size = size;
result.size = prod_size;
result.data = tensor.data;
return result;
@ -510,6 +502,186 @@ qmckl_matrix_of_tensor(const qmckl_tensor tensor,
#define qmckl_ten5(t, i, j, k, l, m) t.data[(i) + m.size[0]*((j) + size[1]*((k) + size[2]*((l) + size[3]*(m))))]
#+end_src
** Copy to/from to ~double*~
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_exit_code
qmckl_double_of_vector(const qmckl_context context,
const qmckl_vector vector,
double* const target,
const int64_t size_max);
#+end_src
Converts a vector to a ~double*~.
#+begin_src c :comments org :tangle (eval c)
qmckl_exit_code
qmckl_double_of_vector(const qmckl_context context,
const qmckl_vector vector,
double* const target,
const int64_t size_max)
{
/* Always true by construction */
assert (qmckl_context_check(context) != QMCKL_NULL_CONTEXT);
assert (vector.size > (int64_t) 0);
assert (target != NULL);
assert (size_max > (int64_t) 0);
assert (size_max >= vector.size);
for (int64_t i=0 ; i<vector.size ; ++i) {
target[i] = vector.data[i];
}
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_exit_code
qmckl_double_of_matrix(const qmckl_context context,
const qmckl_matrix matrix,
double* const target,
const int64_t size_max);
#+end_src
Converts a matrix to a ~double*~.
#+begin_src c :comments org :tangle (eval c)
qmckl_exit_code
qmckl_double_of_matrix(const qmckl_context context,
const qmckl_matrix matrix,
double* const target,
const int64_t size_max)
{
qmckl_vector vector = qmckl_vector_of_matrix(matrix);
return qmckl_double_of_vector(context, vector, target, size_max);
}
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_exit_code
qmckl_double_of_tensor(const qmckl_context context,
const qmckl_tensor tensor,
double* const target,
const int64_t size_max);
#+end_src
Converts a tensor to a ~double*~.
#+begin_src c :comments org :tangle (eval c)
qmckl_exit_code
qmckl_double_of_tensor(const qmckl_context context,
const qmckl_tensor tensor,
double* const target,
const int64_t size_max)
{
qmckl_vector vector = qmckl_vector_of_tensor(tensor);
return qmckl_double_of_vector(context, vector, target, size_max);
}
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_exit_code
qmckl_vector_of_double(const qmckl_context context,
const double* target,
const int64_t size_max,
qmckl_vector* vector);
#+end_src
Converts a ~double*~ to a vector.
#+begin_src c :comments org :tangle (eval c)
qmckl_exit_code
qmckl_vector_of_double(const qmckl_context context,
const double* target,
const int64_t size_max,
qmckl_vector* vector_out)
{
qmckl_vector vector = *vector_out;
/* Always true by construction */
assert (qmckl_context_check(context) != QMCKL_NULL_CONTEXT);
if (vector.size == 0) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_4,
"qmckl_double_of_vector",
"Vector not allocated");
}
if (vector.size != size_max) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_4,
"qmckl_double_of_vector",
"Wrong vector size");
}
for (int64_t i=0 ; i<vector.size ; ++i) {
vector.data[i] = target[i];
}
*vector_out = vector;
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_exit_code
qmckl_matrix_of_double(const qmckl_context context,
const double* target,
const int64_t size_max,
qmckl_matrix* matrix);
#+end_src
Converts a matrix to a ~double*~.
#+begin_src c :comments org :tangle (eval c)
qmckl_exit_code
qmckl_matrix_of_double(const qmckl_context context,
const double* target,
const int64_t size_max,
qmckl_matrix* matrix)
{
qmckl_vector vector = qmckl_vector_of_matrix(*matrix);
qmckl_exit_code rc =
qmckl_vector_of_double(context, target, size_max, &vector);
,*matrix = qmckl_matrix_of_vector(vector, matrix->size[0], matrix->size[1]);
return rc;
}
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_exit_code
qmckl_tensor_of_double(const qmckl_context context,
const double* target,
const int64_t size_max,
qmckl_tensor* tensor);
#+end_src
Converts a matrix to a ~double*~.
#+begin_src c :comments org :tangle (eval c)
qmckl_exit_code
qmckl_tensor_of_double(const qmckl_context context,
const double* target,
const int64_t size_max,
qmckl_tensor* tensor)
{
qmckl_vector vector = qmckl_vector_of_tensor(*tensor);
qmckl_exit_code rc =
qmckl_vector_of_double(context, target, size_max, &vector);
*tensor = qmckl_tensor_of_vector(vector, tensor->order, tensor->size);
return rc;
}
#+end_src
** Tests
#+begin_src c :comments link :tangle (eval c_test)
@ -534,7 +706,7 @@ qmckl_matrix_of_tensor(const qmckl_tensor tensor,
for (int64_t i=0 ; i<m ; ++i)
assert ( qmckl_mat(mat, i, j) == qmckl_vec(vec, i+j*m)) ;
qmckl_vector vec2 = qmckl_vector_of_matrix(mat, p);
qmckl_vector vec2 = qmckl_vector_of_matrix(mat);
assert (vec2.size == p);
assert (vec2.data == vec.data);
for (int64_t i=0 ; i<p ; ++i)
@ -565,14 +737,14 @@ qmckl_matrix_of_tensor(const qmckl_tensor tensor,
| ~m~ | ~int64_t~ | in | Number of rows of the input matrix |
| ~n~ | ~int64_t~ | in | Number of columns of the input matrix |
| ~k~ | ~int64_t~ | in | Number of columns of the input matrix |
| ~alpha~ | ~double~ | in | Number of columns of the input matrix |
| ~alpha~ | ~double~ | in | \alpha |
| ~A~ | ~double[][lda]~ | in | Array containing the matrix $A$ |
| ~lda~ | ~int64_t~ | in | Leading dimension of array ~A~ |
| ~B~ | ~double[][ldb]~ | in | Array containing the matrix $B$ |
| ~ldb~ | ~int64_t~ | in | Leading dimension of array ~B~ |
| ~beta~ | ~double~ | in | Array containing the matrix $B$ |
| ~C~ | ~double[][ldc]~ | out | Array containing the matrix $B$ |
| ~ldc~ | ~int64_t~ | in | Leading dimension of array ~B~ |
| ~beta~ | ~double~ | in | \beta |
| ~C~ | ~double[][ldc]~ | out | Array containing the matrix $C$ |
| ~ldc~ | ~int64_t~ | in | Leading dimension of array ~C~ |
Requirements:
@ -808,6 +980,184 @@ qmckl_exit_code test_qmckl_dgemm(qmckl_context context);
assert(QMCKL_SUCCESS == test_qmckl_dgemm(context));
#+end_src
** ~qmckl_matmul~
Matrix multiplication:
\[
C_{ij} = \beta C_{ij} + \alpha \sum_{k} A_{ik} \cdot B_{kj}
\]
# TODO: Add description about the external library dependence.
#+NAME: qmckl_matmul_args
| Variable | Type | In/Out | Description |
|-----------+-----------------+--------+-------------------|
| ~context~ | ~qmckl_context~ | in | Global state |
| ~TransA~ | ~char~ | in | 'T' is transposed |
| ~TransB~ | ~char~ | in | 'T' is transposed |
| ~alpha~ | ~double~ | in | \alpha |
| ~A~ | ~qmckl_matrix~ | in | Matrix $A$ |
| ~B~ | ~qmckl_matrix~ | in | Matrix $B$ |
| ~beta~ | ~double~ | in | \beta |
| ~C~ | ~qmckl_matrix~ | out | Matrix $C$ |
#+CALL: generate_c_header(table=qmckl_matmul_args,rettyp="qmckl_exit_code",fname="qmckl_matmul")
#+RESULTS:
#+begin_src c :tangle (eval h_private_func) :comments org
qmckl_exit_code
qmckl_matmul (const qmckl_context context,
const char TransA,
const char TransB,
const double alpha,
const qmckl_matrix A,
const qmckl_matrix B,
const double beta,
qmckl_matrix* const C );
#+end_src
#+begin_src c :tangle (eval c) :comments org
qmckl_exit_code
qmckl_matmul (const qmckl_context context,
const char TransA,
const char TransB,
const double alpha,
const qmckl_matrix A,
const qmckl_matrix B,
const double beta,
qmckl_matrix* const C )
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_INVALID_CONTEXT;
}
qmckl_exit_code rc = QMCKL_SUCCESS;
if (TransA != 'N' && TransA != 'T') {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
"qmckl_matmul",
"TransA should be 'N' or 'T'");
}
if (TransB != 'N' && TransB != 'T') {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_matmul",
"TransB should be 'N' or 'T'");
}
if (A.size[0] < 1) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_5,
"qmckl_matmul",
"Invalid size for A");
}
if (B.size[0] < 1) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_6,
"qmckl_matmul",
"Invalid size for B");
}
if (C == NULL) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_8,
"qmckl_matmul",
"Null pointer");
}
int t = 0;
if (TransA == 'T') t +=1;
if (TransB == 'T') t +=2;
/*
| t | TransA | TransB |
+---+--------+--------+
| 0 | N | N |
| 1 | T | N |
| 2 | N | T |
| 3 | T | T |
,*/
switch (t) {
case 0:
if (A.size[1] != B.size[0]) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
"qmckl_matmul",
"A and B have incompatible dimensions");
}
C->size[0] = A.size[0];
C->size[1] = B.size[1];
rc = qmckl_dgemm (context, 'N', 'N',
C->size[0], C->size[1], A.size[1],
alpha,
A.data, A.size[0],
B.data, B.size[0],
beta,
C->data, C->size[0]);
break;
case 1:
if (A.size[0] != B.size[0]) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
"qmckl_matmul",
"A and B have incompatible dimensions");
}
C->size[0] = A.size[1];
C->size[1] = B.size[1];
rc = qmckl_dgemm (context, 'T', 'N',
C->size[0], C->size[1], A.size[0],
alpha,
A.data, A.size[0],
B.data, B.size[0],
beta,
C->data, C->size[0]);
break;
case 2:
if (A.size[1] != B.size[1]) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
"qmckl_matmul",
"A and B have incompatible dimensions");
}
C->size[0] = A.size[0];
C->size[1] = B.size[0];
rc = qmckl_dgemm (context, 'N', 'T',
C->size[0], C->size[1], A.size[1],
alpha,
A.data, A.size[0],
B.data, B.size[0],
beta,
C->data, C->size[0]);
break;
case 3:
if (A.size[0] != B.size[1]) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
"qmckl_matmul",
"A and B have incompatible dimensions");
}
C->size[0] = A.size[1];
C->size[1] = B.size[0];
rc = qmckl_dgemm (context, 'T', 'T',
C->size[0], C->size[1], A.size[0],
alpha,
A.data, A.size[0],
B.data, B.size[0],
beta,
C->data, C->size[0]);
break;
}
return rc;
}
#+end_src
*** TODO Test :noexport:
** ~qmckl_adjugate~
Given a matrix $\mathbf{A}$, the adjugate matrix
@ -1693,6 +2043,90 @@ qmckl_exit_code test_qmckl_adjugate(qmckl_context context);
assert(QMCKL_SUCCESS == test_qmckl_adjugate(context));
#+end_src
** ~qmckl_transpose~
Transposes a matrix: $A^\dagger_{ji} = A_{ij}$.
| Variable | Type | In/Out | Description |
|----------+---------------+--------+-------------------|
| context | qmckl_context | in | Global state |
| A | qmckl_matrix | in | Input matrix |
| At | qmckl_matrix | out | Transposed matrix |
#+begin_src c :tangle (eval h_private_func) :comments org
qmckl_exit_code
qmckl_transpose (qmckl_context context,
const qmckl_matrix A,
qmckl_matrix At );
#+end_src
#+begin_src c :tangle (eval c) :comments org
qmckl_exit_code
qmckl_transpose (qmckl_context context,
const qmckl_matrix A,
qmckl_matrix At )
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_INVALID_CONTEXT;
}
if (A.size[0] < 1) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
"qmckl_transpose",
"Invalid size for A");
}
if (At.data == NULL) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_transpose",
"Output matrix not allocated");
}
if (At.size[0] != A.size[1] || At.size[1] != A.size[0]) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_transpose",
"Invalid size for At");
}
for (int64_t j=0 ; j<At.size[1] ; ++j)
for (int64_t i=0 ; i<At.size[0] ; ++i)
qmckl_mat(At, i, j) = qmckl_mat(A, j, i);
return QMCKL_SUCCESS;
}
#+end_src
*** Test
#+begin_src c :comments link :tangle (eval c_test)
{
qmckl_matrix A;
qmckl_matrix At;
A = qmckl_matrix_alloc(context, 2, 3);
At = qmckl_matrix_alloc(context, 3, 2);
for (int j=0 ; j<3 ; ++j)
for (int i=0 ; i<2 ; ++i)
qmckl_mat(A, i, j) = (double) 10*i+j;
qmckl_exit_code rc = qmckl_transpose(context, A, At);
assert(rc == QMCKL_SUCCESS);
assert(A.size[0] == At.size[1]);
assert(A.size[1] == At.size[0]);
for (int j=0 ; j<3 ; ++j)
for (int i=0 ; i<2 ; ++i)
assert (qmckl_mat(A, i, j) == qmckl_mat(At, j, i));
qmckl_matrix_free(context, A);
qmckl_matrix_free(context, At);
}
#+end_src
* End of files :noexport:

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;

4
org/qmckl_determinant.org
View File

@ -1160,10 +1160,10 @@ assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_num (context, nucl_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]));
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), nucl_num*3);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_charge(context, nucl_charge);
rc = qmckl_set_nucleus_charge(context, nucl_charge, nucl_num);
assert(rc == QMCKL_SUCCESS);
assert(qmckl_nucleus_provided(context));

393
org/qmckl_electron.org
View File

@ -26,6 +26,7 @@ up-spin and down-spin electrons, and the electron coordinates.
#+begin_src c :tangle (eval c_test) :noweb yes
#include "qmckl.h"
#include <assert.h>
#include <stdio.h>
#include <math.h>
#ifdef HAVE_CONFIG_H
#include "config.h"
@ -68,22 +69,22 @@ int main() {
The following data stored in the context:
| Variable | Type | Description |
|---------------------------+----------------------------+-------------------------------------------|
| ~uninitialized~ | ~int32_t~ | Keeps bit set for uninitialized data |
| ~num~ | ~int64_t~ | Total number of electrons |
| ~up_num~ | ~int64_t~ | Number of up-spin electrons |
| ~down_num~ | ~int64_t~ | Number of down-spin electrons |
| ~walk_num~ | ~int64_t~ | Number of walkers |
| ~rescale_factor_kappa_ee~ | ~double~ | The distance scaling factor |
| ~rescale_factor_kappa_en~ | ~double~ | The distance scaling factor |
| ~provided~ | ~bool~ | If true, ~electron~ is valid |
| ~coord_new~ | ~double[walk_num][3][num]~ | New set of electron coordinates |
| ~coord_old~ | ~double[walk_num][3][num]~ | Old set of electron coordinates |
| ~coord_new_date~ | ~uint64_t~ | Last modification date of the coordinates |
| Variable | Type | Description |
|---------------------------+----------------+---------------------------------------------------------------|
| ~uninitialized~ | ~int32_t~ | Keeps bit set for uninitialized data |
| ~num~ | ~int64_t~ | Total number of electrons |
| ~up_num~ | ~int64_t~ | Number of up-spin electrons |
| ~down_num~ | ~int64_t~ | Number of down-spin electrons |
| ~walk_num~ | ~int64_t~ | Number of walkers |
| ~rescale_factor_kappa_ee~ | ~double~ | The distance scaling factor |
| ~rescale_factor_kappa_en~ | ~double~ | The distance scaling factor |
| ~provided~ | ~bool~ | If true, ~electron~ is valid |
| ~coord_new~ | ~qmckl_matrix~ | Current set of electron coordinates. Pointer to ~ctx->points~ |
| ~coord_old~ | ~qmckl_matrix~ | Old set of electron coordinates |
| ~coord_new_date~ | ~uint64_t~ | Last modification date of the coordinates |
Computed data:
| Variable | Type | Description |
|-------------------------------------+--------------------------------------+----------------------------------------------------------------------|
| ~ee_distance~ | ~double[walk_num][num][num]~ | Electron-electron distances |
@ -107,33 +108,33 @@ int main() {
#+begin_src c :comments org :tangle (eval h_private_type)
typedef struct qmckl_electron_struct {
int64_t num;
int64_t up_num;
int64_t down_num;
int64_t walk_num;
double rescale_factor_kappa_ee;
double rescale_factor_kappa_en;
int64_t coord_new_date;
int64_t ee_distance_date;
int64_t en_distance_date;
int64_t ee_pot_date;
int64_t en_pot_date;
int64_t ee_distance_rescaled_date;
int64_t ee_distance_rescaled_deriv_e_date;
int64_t en_distance_rescaled_date;
int64_t en_distance_rescaled_deriv_e_date;
double* coord_new;
double* coord_old;
double* ee_distance;
double* en_distance;
double* ee_pot;
double* en_pot;
double* ee_distance_rescaled;
double* ee_distance_rescaled_deriv_e;
double* en_distance_rescaled;
double* en_distance_rescaled_deriv_e;
int32_t uninitialized;
bool provided;
int64_t num;
int64_t up_num;
int64_t down_num;
int64_t walk_num;
double rescale_factor_kappa_ee;
double rescale_factor_kappa_en;
int64_t coord_new_date;
int64_t ee_distance_date;
int64_t en_distance_date;
int64_t ee_pot_date;
int64_t en_pot_date;
int64_t ee_distance_rescaled_date;
int64_t ee_distance_rescaled_deriv_e_date;
int64_t en_distance_rescaled_date;
int64_t en_distance_rescaled_deriv_e_date;
qmckl_matrix coord_new;
qmckl_matrix coord_old;
double* ee_distance;
double* en_distance;
double* ee_pot;
double* en_pot;
double* ee_distance_rescaled;
double* ee_distance_rescaled_deriv_e;
double* en_distance_rescaled;
double* en_distance_rescaled_deriv_e;
int32_t uninitialized;
bool provided;
} qmckl_electron_struct;
#+end_src
@ -145,10 +146,10 @@ typedef struct qmckl_electron_struct {
Some values are initialized by default, and are not concerned by
this mechanism.
#+begin_src c :comments org :tangle (eval h_private_func)
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_exit_code qmckl_init_electron(qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c)
qmckl_exit_code qmckl_init_electron(qmckl_context context) {
@ -168,8 +169,8 @@ qmckl_exit_code qmckl_init_electron(qmckl_context context) {
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src c :comments org :tangle (eval h_func)
bool qmckl_electron_provided (const qmckl_context context);
#+end_src
@ -203,7 +204,7 @@ if ( (ctx->electron.uninitialized & mask) != 0) {
return NULL;
}
#+end_src
*** Number of electrons
#+begin_src c :comments org :tangle (eval h_func) :exports none
@ -394,27 +395,34 @@ qmckl_get_electron_rescale_factor_en (const qmckl_context context, double* const
*** Electron coordinates
Returns the current electron coordinates. The pointer is assumed
to point on a memory block of size ~size_max~ \ge ~3 * elec_num * walk_num~.
to point on a memory block of size ~size_max~ \ge ~3 * elec_num * walk_num~.
The order of the indices is:
| | Normal | Transposed |
|---------+---------------------------+---------------------------|
| C | ~[walk_num][elec_num][3]~ | ~[walk_num][3][elec_num]~ |
| Fortran | ~(3,elec_num,walk_num)~ | ~(elec_num,3,walk_num)~ |
| | Normal | Transposed |
|---------+--------------------------+--------------------------|
| C | ~[walk_num*elec_num][3]~ | ~[3][walk_num*elec_num]~ |
| Fortran | ~(3,walk_num*elec_num)~ | ~(walk_num*elec_num, 3)~ |
#+begin_src c :comments org :tangle (eval h_func) :exports none
qmckl_exit_code qmckl_get_electron_coord (const qmckl_context context, const char transp, double* const coord);
qmckl_exit_code
qmckl_get_electron_coord (const qmckl_context context,
const char transp,
double* const coord,
const int64_t size_max);
#+end_src
As the ~coord_new~ attribute is a pointer equal to ~points~,
returning the current electron coordinates is equivalent to
returning the current points.
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_electron_coord (const qmckl_context context, const char transp, double* const coord) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_INVALID_CONTEXT;
}
qmckl_get_electron_coord (const qmckl_context context,
const char transp,
double* const coord,
const int64_t size_max)
{
if (transp != 'N' && transp != 'T') {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
@ -429,42 +437,32 @@ qmckl_get_electron_coord (const qmckl_context context, const char transp, double
"coord is a null pointer");
}
if (size_max <= 0) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_4,
"qmckl_get_electron_coord",
"size_max should be > 0");
}
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_INVALID_CONTEXT;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
if ( !(ctx->electron.provided) ) {
if (!ctx->electron.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_get_electron_coord",
"electron data is not provided");
NULL);
}
int64_t elec_num = ctx->electron.num;
int64_t walk_num = ctx->electron.walk_num;
assert (ctx->point.coord.data == ctx->electron.coord_new.data);
assert (ctx->point.coord.size[0] == ctx->electron.coord_new.size[0]);
assert (ctx->point.coord.size[1] == ctx->electron.coord_new.size[1]);
assert (ctx->electron.coord_new != NULL);
double* ptr1 = ctx->electron.coord_new;
double* ptr2 = coord;
if (transp == 'N') {
for (int64_t i=0 ; i<walk_num ; ++i) {
qmckl_exit_code rc;
rc = qmckl_transpose(context, elec_num, 3,
ptr1, elec_num, ptr2, 3);
if (rc != QMCKL_SUCCESS) return rc;
ptr1 += elec_num * 3;
ptr2 += elec_num * 3;
}
} else {
memcpy(ptr2, ptr1, 3*elec_num*walk_num*sizeof(double));
}
return QMCKL_SUCCESS;
return qmckl_get_point(context, transp, coord, size_max);
}
#+end_src
@ -500,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",
@ -521,8 +522,9 @@ if (ctx->electron.provided) {
}
ctx->electron.coord_new = coord_new;
double* coord_old = (double*) qmckl_malloc(context, mem_info);
if (coord_old == NULL) {
qmckl_matrix coord_old = qmckl_matrix_alloc(context, point_num, 3);
if (coord_old.data == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_set_electron_num",
@ -530,6 +532,8 @@ if (ctx->electron.provided) {
}
ctx->electron.coord_old = coord_old;
ctx->point.num = point_num;
ctx->point.coord = coord_new;
}
return QMCKL_SUCCESS;
@ -636,18 +640,19 @@ interface
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(in) , value :: alpha
integer (c_int64_t) , intent(in) , value :: beta
integer (c_int64_t) , intent(in) , value :: beta
end function
end interface
interface
integer(c_int32_t) function qmckl_set_electron_walk_num(context, walk_num) bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(in) , value :: walk_num
end function
end interface
@ -659,6 +664,9 @@ end interface
overwritten. This can be done only when the data relative to
electrons have been set.
~size_max~ should be equal to ~elec_num * walk_num * 3~, to be symmetric
with ~qmckl_get_electron_coord~.
Important: changing the electron coordinates increments the date
in the context.
@ -686,10 +694,7 @@ qmckl_set_electron_coord(qmckl_context context,
"coord is a null pointer");
}
int64_t elec_num;
qmckl_exit_code rc;
rc = qmckl_get_electron_num(context, &elec_num);
if (rc != QMCKL_SUCCESS) return rc;
const int64_t elec_num = ctx->electron.num;
if (elec_num == 0L) {
return qmckl_failwith( context,
@ -698,9 +703,9 @@ qmckl_set_electron_coord(qmckl_context context,
"elec_num is not set");
}
int64_t walk_num;
rc = qmckl_get_electron_walk_num(context, &walk_num);
if (rc != QMCKL_SUCCESS) return rc;
const int64_t walk_num = ctx->electron.walk_num;
const int64_t point_num = ctx->point.num ;
if (walk_num == 0L) {
return qmckl_failwith( context,
@ -709,41 +714,18 @@ qmckl_set_electron_coord(qmckl_context context,
"walk_num is not set");
}
if (size_max < walk_num*3*elec_num) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_4,
"qmckl_set_electron_coord",
"destination array is too small");
}
/* If num and walk_num are set, the arrays should be allocated */
assert (ctx->electron.coord_old != NULL);
assert (ctx->electron.coord_new != NULL);
/* Increment the date of the context */
ctx->date += 1UL;
assert(point_num == elec_num * walk_num);
/* Swap pointers */
double * swap;
swap = ctx->electron.coord_old;
ctx->electron.coord_old = ctx->electron.coord_new;
ctx->electron.coord_new = swap;
ctx->point.coord = ctx->electron.coord_old;
ctx->electron.coord_old = ctx->electron.coord_new ;
double* ptr1 = ctx->electron.coord_new;
if (transp == 'N') {
qmckl_exit_code rc;
rc = qmckl_set_point(context, transp, coord, size_max/3);
assert (rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i<walk_num ; ++i) {
rc = qmckl_transpose(context, 3, elec_num,
&(coord[i*3*elec_num]), 3, ptr1, elec_num);
if (rc != QMCKL_SUCCESS) return rc;
ptr1 += elec_num * 3;
}
} 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;
@ -758,8 +740,8 @@ interface
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
character , intent(in) , value :: transp
integer (c_int64_t) , intent(in) , value :: context
character , intent(in) , value :: transp
double precision , intent(in) :: coord(*)
integer (c_int64_t) , intent(in) , value :: size_max
end function
@ -865,9 +847,10 @@ assert(rc == QMCKL_SUCCESS);
double elec_coord2[walk_num*3*elec_num];
rc = qmckl_get_electron_coord (context, 'N', elec_coord2);
rc = qmckl_get_electron_coord (context, 'N', elec_coord2, walk_num*3*elec_num);
assert(rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i<3*elec_num ; ++i) {
for (int64_t i=0 ; i<3*elec_num*walk_num ; ++i) {
printf("%f %f\n", elec_coord[i], elec_coord2[i]);
assert( elec_coord[i] == elec_coord2[i] );
}
@ -882,7 +865,7 @@ for (int64_t i=0 ; i<3*elec_num ; ++i) {
the dependencies are more recent than the date of the data to
compute. If it is the case, then the data is recomputed and the
current date is stored.
** Electron-electron distances
*** Get
@ -969,7 +952,7 @@ qmckl_exit_code qmckl_provide_ee_distance(qmckl_context context)
qmckl_compute_ee_distance(context,
ctx->electron.num,
ctx->electron.walk_num,
ctx->electron.coord_new,
ctx->electron.coord_new.data,
ctx->electron.ee_distance);
if (rc != QMCKL_SUCCESS) {
return rc;
@ -995,7 +978,7 @@ qmckl_exit_code qmckl_provide_ee_distance(qmckl_context context)
| ~context~ | ~qmckl_context~ | in | Global state |
| ~elec_num~ | ~int64_t~ | in | Number of electrons |
| ~walk_num~ | ~int64_t~ | in | Number of walkers |
| ~coord~ | ~double[walk_num][3][elec_num]~ | in | Electron coordinates |
| ~coord~ | ~double[3][walk_num][elec_num]~ | in | Electron coordinates |
| ~ee_distance~ | ~double[walk_num][elec_num][elec_num]~ | out | Electron-electron distances |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
@ -1006,10 +989,11 @@ integer function qmckl_compute_ee_distance_f(context, elec_num, walk_num, coord,
integer(qmckl_context), intent(in) :: context
integer*8 , intent(in) :: elec_num
integer*8 , intent(in) :: walk_num
double precision , intent(in) :: coord(elec_num,3,walk_num)
double precision , intent(in) :: coord(elec_num,walk_num,3)
double precision , intent(out) :: ee_distance(elec_num,elec_num,walk_num)
integer*8 :: k
integer*8 :: k, i, j
double precision :: x, y, z
info = QMCKL_SUCCESS
@ -1030,8 +1014,8 @@ integer function qmckl_compute_ee_distance_f(context, elec_num, walk_num, coord,
do k=1,walk_num
info = qmckl_distance(context, 'T', 'T', elec_num, elec_num, &
coord(1,1,k), elec_num, &
coord(1,1,k), elec_num, &
coord(1,k,1), elec_num * walk_num, &
coord(1,k,1), elec_num * walk_num, &
ee_distance(1,1,k), elec_num)
if (info /= QMCKL_SUCCESS) then
exit
@ -1075,7 +1059,7 @@ qmckl_exit_code qmckl_compute_ee_distance (
#+end_src
*** Test
#+begin_src python :results output :exports none
import numpy as np
@ -1210,7 +1194,7 @@ qmckl_exit_code qmckl_provide_ee_distance_rescaled(qmckl_context context)
ctx->electron.num,
ctx->electron.rescale_factor_kappa_en,
ctx->electron.walk_num,
ctx->electron.coord_new,
ctx->electron.coord_new.data,
ctx->electron.ee_distance_rescaled);
if (rc != QMCKL_SUCCESS) {
return rc;
@ -1274,8 +1258,8 @@ integer function qmckl_compute_ee_distance_rescaled_f(context, elec_num, rescale
do k=1,walk_num
info = qmckl_distance_rescaled(context, 'T', 'T', elec_num, elec_num, &
coord(1,1,k), elec_num, &
coord(1,1,k), elec_num, &
coord(1,k,1), elec_num * walk_num, &
coord(1,k,1), elec_num * walk_num, &
ee_distance_rescaled(1,1,k), elec_num, rescale_factor_kappa_ee)
if (info /= QMCKL_SUCCESS) then
exit
@ -1376,12 +1360,12 @@ assert(fabs(ee_distance_rescaled[elec_num*elec_num+1]-0.9985724058042633) < 1.e-
#+end_src
** Electron-electron rescaled distance gradients and laplacian with respect to electron coords
** Electron-electron rescaled distance gradients and laplacian with respect to electron coords
The rescaled distances which is given as $R = (1 - \exp{-\kappa r})/\kappa$
The rescaled distances which is given as $R = (1 - \exp{-\kappa r})/\kappa$
needs to be perturbed with respect to the electorn coordinates.
This data is stored in the ~ee_distance_rescaled_deriv_e~ tensor. The
The first three elements of this three index tensor ~[4][num][num]~ gives the
This data is stored in the ~ee_distance_rescaled_deriv_e~ tensor. The
The first three elements of this three index tensor ~[4][num][num]~ gives the
derivatives in the x, y, and z directions $dx, dy, dz$ and the last index
gives the Laplacian $\partial x^2 + \partial y^2 + \partial z^2$.
@ -1414,7 +1398,7 @@ qmckl_exit_code qmckl_get_electron_ee_distance_rescaled_deriv_e(qmckl_context co
#+end_src
*** Provide :noexport:
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_provide_ee_distance_rescaled_deriv_e(qmckl_context context);
#+end_src
@ -1456,7 +1440,7 @@ qmckl_exit_code qmckl_provide_ee_distance_rescaled_deriv_e(qmckl_context context
ctx->electron.num,
ctx->electron.rescale_factor_kappa_en,
ctx->electron.walk_num,
ctx->electron.coord_new,
ctx->electron.coord_new.data,
ctx->electron.ee_distance_rescaled_deriv_e);
if (rc != QMCKL_SUCCESS) {
return rc;
@ -1604,7 +1588,7 @@ rc = qmckl_get_electron_ee_distance_rescaled_deriv_e(context, ee_distance_rescal
#+end_src
** Electron-electron potential
~ee_pot~ calculates the ~ee~ potential energy.
\[
@ -1641,7 +1625,7 @@ qmckl_exit_code qmckl_get_electron_ee_potential(qmckl_context context, double* c
return QMCKL_SUCCESS;
}
#+end_src
*** Provide :noexport:
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
@ -1770,7 +1754,7 @@ end function qmckl_compute_ee_potential_f
const int64_t elec_num,
const int64_t walk_num,
const double* ee_distance,
double* const ee_pot );
double* const ee_pot );
#+end_src
#+CALL: generate_c_interface(table=qmckl_ee_potential_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
@ -1897,8 +1881,8 @@ qmckl_exit_code qmckl_provide_en_distance(qmckl_context context)
ctx->electron.num,
ctx->nucleus.num,
ctx->electron.walk_num,
ctx->electron.coord_new,
ctx->nucleus.coord,
ctx->electron.coord_new.data,
ctx->nucleus.coord.data,
ctx->electron.en_distance);
if (rc != QMCKL_SUCCESS) {
return rc;
@ -1938,7 +1922,7 @@ integer function qmckl_compute_en_distance_f(context, elec_num, nucl_num, walk_n
integer*8 , intent(in) :: elec_num
integer*8 , intent(in) :: nucl_num
integer*8 , intent(in) :: walk_num
double precision , intent(in) :: elec_coord(elec_num,3,walk_num)
double precision , intent(in) :: elec_coord(elec_num,walk_num,3)
double precision , intent(in) :: nucl_coord(nucl_num,3)
double precision , intent(out) :: en_distance(elec_num,nucl_num,walk_num)
@ -1968,7 +1952,7 @@ integer function qmckl_compute_en_distance_f(context, elec_num, nucl_num, walk_n
do k=1,walk_num
info = qmckl_distance(context, 'T', 'T', elec_num, nucl_num, &
elec_coord(1,1,k), elec_num, &
elec_coord(1,k,1), elec_num * walk_num, &
nucl_coord, nucl_num, &
en_distance(1,1,k), elec_num)
if (info /= QMCKL_SUCCESS) then
@ -2005,7 +1989,7 @@ qmckl_exit_code qmckl_compute_en_distance (
integer (c_int64_t) , intent(in) , value :: elec_num
integer (c_int64_t) , intent(in) , value :: nucl_num
integer (c_int64_t) , intent(in) , value :: walk_num
real (c_double ) , intent(in) :: elec_coord(elec_num,3,walk_num)
real (c_double ) , intent(in) :: elec_coord(elec_num,walk_num,3)
real (c_double ) , intent(in) :: nucl_coord(elec_num,3)
real (c_double ) , intent(out) :: en_distance(elec_num,nucl_num,walk_num)
@ -2052,10 +2036,10 @@ assert(qmckl_electron_provided(context));
rc = qmckl_set_nucleus_num (context, nucl_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_charge (context, charge);
rc = qmckl_set_nucleus_charge (context, charge, nucl_num);
assert (rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_coord (context, 'T', nucl_coord);
rc = qmckl_set_nucleus_coord (context, 'T', nucl_coord, 3*nucl_num);
assert (rc == QMCKL_SUCCESS);
assert(qmckl_nucleus_provided(context));
@ -2088,7 +2072,7 @@ assert(fabs(en_distance[1][0][1] - 3.1804527583077356) < 1.e-12);
** Electron-nucleus rescaled distances
~en_distance_rescaled~ stores the matrix of the rescaled distances between
~en_distance_rescaled~ stores the matrix of the rescaled distances between
electrons and nuclei.
\[
@ -2103,6 +2087,7 @@ assert(fabs(en_distance[1][0][1] - 3.1804527583077356) < 1.e-12);
qmckl_exit_code qmckl_get_electron_en_distance_rescaled(qmckl_context context, double* distance_rescaled);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_electron_en_distance_rescaled(qmckl_context context, double* distance_rescaled)
{
@ -2173,8 +2158,8 @@ qmckl_exit_code qmckl_provide_en_distance_rescaled(qmckl_context context)
ctx->nucleus.num,
ctx->electron.rescale_factor_kappa_en,
ctx->electron.walk_num,
ctx->electron.coord_new,
ctx->nucleus.coord,
ctx->electron.coord_new.data,
ctx->nucleus.coord.data,
ctx->electron.en_distance_rescaled);
if (rc != QMCKL_SUCCESS) {
return rc;
@ -2217,7 +2202,7 @@ integer function qmckl_compute_en_distance_rescaled_f(context, elec_num, nucl_nu
integer*8 , intent(in) :: nucl_num
double precision , intent(in) :: rescale_factor_kappa_en
integer*8 , intent(in) :: walk_num
double precision , intent(in) :: elec_coord(elec_num,3,walk_num)
double precision , intent(in) :: elec_coord(elec_num,walk_num,3)
double precision , intent(in) :: nucl_coord(nucl_num,3)
double precision , intent(out) :: en_distance_rescaled(elec_num,nucl_num,walk_num)
@ -2253,7 +2238,7 @@ integer function qmckl_compute_en_distance_rescaled_f(context, elec_num, nucl_nu
do k=1,walk_num
info = qmckl_distance_rescaled(context, 'T', 'T', elec_num, nucl_num, &
elec_coord(1,1,k), elec_num, &
elec_coord(1,k,1), elec_num*walk_num, &
nucl_coord, nucl_num, &
en_distance_rescaled(1,1,k), elec_num, rescale_factor_kappa_en)
if (info /= QMCKL_SUCCESS) then
@ -2292,7 +2277,7 @@ qmckl_exit_code qmckl_compute_en_distance_rescaled (
integer (c_int64_t) , intent(in) , value :: nucl_num
real (c_double ) , intent(in) , value :: rescale_factor_kappa_en
integer (c_int64_t) , intent(in) , value :: walk_num
real (c_double ) , intent(in) :: elec_coord(elec_num,3,walk_num)
real (c_double ) , intent(in) :: elec_coord(elec_num,walk_num,3)
real (c_double ) , intent(in) :: nucl_coord(elec_num,3)
real (c_double ) , intent(out) :: en_distance_rescaled(elec_num,nucl_num,walk_num)
@ -2341,10 +2326,10 @@ assert(qmckl_electron_provided(context));
rc = qmckl_set_nucleus_num (context, nucl_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_charge (context, charge);
rc = qmckl_set_nucleus_charge (context, charge, nucl_num);
assert (rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_coord (context, 'T', nucl_coord);
rc = qmckl_set_nucleus_coord (context, 'T', nucl_coord, 3*nucl_num);
assert (rc == QMCKL_SUCCESS);
assert(qmckl_nucleus_provided(context));
@ -2378,10 +2363,10 @@ assert(fabs(en_distance_rescaled[1][0][1] - 0.9584331688679852) < 1.e-12);
** Electron-nucleus rescaled distance gradients and laplacian with respect to electron coords
The rescaled distances which is given as $R = (1 - \exp{-\kappa r})/\kappa$
The rescaled distances which is given as $R = (1 - \exp{-\kappa r})/\kappa$
needs to be perturbed with respect to the nuclear coordinates.
This data is stored in the ~en_distance_rescaled_deriv_e~ tensor. The
The first three elements of this three index tensor ~[4][nucl_num][elec_num]~ gives the
This data is stored in the ~en_distance_rescaled_deriv_e~ tensor. The
The first three elements of this three index tensor ~[4][nucl_num][elec_num]~ gives the
derivatives in the x, y, and z directions $dx, dy, dz$ and the last index
gives the Laplacian $\partial x^2 + \partial y^2 + \partial z^2$.
@ -2461,8 +2446,8 @@ qmckl_exit_code qmckl_provide_en_distance_rescaled_deriv_e(qmckl_context context
ctx->nucleus.num,
ctx->electron.rescale_factor_kappa_en,
ctx->electron.walk_num,
ctx->electron.coord_new,
ctx->nucleus.coord,
ctx->electron.coord_new.data,
ctx->nucleus.coord.data,
ctx->electron.en_distance_rescaled_deriv_e);
if (rc != QMCKL_SUCCESS) {
return rc;
@ -2506,7 +2491,7 @@ integer function qmckl_compute_en_distance_rescaled_deriv_e_f(context, elec_num,
integer*8 , intent(in) :: nucl_num
double precision , intent(in) :: rescale_factor_kappa_en
integer*8 , intent(in) :: walk_num
double precision , intent(in) :: elec_coord(elec_num,3,walk_num)
double precision , intent(in) :: elec_coord(elec_num,walk_num,3)
double precision , intent(in) :: nucl_coord(nucl_num,3)
double precision , intent(out) :: en_distance_rescaled_deriv_e(elec_num,nucl_num,walk_num)
@ -2542,7 +2527,7 @@ integer function qmckl_compute_en_distance_rescaled_deriv_e_f(context, elec_num,
do k=1,walk_num
info = qmckl_distance_rescaled_deriv_e(context, 'T', 'T', elec_num, nucl_num, &
elec_coord(1,1,k), elec_num, &
elec_coord(1,k,1), elec_num*walk_num, &
nucl_coord, nucl_num, &
en_distance_rescaled_deriv_e(1,1,k), elec_num, rescale_factor_kappa_en)
if (info /= QMCKL_SUCCESS) then
@ -2581,7 +2566,7 @@ qmckl_exit_code qmckl_compute_en_distance_rescaled_deriv_e (
integer (c_int64_t) , intent(in) , value :: nucl_num
real (c_double ) , intent(in) , value :: rescale_factor_kappa_en
integer (c_int64_t) , intent(in) , value :: walk_num
real (c_double ) , intent(in) :: elec_coord(elec_num,3,walk_num)
real (c_double ) , intent(in) :: elec_coord(elec_num,walk_num,3)
real (c_double ) , intent(in) :: nucl_coord(elec_num,3)
real (c_double ) , intent(out) :: en_distance_rescaled_deriv_e(elec_num,nucl_num,walk_num)
@ -2593,7 +2578,7 @@ qmckl_exit_code qmckl_compute_en_distance_rescaled_deriv_e (
#+end_src
*** Test
#+begin_src python :results output :exports none
import numpy as np
@ -2611,10 +2596,10 @@ assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_rescale_factor (context, nucl_rescale_factor_kappa);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_charge (context, charge);
rc = qmckl_set_nucleus_charge (context, charge, nucl_num);
assert (rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_coord (context, 'T', nucl_coord);
rc = qmckl_set_nucleus_coord (context, 'T', nucl_coord, 3*nucl_num);
assert (rc == QMCKL_SUCCESS);
assert(qmckl_nucleus_provided(context));
@ -2656,7 +2641,7 @@ assert (rc == QMCKL_SUCCESS);
where \(\mathcal{V}_{en}\) is the ~en~ potential, \[r_{iA}\] the ~en~
distance and \[Z_A\] is the nuclear charge.
*** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
@ -2684,7 +2669,7 @@ qmckl_exit_code qmckl_get_electron_en_potential(qmckl_context context, double* c
return QMCKL_SUCCESS;
}
#+end_src
*** Provide :noexport:
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
@ -2732,7 +2717,7 @@ qmckl_exit_code qmckl_provide_en_potential(qmckl_context context)
ctx->electron.num,
ctx->nucleus.num,
ctx->electron.walk_num,
ctx->nucleus.charge,
ctx->nucleus.charge.data,
ctx->electron.en_distance,
ctx->electron.en_pot);
if (rc != QMCKL_SUCCESS) {
@ -2822,7 +2807,7 @@ end function qmckl_compute_en_potential_f
const int64_t walk_num,
const double* charge,
const double* en_distance,
double* const en_pot );
double* const en_pot );
#+end_src
#+CALL: generate_c_interface(table=qmckl_en_potential_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
@ -2858,7 +2843,7 @@ double en_pot[walk_num];
rc = qmckl_get_electron_en_potential(context, &(en_pot[0]));
assert (rc == QMCKL_SUCCESS);
#+end_src
** Generate initial coordinates
*** Compute :noexport:
@ -2872,11 +2857,11 @@ subroutine draw_init_points
integer :: iwalk
logical, allocatable :: do_elec(:)
integer :: acc_num
real, allocatable :: xmin(:,:)
integer :: i, j, k, l, kk
real :: norm
allocate (do_elec(elec_num), xmin(3,elec_num))
xmin = -huge(1.)
@ -2890,7 +2875,7 @@ subroutine draw_init_points
call rinfo( irp_here, 'Norm : ', norm )
call rinfo( irp_here, 'mo_scale: ' , mo_scale )
mo_coef_transp = mo_coef_transp/norm
double precision :: qmc_ranf
real :: mo_max
do i=1,elec_alpha_num
@ -2907,7 +2892,7 @@ subroutine draw_init_points
xmin(2,i) = nucl_coord(l,2)
xmin(3,i) = nucl_coord(l,3)
enddo
call iinfo(irp_here, 'Det num = ', det_num )
do k=1,elec_beta_num
i = k+elec_alpha_num
@ -2924,10 +2909,10 @@ subroutine draw_init_points
xmin(2,i) = nucl_coord(l,2)
xmin(3,i) = nucl_coord(l,3)
enddo
call rinfo( irp_here, 'time step =', time_step )
do iwalk=1,walk_num
print *, 'Generating initial positions for walker', iwalk
print *, 'Generating initial positions for walker', iwalk
acc_num = 0
do_elec = .True.
integer :: iter
@ -2949,7 +2934,7 @@ subroutine draw_init_points
TOUCH elec_coord
re_compute = minval(nucl_elec_dist(1:nucl_num,1:elec_num))
enddo
do i=1,elec_alpha_num
if (do_elec(i)) then
if ( mo_value_transp(i,i)**2 >= qmc_ranf()) then
@ -2958,7 +2943,7 @@ subroutine draw_init_points
endif
endif
enddo
do i=1,elec_beta_num
if (do_elec(i+elec_alpha_num)) then
if ( mo_value_transp(i,i+elec_alpha_num)**2 >= qmc_ranf()) then
@ -2967,9 +2952,9 @@ subroutine draw_init_points
endif
endif
enddo
enddo
do l=1,3
do i=1,elec_num+1
elec_coord_full(i,l,iwalk) = elec_coord(i,l)
@ -2982,7 +2967,7 @@ subroutine draw_init_points
endif
SOFT_TOUCH elec_coord elec_coord_full
deallocate (do_elec, xmin)
end
# end_src

64
org/qmckl_jastrow.org
View File

@ -1,4 +1,5 @@
#+TITLE: Jastrow Factor
#+SETUPFILE: ../tools/theme.setup
#+INCLUDE: ../tools/lib.org
@ -103,26 +104,25 @@ int main() {
computed data:
| Variable | Type | In/Out | Description |
|-------------------------------+---------------------------------------------------------------------+---------------------------------------------------------------------------------------------------------+---------------------------------|
| ~dim_cord_vect~ | ~int64_t~ | Number of unique C coefficients | |
| ~dim_cord_vect_date~ | ~uint64_t~ | Number of unique C coefficients | |
| ~asymp_jasb~ | ~double[2]~ | Asymptotic component | |
| ~asymp_jasb_date~ | ~uint64_t~ | Asymptotic component | |
| ~cord_vect_full~ | ~double[dim_cord_vect][nucl_num]~ | vector of non-zero coefficients | |
| ~cord_vect_full_date~ | ~uint64_t~ | Keep track of changes here | |
| ~lkpm_combined_index~ | ~int64_t[4][dim_cord_vect]~ | Transform l,k,p, and m into consecutive indices | |
| ~lkpm_combined_index_date~ | ~uint64_t~ | Transform l,k,p, and m into consecutive indices | |
| ~tmp_c~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]~ | out | vector of non-zero coefficients |
| ~dtmp_c~ | ~double[walk_num][elec_num][4][nucl_num][0:cord_num][0:cord_num-1]~ | vector of non-zero coefficients | |
| ~een_rescaled_e~ | ~double[walk_num][walk_num][elec_num][elec_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord | |
| ~een_rescaled_e_date~ | ~uint64_t~ | Keep track of the date of creation | |
| ~een_rescaled_n~ | ~double[walk_num][elec_num][nucl_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord | |
| ~een_rescaled_n_date~ | ~uint64_t~ | Keep track of the date of creation | |
| ~een_rescaled_e_deriv_e~ | ~double[walk_num][elec_num][4][elec_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons | |
| ~een_rescaled_e_deriv_e_date~ | ~uint64_t~ | Keep track of the date of creation | |
| ~een_rescaled_n_deriv_e~ | ~double[walk_num][elec_num][4][nucl_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons | |
| ~een_rescaled_n_deriv_e_date~ | ~uint64_t~ | Keep track of the date of creation | |
| Variable | Type | In/Out | Description |
|----------------------------+---------------------------------------------------------------------+-------------------------------------------------+---------------------------------|
| ~dim_cord_vect~ | ~int64_t~ | Number of unique C coefficients | |
| ~dim_cord_vect_date~ | ~uint64_t~ | Number of unique C coefficients | |
| ~asymp_jasb~ | ~double[2]~ | Asymptotic component | |
| ~asymp_jasb_date~ | ~uint64_t~ | Asymptotic component | |
| ~cord_vect_full~ | ~double[dim_cord_vect][nucl_num]~ | vector of non-zero coefficients | |
| ~cord_vect_full_date~ | ~uint64_t~ | Keep track of changes here | |
| ~lkpm_combined_index~ | ~int64_t[4][dim_cord_vect]~ | Transform l,k,p, and m into consecutive indices | |
| ~lkpm_combined_index_date~ | ~uint64_t~ | Transform l,k,p, and m into consecutive indices | |
| ~tmp_c~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]~ | out | vector of non-zero coefficients |
| ~dtmp_c~ | ~double[walk_num][elec_num][4][nucl_num][0:cord_num][0:cord_num-1]~ | vector of non-zero coefficients | |
| ~een_rescaled_n~ | ~double[walk_num][elec_num][nucl_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord | |
| ~een_rescaled_n_date~ | ~uint64_t~ | Keep track of the date of creation | |
| ~een_rescaled_e_deriv_e~ | ~double[walk_num][elec_num][4][elec_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons | |
| ~een_rescaled_e_deriv_e_date~ | ~uint64_t~ | Keep track of the date of creation | |
| ~een_rescaled_n_deriv_e~ | ~double[walk_num][elec_num][4][nucl_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons | |
| ~een_rescaled_n_deriv_e_date~ | ~uint64_t~ | Keep track of the date of creation | |
For H2O we have the following data:
@ -1093,7 +1093,7 @@ assert(rc == QMCKL_SUCCESS);
double elec_coord2[walk_num*3*elec_num];
rc = qmckl_get_electron_coord (context, 'N', elec_coord2);
rc = qmckl_get_electron_coord (context, 'N', elec_coord2, walk_num*3*elec_num);
assert(rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i<3*elec_num ; ++i) {
assert( elec_coord[i] == elec_coord2[i] );
@ -1131,15 +1131,15 @@ assert(k == nucl_rescale_factor_kappa);
double nucl_coord2[3*nucl_num];
rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2);
rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2, 3*nucl_num);
assert(rc == QMCKL_NOT_PROVIDED);
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]));
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), 3*nucl_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_nucleus_provided(context));
rc = qmckl_get_nucleus_coord (context, 'N', nucl_coord2);
rc = qmckl_get_nucleus_coord (context, 'N', nucl_coord2, nucl_num*3);
assert(rc == QMCKL_SUCCESS);
for (int64_t k=0 ; k<3 ; ++k) {
for (int64_t i=0 ; i<nucl_num ; ++i) {
@ -1147,7 +1147,7 @@ for (int64_t k=0 ; k<3 ; ++k) {
}
}
rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2);
rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2, nucl_num*3);
assert(rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i<3*nucl_num ; ++i) {
assert( nucl_coord[i] == nucl_coord2[i] );
@ -1155,13 +1155,13 @@ for (int64_t i=0 ; i<3*nucl_num ; ++i) {
double nucl_charge2[nucl_num];
rc = qmckl_get_nucleus_charge(context, nucl_charge2);
rc = qmckl_get_nucleus_charge(context, nucl_charge2, nucl_num);
assert(rc == QMCKL_NOT_PROVIDED);
rc = qmckl_set_nucleus_charge(context, nucl_charge);
rc = qmckl_set_nucleus_charge(context, nucl_charge, nucl_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_get_nucleus_charge(context, nucl_charge2);
rc = qmckl_get_nucleus_charge(context, nucl_charge2, nucl_num);
assert(rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i<nucl_num ; ++i) {
assert( nucl_charge[i] == nucl_charge2[i] );
@ -1169,7 +1169,7 @@ for (int64_t i=0 ; i<nucl_num ; ++i) {
assert(qmckl_nucleus_provided(context));
#+end_src
* Computation
The computed data is stored in the context so that it can be reused
@ -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,8 +3924,8 @@ qmckl_exit_code qmckl_provide_een_rescaled_n_deriv_e(qmckl_context context)
ctx->nucleus.num,
ctx->jastrow.cord_num,
ctx->electron.rescale_factor_kappa_en,
ctx->electron.coord_new,
ctx->nucleus.coord,
ctx->electron.coord_new.data,
ctx->nucleus.coord.data,
ctx->electron.en_distance,
ctx->jastrow.een_rescaled_n,
ctx->jastrow.een_rescaled_n_deriv_e);

4
org/qmckl_local_energy.org
View File

@ -576,10 +576,10 @@ assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_num (context, nucl_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]));
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), nucl_num*3);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_charge(context, nucl_charge);
rc = qmckl_set_nucleus_charge(context, nucl_charge, nucl_num);
assert(rc == QMCKL_SUCCESS);
assert(qmckl_nucleus_provided(context));

4
org/qmckl_mo.org
View File

@ -720,10 +720,10 @@ assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_num (context, nucl_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]));
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), nucl_num*3);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_charge(context, nucl_charge);
rc = qmckl_set_nucleus_charge(context, nucl_charge, nucl_num);
assert(rc == QMCKL_SUCCESS);
assert(qmckl_nucleus_provided(context));

534
org/qmckl_nucleus.org
View File

@ -67,42 +67,42 @@ int main() {
The following data stored in the context:
|------------------------+----------------+-------------------------------------------|
| ~uninitialized~ | int32_t | Keeps bit set for uninitialized data |
| ~num~ | int64_t | Total number of nuclei |
| ~provided~ | bool | If true, ~nucleus~ is valid |
| ~charge~ | double[num] | Nuclear charges |
| ~coord~ | double[3][num] | Nuclear coordinates, in transposed format |
| ~coord_date~ | int64_t | Nuclear coordinates, date if modified |
| ~rescale_factor_kappa~ | double | The distance scaling factor |
|------------------------+--------------+-------------------------------------------|
| ~uninitialized~ | int32_t | Keeps bit set for uninitialized data |
| ~num~ | int64_t | Total number of nuclei |
| ~provided~ | bool | If true, ~nucleus~ is valid |
| ~charge~ | qmckl_vector | Nuclear charges |
| ~coord~ | qmckl_matrix | Nuclear coordinates, in transposed format |
| ~coord_date~ | int64_t | Nuclear coordinates, date if modified |
| ~rescale_factor_kappa~ | double | The distance scaling factor |
Computed data:
|-----------------------------+------------------+------------------------------------------------------------|
| ~nn_distance~ | double[num][num] | Nucleus-nucleus distances |
| ~nn_distance_date~ | int64_t | Date when Nucleus-nucleus distances were computed |
| ~nn_distance_rescaled~ | double[num][num] | Nucleus-nucleus rescaled distances |
| ~nn_distance_rescaled_date~ | int64_t | Date when Nucleus-nucleus rescaled distances were computed |
| ~repulsion~ | double | Nuclear repulsion energy |
| ~repulsion_date~ | int64_t | Date when the nuclear repulsion energy was computed |
|-----------------------------+--------------+------------------------------------------------------------|
| ~nn_distance~ | qmckl_matrix | Nucleus-nucleus distances |
| ~nn_distance_date~ | int64_t | Date when Nucleus-nucleus distances were computed |
| ~nn_distance_rescaled~ | qmckl_matrix | Nucleus-nucleus rescaled distances |
| ~nn_distance_rescaled_date~ | int64_t | Date when Nucleus-nucleus rescaled distances were computed |
| ~repulsion~ | double | Nuclear repulsion energy |
| ~repulsion_date~ | int64_t | Date when the nuclear repulsion energy was computed |
** Data structure
#+begin_src c :comments org :tangle (eval h_private_type)
typedef struct qmckl_nucleus_struct {
int64_t num;
int64_t repulsion_date;
int64_t nn_distance_date;
int64_t nn_distance_rescaled_date;
int64_t coord_date;
double* coord;
double* charge;
double* nn_distance;
double* nn_distance_rescaled;
double repulsion;
double rescale_factor_kappa;
int32_t uninitialized;
bool provided;
int64_t num;
int64_t repulsion_date;
int64_t nn_distance_date;
int64_t nn_distance_rescaled_date;
int64_t coord_date;
qmckl_vector charge;
qmckl_matrix coord;
qmckl_matrix nn_distance;
qmckl_matrix nn_distance_rescaled;
double repulsion;
double rescale_factor_kappa;
int32_t uninitialized;
bool provided;
} qmckl_nucleus_struct;
#+end_src
@ -138,15 +138,6 @@ qmckl_exit_code qmckl_init_nucleus(qmckl_context context) {
#+end_src
** Access functions
#+begin_src c :comments org :tangle (eval h_func) :exports none
qmckl_exit_code qmckl_get_nucleus_num (const qmckl_context context, int64_t* const num);
qmckl_exit_code qmckl_get_nucleus_charge (const qmckl_context context, double* const charge);
qmckl_exit_code qmckl_get_nucleus_coord (const qmckl_context context, const char transp, double* const coord);
qmckl_exit_code qmckl_get_nucleus_rescale_factor (const qmckl_context context, double* const rescale_factor_kappa);
#+end_src
#+NAME:post
#+begin_src c :exports none
if ( (ctx->nucleus.uninitialized & mask) != 0) {
@ -154,6 +145,13 @@ if ( (ctx->nucleus.uninitialized & mask) != 0) {
}
#+end_src
#+begin_src c :comments org :tangle (eval h_func) :exports none
qmckl_exit_code
qmckl_get_nucleus_num(const qmckl_context context,
int64_t* const num);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_nucleus_num (const qmckl_context context, int64_t* const num) {
@ -187,10 +185,22 @@ qmckl_get_nucleus_num (const qmckl_context context, int64_t* const num) {
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src c :comments org :tangle (eval h_func) :exports none
qmckl_exit_code
qmckl_get_nucleus_charge (const qmckl_context context, double* const charge) {
qmckl_get_nucleus_charge(const qmckl_context context,
double* const charge,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_nucleus_charge (const qmckl_context context,
double* const charge,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_INVALID_CONTEXT;
@ -215,16 +225,31 @@ qmckl_get_nucleus_charge (const qmckl_context context, double* const charge) {
"nucleus data is not provided");
}
assert (ctx->nucleus.charge != NULL);
assert (ctx->nucleus.charge.data != NULL);
int64_t nucl_num = ctx->nucleus.num;
if (ctx->nucleus.num > size_max) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_nucleus_charge",
"Array too small");
}
memcpy(charge, ctx->nucleus.charge, nucl_num*sizeof(double));
qmckl_exit_code rc;
rc = qmckl_double_of_vector(context, ctx->nucleus.charge, charge, size_max);
return QMCKL_SUCCESS;
return rc;
}
#+end_src
#+begin_src c :comments org :tangle (eval h_func) :exports none
qmckl_exit_code
qmckl_get_nucleus_rescale_factor(const qmckl_context context,
double* const rescale_factor_kappa);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_nucleus_rescale_factor (const qmckl_context context,
double* const rescale_factor_kappa)
@ -250,10 +275,24 @@ qmckl_get_nucleus_rescale_factor (const qmckl_context context,
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src c :comments org :tangle (eval h_func) :exports none
qmckl_exit_code
qmckl_get_nucleus_coord (const qmckl_context context, const char transp, double* const coord) {
qmckl_get_nucleus_coord(const qmckl_context context,
const char transp,
double* const coord,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_nucleus_coord (const qmckl_context context,
const char transp,
double* const coord,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_INVALID_CONTEXT;
@ -285,25 +324,21 @@ qmckl_get_nucleus_coord (const qmckl_context context, const char transp, double*
"nucleus data is not provided");
}
int64_t nucl_num = ctx->nucleus.num;
assert (ctx->nucleus.coord.data != NULL);
assert (ctx->nucleus.coord != NULL);
qmckl_exit_code rc;
if (transp == 'N') {
qmckl_exit_code rc;
rc = qmckl_transpose(context, nucl_num, 3,
ctx->nucleus.coord, nucl_num,
coord, 3);
qmckl_matrix At = qmckl_matrix_alloc(context, 3, ctx->nucleus.coord.size[0]);
rc = qmckl_transpose(context, ctx->nucleus.coord, At);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_double_of_matrix(context, At, coord, size_max);
qmckl_matrix_free(context, At);
} else {
memcpy(coord, ctx->nucleus.coord, 3*nucl_num*sizeof(double));
rc = qmckl_double_of_matrix(context, ctx->nucleus.coord, coord, size_max);
}
return QMCKL_SUCCESS;
return rc;
}
#+end_src
@ -330,51 +365,6 @@ bool qmckl_nucleus_provided(const qmckl_context context) {
** Initialization functions
To set the data relative to the nuclei in the context, the
following functions need to be called.
#+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code qmckl_set_nucleus_num (qmckl_context context, const int64_t num);
qmckl_exit_code qmckl_set_nucleus_charge (qmckl_context context, const double* charge);
qmckl_exit_code qmckl_set_nucleus_coord (qmckl_context context, const char transp, const double* coord);
qmckl_exit_code qmckl_set_nucleus_rescale_factor (qmckl_context context, const double rescale_factor_kappa);
#+end_src
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_set_nucleus_num(context, num) &
bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(in) , value :: num
end function
end interface
interface
integer(c_int32_t) function qmckl_set_nucleus_charge(context, charge) &
bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
real (c_double) , intent(in) :: charge(*)
end function
end interface
interface
integer(c_int32_t) function qmckl_set_nucleus_coord(context, transp, coord) &
bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
character(c_char) , intent(in) , value :: transp
real (c_double) , intent(in) :: coord(*)
end function
end interface
#+end_src
#+NAME:pre2
#+begin_src c :exports none
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
@ -392,12 +382,22 @@ ctx->nucleus.provided = (ctx->nucleus.uninitialized == 0);
return QMCKL_SUCCESS;
#+end_src
To set the data relative to the nuclei in the context, the
following functions need to be called.
To set the number of nuclei, use
#+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_set_nucleus_num(qmckl_context context,
const int64_t num);
#+end_src
Sets the number of nuclei.
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_set_nucleus_num(qmckl_context context, const int64_t num) {
qmckl_set_nucleus_num(qmckl_context context,
const int64_t num)
{
<<pre2>>
if (num <= 0) {
@ -415,11 +415,35 @@ qmckl_set_nucleus_num(qmckl_context context, const int64_t num) {
}
#+end_src
The following function sets the nuclear charges of all the atoms.
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_set_nucleus_num(context, num) &
bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(in) , value :: num
end function
end interface
#+end_src
#+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_set_nucleus_charge(qmckl_context context,
const double* charge,
const int64_t size_max);
#+end_src
Sets the nuclear charges of all the atoms.
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_set_nucleus_charge(qmckl_context context, const double* charge) {
qmckl_set_nucleus_charge(qmckl_context context,
const double* charge,
const int64_t size_max)
{
<<pre2>>
if (charge == NULL) {
@ -437,34 +461,126 @@ qmckl_set_nucleus_charge(qmckl_context context, const double* charge) {
rc = qmckl_get_nucleus_num(context, &num);
if (rc != QMCKL_SUCCESS) return rc;
if (ctx->nucleus.charge != NULL) {
qmckl_free(context, ctx->nucleus.charge);
ctx->nucleus.charge= NULL;
}
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = num*sizeof(double);
assert (ctx->nucleus.charge == NULL);
ctx->nucleus.charge = (double*) qmckl_malloc(context, mem_info);
if (ctx->nucleus.charge == NULL) {
if (num > size_max) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
QMCKL_INVALID_ARG_3,
"qmckl_set_nucleus_charge",
NULL);
"Array too small");
}
ctx->nucleus.charge= memcpy(ctx->nucleus.charge, charge, num*sizeof(double));
assert (ctx->nucleus.charge != NULL);
ctx->nucleus.charge = qmckl_vector_alloc(context, num);
rc = qmckl_vector_of_double(context, charge, num, &(ctx->nucleus.charge));
<<post2>>
}
#+end_src
The following function sets the rescale parameter for the nuclear distances.
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_set_nucleus_charge(context, charge, size_max) &
bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
real (c_double) , intent(in) :: charge(*)
integer (c_int64_t) , intent(in) , value :: size_max
end function
end interface
#+end_src
#+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_set_nucleus_coord(qmckl_context context,
const char transp,
const double* coord,
const int64_t size_max);
#+end_src
Sets the nuclear coordinates of all the atoms. The coordinates
are be given in atomic units.
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_set_nucleus_rescale_factor(qmckl_context context, const double rescale_factor_kappa) {
qmckl_set_nucleus_coord(qmckl_context context,
const char transp,
const double* coord,
const int64_t size_max)
{
<<pre2>>
qmckl_exit_code rc;
int32_t mask = 1 << 2;
const int64_t nucl_num = (int64_t) ctx->nucleus.num;
if (ctx->nucleus.coord.data != NULL) {
rc = qmckl_matrix_free(context, ctx->nucleus.coord);
if (rc != QMCKL_SUCCESS) return rc;
}
ctx->nucleus.coord = qmckl_matrix_alloc(context, nucl_num, 3);
if (ctx->nucleus.coord.data == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_set_nucleus_coord",
NULL);
}
if (size_max < 3*nucl_num) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_4,
"qmckl_set_nucleus_coord",
"Array too small");
}
if (transp == 'N') {
qmckl_matrix At;
At = qmckl_matrix_alloc(context, 3, nucl_num);
rc = qmckl_matrix_of_double(context, coord, 3*nucl_num, &At);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_transpose(context, At, ctx->nucleus.coord);
} else {
rc = qmckl_matrix_of_double(context, coord, nucl_num*3, &(ctx->nucleus.coord));
}
if (rc != QMCKL_SUCCESS) return rc;
<<post2>>
}
#+end_src
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_set_nucleus_coord(context, transp, coord, size_max) &
bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
character(c_char) , intent(in) , value :: transp
real (c_double) , intent(in) :: coord(*)
integer (c_int64_t) , intent(in) , value :: size_max
end function
end interface
#+end_src
#+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_set_nucleus_rescale_factor(qmckl_context context,
const double kappa);
#+end_src
Sets the rescale parameter for the nuclear distances.
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_set_nucleus_rescale_factor(qmckl_context context,
const double rescale_factor_kappa)
{
<<pre2>>
if (rescale_factor_kappa <= 0.0) {
@ -480,50 +596,17 @@ qmckl_set_nucleus_rescale_factor(qmckl_context context, const double rescale_fac
}
#+end_src
The following function sets the nuclear coordinates of all the
atoms. The coordinates should be given in atomic units.
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_set_nucleus_coord(qmckl_context context, const char transp, const double* coord) {
<<pre2>>
int64_t nucl_num = (int64_t) 0;
qmckl_exit_code rc;
int32_t mask = 1 << 2;
rc = qmckl_get_nucleus_num(context, &nucl_num);
if (rc != QMCKL_SUCCESS) return rc;
if (ctx->nucleus.coord != NULL) {
qmckl_free(context, ctx->nucleus.coord);
ctx->nucleus.coord = NULL;
}
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = 3*nucl_num*sizeof(double);
assert(ctx->nucleus.coord == NULL);
ctx->nucleus.coord = (double*) qmckl_malloc(context, mem_info);
if (coord == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_set_nucleus_coord",
NULL);
}
if (transp == 'N') {
rc = qmckl_transpose(context, 3, nucl_num,
coord, 3,
ctx->nucleus.coord, nucl_num);
if (rc != QMCKL_SUCCESS) return rc;
} else {
memcpy(ctx->nucleus.coord, coord, 3*nucl_num*sizeof(double));
}
<<post2>>
}
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_set_rescale_factor(context, kappa) &
bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
real (c_double) , intent(in) , value :: kappa
end function
end interface
#+end_src
** Test
@ -568,15 +651,15 @@ assert(k == nucl_rescale_factor_kappa);
double nucl_coord2[3*nucl_num];
rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2);
rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2, 3*nucl_num);
assert(rc == QMCKL_NOT_PROVIDED);
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]));
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), 3*nucl_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_nucleus_provided(context));
rc = qmckl_get_nucleus_coord (context, 'N', nucl_coord2);
rc = qmckl_get_nucleus_coord (context, 'N', nucl_coord2, 3*nucl_num);
assert(rc == QMCKL_SUCCESS);
for (size_t k=0 ; k<3 ; ++k) {
for (size_t i=0 ; i<nucl_num ; ++i) {
@ -584,7 +667,7 @@ for (size_t k=0 ; k<3 ; ++k) {
}
}
rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2);
rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2, 3*nucl_num);
assert(rc == QMCKL_SUCCESS);
for (size_t i=0 ; i<3*nucl_num ; ++i) {
assert( nucl_coord[i] == nucl_coord2[i] );
@ -592,13 +675,13 @@ for (size_t i=0 ; i<3*nucl_num ; ++i) {
double nucl_charge2[nucl_num];
rc = qmckl_get_nucleus_charge(context, nucl_charge2);
rc = qmckl_get_nucleus_charge(context, nucl_charge2, nucl_num);
assert(rc == QMCKL_NOT_PROVIDED);
rc = qmckl_set_nucleus_charge(context, nucl_charge);
rc = qmckl_set_nucleus_charge(context, nucl_charge, nucl_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_get_nucleus_charge(context, nucl_charge2);
rc = qmckl_get_nucleus_charge(context, nucl_charge2, nucl_num);
assert(rc == QMCKL_SUCCESS);
for (size_t i=0 ; i<nucl_num ; ++i) {
assert( nucl_charge[i] == nucl_charge2[i] );
@ -621,11 +704,17 @@ assert(qmckl_nucleus_provided(context));
*** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_nucleus_nn_distance(qmckl_context context, double* distance);
qmckl_exit_code
qmckl_get_nucleus_nn_distance(qmckl_context context,
double* distance,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_nucleus_nn_distance(qmckl_context context, double* distance)
qmckl_exit_code
qmckl_get_nucleus_nn_distance(qmckl_context context,
double* distance,
const int64_t size_max)
{
/* Check input parameters */
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
@ -638,22 +727,29 @@ qmckl_exit_code qmckl_get_nucleus_nn_distance(qmckl_context context, double* dis
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
size_t sze = ctx->nucleus.num * ctx->nucleus.num;
memcpy(distance, ctx->nucleus.nn_distance, sze * sizeof(double));
const int64_t sze = ctx->nucleus.num * ctx->nucleus.num;
if (sze > size_max) {
return qmckl_failwith(context,
QMCKL_INVALID_ARG_3,
"qmckl_get_nucleus_nn_distance",
"Array too small");
}
rc = qmckl_double_of_matrix(context, ctx->nucleus.nn_distance, distance, size_max);
return QMCKL_SUCCESS;
return rc;
}
#+end_src
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_get_nucleus_nn_distance(context, distance) &
integer(c_int32_t) function qmckl_get_nucleus_nn_distance(context, distance, size_max) &
bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
real (c_double ) , intent(out) :: distance(*)
integer (c_int64_t) , intent(in) , value :: size_max
end function
end interface
#+end_src
@ -678,26 +774,23 @@ qmckl_exit_code qmckl_provide_nn_distance(qmckl_context context)
if (!ctx->nucleus.provided) return QMCKL_NOT_PROVIDED;
/* Allocate array */
if (ctx->nucleus.nn_distance == NULL) {
if (ctx->nucleus.nn_distance.data == NULL) {
ctx->nucleus.nn_distance =
qmckl_matrix_alloc(context, ctx->nucleus.num, ctx->nucleus.num);
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->nucleus.num * ctx->nucleus.num * sizeof(double);
double* nn_distance = (double*) qmckl_malloc(context, mem_info);
if (nn_distance == NULL) {
if (ctx->nucleus.nn_distance.data == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_nn_distance",
NULL);
}
ctx->nucleus.nn_distance = nn_distance;
}
qmckl_exit_code rc =
qmckl_compute_nn_distance(context,
ctx->nucleus.num,
ctx->nucleus.coord,
ctx->nucleus.nn_distance);
ctx->nucleus.coord.data,
ctx->nucleus.nn_distance.data);
if (rc != QMCKL_SUCCESS) {
return rc;
}
@ -788,7 +881,7 @@ qmckl_exit_code qmckl_compute_nn_distance (
assert(qmckl_nucleus_provided(context));
double distance[nucl_num*nucl_num];
rc = qmckl_get_nucleus_nn_distance(context, distance);
rc = qmckl_get_nucleus_nn_distance(context, distance, nucl_num*nucl_num);
assert(distance[0] == 0.);
assert(distance[1] == distance[nucl_num]);
assert(fabs(distance[1]-2.070304721365169) < 1.e-12);
@ -800,11 +893,17 @@ assert(fabs(distance[1]-2.070304721365169) < 1.e-12);
*** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_nucleus_nn_distance_rescaled(qmckl_context context, double* distance_rescaled);
qmckl_exit_code
qmckl_get_nucleus_nn_distance_rescaled(qmckl_context context,
double* distance_rescaled,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_nucleus_nn_distance_rescaled(qmckl_context context, double* distance_rescaled)
qmckl_exit_code
qmckl_get_nucleus_nn_distance_rescaled(qmckl_context context,
double* distance_rescaled,
const int64_t size_max)
{
/* Check input parameters */
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
@ -817,13 +916,35 @@ qmckl_exit_code qmckl_get_nucleus_nn_distance_rescaled(qmckl_context context, do
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
size_t sze = ctx->nucleus.num * ctx->nucleus.num;
memcpy(distance_rescaled, ctx->nucleus.nn_distance_rescaled, sze * sizeof(double));
return QMCKL_SUCCESS;
const int64_t sze = ctx->nucleus.num * ctx->nucleus.num;
if (sze > size_max) {
return qmckl_failwith(context,
QMCKL_INVALID_ARG_3,
"qmckl_get_nucleus_nn_distance_rescaled",
"Array too small");
}
rc = qmckl_double_of_matrix(context,
ctx->nucleus.nn_distance_rescaled,
distance_rescaled,
size_max);
return rc;
}
#+end_src
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_get_nucleus_nn_distance_rescaled(context, distance_rescaled, size_max) &
bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
real (c_double ) , intent(out) :: distance_rescaled(*)
integer (c_int64_t) , intent(in) , value :: size_max
end function
end interface
#+end_src
*** Provide :noexport:
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
@ -844,27 +965,24 @@ qmckl_exit_code qmckl_provide_nn_distance_rescaled(qmckl_context context)
if (!ctx->nucleus.provided) return QMCKL_NOT_PROVIDED;
/* Allocate array */
if (ctx->nucleus.nn_distance_rescaled == NULL) {
if (ctx->nucleus.nn_distance_rescaled.data == NULL) {
ctx->nucleus.nn_distance_rescaled =
qmckl_matrix_alloc(context, ctx->nucleus.num, ctx->nucleus.num);
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->nucleus.num * ctx->nucleus.num * sizeof(double);
double* nn_distance_rescaled = (double*) qmckl_malloc(context, mem_info);
if (nn_distance_rescaled == NULL) {
if (ctx->nucleus.nn_distance_rescaled.data == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_nn_distance_rescaled",
NULL);
}
ctx->nucleus.nn_distance_rescaled = nn_distance_rescaled;
}
qmckl_exit_code rc =
qmckl_compute_nn_distance_rescaled(context,
ctx->nucleus.num,
ctx->nucleus.rescale_factor_kappa,
ctx->nucleus.coord,
ctx->nucleus.nn_distance_rescaled);
ctx->nucleus.coord.data,
ctx->nucleus.nn_distance_rescaled.data);
if (rc != QMCKL_SUCCESS) {
return rc;
}
@ -959,7 +1077,7 @@ qmckl_exit_code qmckl_compute_nn_distance_rescaled (
//assert(qmckl_nucleus_provided(context));
//
//double distance[nucl_num*nucl_num];
//rc = qmckl_get_nucleus_nn_distance(context, distance);
//rc = qmckl_get_nucleus_nn_distance(context, distance, nucl_num*nucl_num);
//assert(distance[0] == 0.);
//assert(distance[1] == distance[nucl_num]);
//assert(fabs(distance[1]-2.070304721365169) < 1.e-12);
@ -975,11 +1093,11 @@ qmckl_exit_code qmckl_compute_nn_distance_rescaled (
*** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_nucleus_repulsion(qmckl_context context, double* energy);
qmckl_exit_code qmckl_get_nucleus_repulsion(qmckl_context context, double* const energy);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_nucleus_repulsion(qmckl_context context, double* energy)
qmckl_exit_code qmckl_get_nucleus_repulsion(qmckl_context context, double* const energy)
{
/* Check input parameters */
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
@ -1037,8 +1155,8 @@ qmckl_exit_code qmckl_provide_nucleus_repulsion(qmckl_context context)
rc = qmckl_compute_nucleus_repulsion(context,
ctx->nucleus.num,
ctx->nucleus.charge,
ctx->nucleus.nn_distance,
ctx->nucleus.charge.data,
ctx->nucleus.nn_distance.data,
&(ctx->nucleus.repulsion));
if (rc != QMCKL_SUCCESS) {
return rc;

401
org/qmckl_point.org
View File

@ -17,11 +17,14 @@ walkers.
#ifndef QMCKL_POINT_HPT
#define QMCKL_POINT_HPT
#include <stdbool.h>
#include "qmckl_blas_private_type.h"
#+end_src
#+begin_src c :tangle (eval h_private_func)
#ifndef QMCKL_POINT_HPF
#define QMCKL_POINT_HPF
#include "qmckl_blas_private_type.h"
#include "qmckl_blas_private_func.h"
#+end_src
#+begin_src c :tangle (eval c_test) :noweb yes
@ -33,6 +36,8 @@ walkers.
#endif
#include "chbrclf.h"
#include "qmckl_blas_private_type.h"
#include "qmckl_blas_private_func.h"
int main() {
qmckl_context context;
@ -61,44 +66,41 @@ int main() {
#include "qmckl.h"
#include "qmckl_context_private_type.h"
#include "qmckl_memory_private_type.h"
#include "qmckl_memory_private_func.h"
#include "qmckl_blas_private_type.h"
#include "qmckl_point_private_func.h"
#include "qmckl_memory_private_func.h"
#include "qmckl_blas_private_func.h"
#+end_src
* Context
The following data stored in the context:
| Variable | Type | Description |
|-----------+---------------+------------------------|
| ~num~ | ~int64_t~ | Total number of points |
| ~coord_x~ | ~double[num]~ | X coordinates |
| ~coord_y~ | ~double[num]~ | Y coordinates |
| ~coord_z~ | ~double[num]~ | Z coordinates |
| Variable | Type | Description |
|----------+----------------+-------------------------------------------|
| ~num~ | ~int64_t~ | Total number of points |
| ~date~ | ~uint64_t~ | Last modification date of the coordinates |
| ~coord~ | ~qmckl_matrix~ | ~num~ \times 3 matrix3 |
We consider that 'transposed' and 'normal' storage follows the convention:
| | Normal | Transposed |
|---------+------------------+------------------|
| C | ~[point_num][3]~ | ~[3][point_num]~ |
| Fortran | ~(3,point_num)~ | ~(point_num,3)~ |
We consider that the matrix is stored 'transposed' and 'normal'
corresponds to the 3 \times ~num~ matrix.
** Data structure
#+begin_src c :comments org :tangle (eval h_private_type)
typedef struct qmckl_point_struct {
double* coord_x;
double* coord_y;
double* coord_z;
int64_t num;
int64_t num;
uint64_t date;
qmckl_matrix coord;
} qmckl_point_struct;
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func)
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_exit_code qmckl_init_point(qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c)
qmckl_exit_code qmckl_init_point(qmckl_context context) {
@ -109,16 +111,7 @@ qmckl_exit_code qmckl_init_point(qmckl_context context) {
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = sizeof(qmckl_point_struct);
ctx->point = (qmckl_point_struct*) qmckl_malloc(context, mem_info);
if (ctx->point == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_init_point",
NULL);
}
memset(ctx->point, 0, sizeof(qmckl_point_struct));
memset(&(ctx->point), 0, sizeof(qmckl_point_struct));
return QMCKL_SUCCESS;
}
@ -139,7 +132,7 @@ qmckl_exit_code qmckl_get_point_num (const qmckl_context context, int64_t* const
#+end_src
Returns the number of points stored in the context.
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_point_num (const qmckl_context context, int64_t* const num) {
@ -157,10 +150,9 @@ qmckl_get_point_num (const qmckl_context context, int64_t* const num) {
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
assert (ctx->point != NULL);
assert (ctx->point->num > (int64_t) 0);
,*num = ctx->point->num;
assert (ctx->point.num >= (int64_t) 0);
,*num = ctx->point.num;
return QMCKL_SUCCESS;
}
#+end_src
@ -172,16 +164,17 @@ interface
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(out) :: num
integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(out) :: num
end function
end interface
#+end_src
*** Point coordinates
#+begin_src c :comments org :tangle (eval h_func) :exports none
qmckl_exit_code qmckl_get_point(const qmckl_context context,
const char transp,
double* const coord,
const int64_t size_max);
#+end_src
@ -193,6 +186,7 @@ qmckl_exit_code qmckl_get_point(const qmckl_context context,
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_point(const qmckl_context context,
const char transp,
double* const coord,
const int64_t size_max)
{
@ -210,13 +204,11 @@ qmckl_get_point(const qmckl_context context,
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
assert (ctx->point != NULL);
int64_t point_num = ctx->point->num;
int64_t point_num = ctx->point.num;
if (point_num == 0) return QMCKL_NOT_PROVIDED;
assert (ctx->point->coord_x != NULL);
assert (ctx->point->coord_y != NULL);
assert (ctx->point->coord_z != NULL);
assert (ctx->point.coord.data != NULL);
if (size_max < 3*point_num) {
return qmckl_failwith( context,
@ -225,212 +217,120 @@ qmckl_get_point(const qmckl_context context,
"size_max too small");
}
double * ptr = coord;
for (int64_t i=0 ; i<point_num ; ++i) {
,*ptr = ctx->point->coord_x[i]; ++ptr;
,*ptr = ctx->point->coord_y[i]; ++ptr;
,*ptr = ctx->point->coord_z[i]; ++ptr;
qmckl_exit_code rc;
if (transp == 'N') {
qmckl_matrix At = qmckl_matrix_alloc( context, 3, point_num);
rc = qmckl_transpose( context, ctx->point.coord, At );
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_double_of_matrix( context, At, coord, size_max);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_matrix_free(context, At);
} else {
rc = qmckl_double_of_matrix( context, ctx->point.coord, coord, size_max);
}
return QMCKL_SUCCESS;
if (rc != QMCKL_SUCCESS) return rc;
return rc;
}
#+end_src
#+begin_src f90 :comments org :tangle (eval fh_func) :noweb yes
interface
integer(c_int32_t) function qmckl_get_point(context, coord, size_max) bind(C)
integer(c_int32_t) function qmckl_get_point(context, transp, coord, size_max) bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
character(c_char) , intent(in) , value :: transp
real (c_double ) , intent(out) :: coord(*)
integer (c_int64_t) , intent(in) :: size_max
end function
end interface
#+end_src
#+begin_src c :comments org :tangle (eval h_func) :exports none
qmckl_exit_code qmckl_get_point_xyz (const qmckl_context context,
double* const coord_x,
double* const coord_y,
double* const coord_z,
const int64_t size_max);
#+end_src
Returns the point coordinates in three different arrays, one for
each component x,y,z.
The pointers are assumed to point on a memory block of size
~size_max~ \ge ~point_num~.
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_point_xyz (const qmckl_context context,
double* const coord_x,
double* const coord_y,
double* const coord_z,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_INVALID_CONTEXT;
}
if (coord_x == NULL) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
"qmckl_get_point_coord_xyz",
"coord_x is a null pointer");
}
if (coord_y == NULL) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_point_coord_xyz",
"coord_y is a null pointer");
}
if (coord_z == NULL) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_4,
"qmckl_get_point_coord_xyz",
"coord_z is a null pointer");
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
assert (ctx->point != NULL);
int64_t point_num = ctx->point->num;
assert (ctx->point->coord_x != NULL);
assert (ctx->point->coord_y != NULL);
assert (ctx->point->coord_z != NULL);
if (size_max < point_num) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_5,
"qmckl_get_point_coord_xyz",
"size_max too small");
}
memcpy(coord_x, ctx->point->coord_x, point_num*sizeof(double));
memcpy(coord_y, ctx->point->coord_y, point_num*sizeof(double));
memcpy(coord_z, ctx->point->coord_z, point_num*sizeof(double));
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src f90 :comments org :tangle (eval fh_func) :noweb yes
interface
integer(c_int32_t) function qmckl_get_point_xyz(context, &
coord_x, coord_y, coord_z, size_max) bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
real (c_double ) , intent(out) :: coord_x(*)
real (c_double ) , intent(out) :: coord_y(*)
real (c_double ) , intent(out) :: coord_z(*)
integer (c_int64_t) , intent(in) :: size_max
end function
end interface
#+end_src
** Initialization functions
When the data is set in the context, if the arrays are large
enough, we overwrite the data contained in them.
#+NAME: check_alloc
#+begin_src c :exports none
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
assert (ctx->point != NULL);
if (ctx->point->num < num) {
if (ctx->point->coord_x != NULL) {
qmckl_free(context, ctx->point->coord_x);
ctx->point->coord_x = NULL;
}
if (ctx->point->coord_y != NULL) {
qmckl_free(context, ctx->point->coord_y);
ctx->point->coord_y = NULL;
}
if (ctx->point->coord_z != NULL) {
qmckl_free(context, ctx->point->coord_z);
ctx->point->coord_z = NULL;
}
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = num*sizeof(double);
ctx->point->coord_x = (double*) qmckl_malloc(context, mem_info);
if (ctx->point->coord_x == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_set_point",
NULL);
}
ctx->point->coord_y = (double*) qmckl_malloc(context, mem_info);
if (ctx->point->coord_y == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_set_point",
NULL);
}
ctx->point->coord_z = (double*) qmckl_malloc(context, mem_info);
if (ctx->point->coord_z == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_set_point",
NULL);
}
};
ctx->point->num = num;
#+end_src
To set the data relative to the points in the context, one of the
following functions need to be called.
following functions need to be called.
#+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code qmckl_set_point (qmckl_context context,
const char transp,
const double* coord,
const int64_t num);
#+end_src
Copy a sequence of (x,y,z) into the context.
Copy a sequence of ~num~ points $(x,y,z)$ into the context.
#+begin_src c :comments org :tangle (eval c) :noweb yes
qmckl_exit_code
qmckl_set_point (qmckl_context context,
const char transp,
const double* coord,
const int64_t num)
{
<<check_alloc>>
for (int64_t i=0 ; i<num ; ++i) {
ctx->point->coord_x[i] = coord[3*i ];
ctx->point->coord_y[i] = coord[3*i+1];
ctx->point->coord_z[i] = coord[3*i+2];
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
if (transp != 'N' && transp != 'T') {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
"qmckl_set_point",
"transp should be 'N' or 'T'");
}
if (coord == NULL) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_set_point",
"coord is a NULL pointer");
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
qmckl_exit_code rc;
if (ctx->point.num < num) {
if (ctx->point.coord.data != NULL) {
rc = qmckl_matrix_free(context, ctx->point.coord);
assert (rc == QMCKL_SUCCESS);
}
ctx->point.coord = qmckl_matrix_alloc(context, num, 3);
if (ctx->point.coord.data == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_set_point",
NULL);
}
};
ctx->point.num = num;
if (transp == 'T') {
memcpy(ctx->point.coord.data, coord, 3*num*sizeof(double));
} else {
for (int64_t i=0 ; 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];
}
}
/* Increment the date of the context */
ctx->date += 1UL;
ctx->point.date = ctx->date;
return QMCKL_SUCCESS;
}
@ -439,76 +339,38 @@ qmckl_set_point (qmckl_context context,
#+begin_src f90 :comments org :tangle (eval fh_func) :noweb yes
interface
integer(c_int32_t) function qmckl_set_point(context, &
coord_x, coord_y, coord_z, size_max) bind(C)
transp, coord, size_max) bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
real (c_double ) , intent(in) :: coord_x(*)
real (c_double ) , intent(in) :: coord_y(*)
real (c_double ) , intent(in) :: coord_z(*)
character(c_char) , intent(in) , value :: transp
real (c_double ) , intent(in) :: coord(*)
integer (c_int64_t) , intent(in) , value :: size_max
end function
end interface
#+end_src
#+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code qmckl_set_point_xyz (qmckl_context context,
const double* coord_x,
const double* coord_y,
const double* coord_z,
const int64_t num);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes
qmckl_exit_code
qmckl_set_point_xyz (qmckl_context context,
const double* coord_x,
const double* coord_y,
const double* coord_z,
const int64_t num)
{
<<check_alloc>>
memcpy(ctx->point->coord_x, coord_x, num*sizeof(double));
memcpy(ctx->point->coord_y, coord_y, num*sizeof(double));
memcpy(ctx->point->coord_z, coord_z, num*sizeof(double));
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src f90 :comments org :tangle (eval fh_func) :noweb yes
interface
integer(c_int32_t) function qmckl_set_point_xyz(context, &
coord_x, coord_y, coord_z, size_max) bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
real (c_double ) , intent(in) :: coord_x(*)
real (c_double ) , intent(in) :: coord_y(*)
real (c_double ) , intent(in) :: coord_z(*)
integer (c_int64_t) , intent(in) , value :: size_max
end function
end interface
#+end_src
** Test
#+begin_src c :tangle (eval c_test)
/* Reference input data */
int64_t point_num = chbrclf_elec_num;
double* coord = &(chbrclf_elec_coord[0][0][0]);
const double* coord = &(chbrclf_elec_coord[0][0][0]);
/* --- */
qmckl_exit_code rc;
double coord2[point_num*3];
double coord3[point_num*3];
rc = qmckl_set_point (context, coord, point_num);
rc = qmckl_get_point (context, 'N', coord2, (point_num*3));
assert(rc == QMCKL_NOT_PROVIDED);
rc = qmckl_set_point (context, 'N', coord, point_num);
assert(rc == QMCKL_SUCCESS);
int64_t n;
@ -516,25 +378,30 @@ rc = qmckl_get_point_num (context, &n);
assert(rc == QMCKL_SUCCESS);
assert(n == point_num);
double coord2[point_num*3];
double coord_x[point_num];
double coord_y[point_num];
double coord_z[point_num];
rc = qmckl_get_point_xyz (context, coord_x, coord_y, coord_z, point_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_get_point (context, coord2, (point_num*3));
rc = qmckl_get_point (context, 'N', coord2, (point_num*3));
assert(rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i<3*point_num ; ++i) {
assert( coord[i] == coord2[i] );
}
rc = qmckl_get_point (context, 'T', coord2, (point_num*3));
assert(rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i<point_num ; ++i) {
assert( coord[3*i+0] == coord_x[i] );
assert( coord[3*i+1] == coord_y[i] );
assert( coord[3*i+2] == coord_z[i] );
assert( coord[3*i+0] == coord2[i] );
assert( coord[3*i+1] == coord2[i+point_num] );
assert( coord[3*i+2] == coord2[i+point_num*2] );
}
rc = qmckl_set_point (context, 'T', coord2, point_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_get_point (context, 'N', coord3, (point_num*3));
assert(rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i<3*point_num ; ++i) {
assert( coord[i] == coord3[i] );
}
#+end_src

10
org/qmckl_trexio.org
View File

@ -208,8 +208,7 @@ qmckl_trexio_read_nucleus_X(qmckl_context context, trexio_t* const file)
trexio_string_of_error(rcio));
}
//rc = qmckl_set_nucleus_charge(context, nucl_charge, nucleus_num);
rc = qmckl_set_nucleus_charge(context, nucl_charge);
rc = qmckl_set_nucleus_charge(context, nucl_charge, nucleus_num);
qmckl_free(context, nucl_charge);
nucl_charge = NULL;
@ -249,8 +248,7 @@ qmckl_trexio_read_nucleus_X(qmckl_context context, trexio_t* const file)
trexio_string_of_error(rcio));
}
//rc = qmckl_set_nucleus_coord(context, 'N', nucl_coord, 3*nucleus_num);
rc = qmckl_set_nucleus_coord(context, 'N', nucl_coord);
rc = qmckl_set_nucleus_coord(context, 'N', nucl_coord, 3*nucleus_num);
qmckl_free(context, nucl_coord);
nucl_coord = NULL;
@ -1183,7 +1181,7 @@ assert (nucl_num == chbrclf_nucl_num);
printf("Nuclear charges\n");
double * charge = (double*) malloc (nucl_num * sizeof(double));
rc = qmckl_get_nucleus_charge(context, charge);
rc = qmckl_get_nucleus_charge(context, charge, nucl_num);
assert (rc == QMCKL_SUCCESS);
for (int i=0 ; i<nucl_num ; i++) {
assert (charge[i] == chbrclf_charge[i]);
@ -1193,7 +1191,7 @@ charge = NULL;
printf("Nuclear coordinates\n");
double * coord = (double*) malloc (nucl_num * 3 * sizeof(double));
rc = qmckl_get_nucleus_coord(context, 'T', coord);
rc = qmckl_get_nucleus_coord(context, 'T', coord, 3*nucl_num);
assert (rc == QMCKL_SUCCESS);
int k=0;
for (int j=0 ; j<3 ; ++j) {

229
org/qmckl_utils.org
View File

@ -1,229 +0,0 @@
#+TITLE: Utility functions
#+SETUPFILE: ../tools/theme.setup
#+INCLUDE: ../tools/lib.org
* Headers :noexport:
#+begin_src elisp :noexport :results none
(org-babel-lob-ingest "../tools/lib.org")
#+end_src
#+begin_src c :comments link :tangle (eval c_test) :noweb yes
#include "qmckl.h"
#include "assert.h"
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
int main() {
qmckl_context context;
context = qmckl_context_create();
#+end_src
* Matrix operations
** ~qmckl_transpose~
Transposes a matrix: $B_{ji} = A_{ij}$
#+NAME: qmckl_transpose_args
| qmckl_context | context | in | Global state |
| int64_t | m | in | Number of rows of the input matrix |
| int64_t | n | in | Number of columns of the input matrix |
| double | A[][lda] | in | Array containing the $m \times n$ matrix $A$ |
| int64_t | lda | in | Leading dimension of array ~A~ |
| double | B[][ldb] | out | Array containing the $n \times m$ matrix $B$ |
| int64_t | ldb | in | Leading dimension of array ~B~ |
*** Requirements
- ~context~ is not ~QMCKL_NULL_CONTEXT~
- ~m > 0~
- ~n > 0~
- ~lda >= m~
- ~ldb >= n~
- ~A~ is allocated with at least $m \times n \times 8$ bytes
- ~B~ is allocated with at least $n \times m \times 8$ bytes
*** C header
#+CALL: generate_c_header(table=qmckl_transpose_args,rettyp="qmckl_exit_code",fname="qmckl_transpose")
#+RESULTS:
#+begin_src c :tangle (eval h_func) :comments org
qmckl_exit_code qmckl_transpose (
const qmckl_context context,
const int64_t m,
const int64_t n,
const double* A,
const int64_t lda,
double* const B,
const int64_t ldb );
#+end_src
*** Source
#+begin_src f90 :tangle (eval f)
integer function qmckl_transpose_f(context, m, n, A, LDA, B, LDB) &
result(info)
use qmckl
implicit none
integer(qmckl_context) , intent(in) :: context
integer*8 , intent(in) :: m, n
integer*8 , intent(in) :: lda
real*8 , intent(in) :: A(lda,*)
integer*8 , intent(in) :: ldb
real*8 , intent(out) :: B(ldb,*)
integer*8 :: i,j
info = QMCKL_SUCCESS
if (context == QMCKL_NULL_CONTEXT) then
info = QMCKL_INVALID_CONTEXT
return
endif
if (m <= 0_8) then
info = QMCKL_INVALID_ARG_2
return
endif
if (n <= 0_8) then
info = QMCKL_INVALID_ARG_3
return
endif
if (LDA < m) then
info = QMCKL_INVALID_ARG_5
return
endif
if (LDB < n) then
info = QMCKL_INVALID_ARG_7
return
endif
do j=1,m
do i=1,n
B(i,j) = A(j,i)
end do
end do
end function qmckl_transpose_f
#+end_src
*** C interface :noexport:
#+CALL: generate_c_interface(table=qmckl_transpose_args,rettyp="qmckl_exit_code",fname="qmckl_transpose")
#+RESULTS:
#+begin_src f90 :tangle (eval f) :comments org :exports none
integer(c_int32_t) function qmckl_transpose &
(context, m, n, A, lda, B, ldb) &
bind(C) result(info)
use, intrinsic :: iso_c_binding
implicit none
integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(in) , value :: m
integer (c_int64_t) , intent(in) , value :: n
real (c_double ) , intent(in) :: A(lda,*)
integer (c_int64_t) , intent(in) , value :: lda
real (c_double ) , intent(out) :: B(ldb,*)
integer (c_int64_t) , intent(in) , value :: ldb
integer(c_int32_t), external :: qmckl_transpose_f
info = qmckl_transpose_f &
(context, m, n, A, lda, B, ldb)
end function qmckl_transpose
#+end_src
#+CALL: generate_f_interface(table=qmckl_transpose_args,rettyp="qmckl_exit_code",fname="qmckl_transpose")
#+RESULTS:
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_transpose &
(context, m, n, A, lda, B, ldb) &
bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(in) , value :: m
integer (c_int64_t) , intent(in) , value :: n
real (c_double ) , intent(in) :: A(lda,*)
integer (c_int64_t) , intent(in) , value :: lda
real (c_double ) , intent(out) :: B(ldb,*)
integer (c_int64_t) , intent(in) , value :: ldb
end function qmckl_transpose
end interface
#+end_src
*** Test :noexport:
#+begin_src f90 :tangle (eval f_test)
integer(qmckl_exit_code) function test_qmckl_transpose(context) bind(C)
use qmckl
implicit none
integer(qmckl_context), intent(in), value :: context
double precision, allocatable :: A(:,:), B(:,:)
integer*8 :: m, n, LDA, LDB
integer*8 :: i,j
double precision :: x
m = 5
n = 6
LDA = m+3
LDB = n+1
allocate( A(LDA,n), B(LDB,m) )
A = 0.d0
B = 0.d0
do j=1,n
do i=1,m
A(i,j) = -10.d0 + dble(i+j)
end do
end do
test_qmckl_transpose = qmckl_transpose(context, m, n, A, LDA, B, LDB)
if (test_qmckl_transpose /= QMCKL_SUCCESS) return
test_qmckl_transpose = QMCKL_FAILURE
x = 0.d0
do j=1,n
do i=1,m
x = x + (A(i,j)-B(j,i))**2
end do
end do
if (dabs(x) <= 1.d-15) then
test_qmckl_transpose = QMCKL_SUCCESS
endif
deallocate(A,B)
end function test_qmckl_transpose
#+end_src
#+begin_src c :comments link :tangle (eval c_test)
qmckl_exit_code test_qmckl_transpose(qmckl_context context);
assert(QMCKL_SUCCESS == test_qmckl_transpose(context));
#+end_src
* End of files :noexport:
#+begin_src c :comments link :tangle (eval c_test)
assert (qmckl_context_destroy(context) == QMCKL_SUCCESS);
return 0;
}
#+end_src
# -*- mode: org -*-
# vim: syntax=c