TREX/qmckl
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2025-01-03 10:06:09 +01:00
Code Issues Projects Releases Wiki Activity

More compact error checking in Jastrow

Browse Source
This commit is contained in:
Anthony Scemama 2023-05-26 09:51:15 +02:00
parent e3f99d0030
commit dee0054c34
1 changed files with 359 additions and 471 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

626
org/qmckl_jastrow_champ.org
View File

@ -146,7 +146,7 @@ int main() {
| ~type_nucl_vector~ | ~int64_t[nucl_num]~ | IDs of types of Nuclei |
| ~a_vector~ | ~double[aord_num + 1][type_nucl_num]~ | a polynomial coefficients |
| ~b_vector~ | ~double[bord_num + 1]~ | b polynomial coefficients |
| ~c_vector~ | ~double[cord_num][type_nucl_num]~ | c polynomial coefficients |
| ~c_vector~ | ~double[dim_c_vector][type_nucl_num]~ | c polynomial coefficients |
Computed data:
@ -777,13 +777,16 @@ qmckl_set_jastrow_champ_c_vector(qmckl_context context,
}
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = dim_c_vector * type_nucl_num * sizeof(double);
mem_info.size = dim_c_vector*type_nucl_num * sizeof(double);
if ((size_t) size_max < mem_info.size/sizeof(double)) {
if (size_max < dim_c_vector*type_nucl_num) {
char msg[256];
sprintf(msg, "Array too small. Expected dim_c_vector*type_nucl_num = %ld",
dim_c_vector*type_nucl_num );
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_set_jastrow_champ_c_vector",
"Array too small. Expected dim_c_vector * type_nucl_num");
msg);
}
double* new_array = (double*) qmckl_malloc(context, mem_info);
@ -1641,7 +1644,6 @@ assert(qmckl_nucleus_provided(context));
compute. If it is the case, then the data is recomputed and the
current date is stored.
** Electron-electron component
*** Asymptotic component
@ -2815,11 +2817,21 @@ qmckl_exit_code qmckl_compute_jastrow_champ_factor_ee_gl_hpc(
memset(factor_ee_gl, 0, elec_num*4*walk_num*sizeof(double));
double kf[bord_num+1];
for (int k=0 ; k<=bord_num ; ++k) {
kf[k] = (double) k;
}
for (int nw = 0; nw < walk_num; ++nw) {
for (int j = 0; j < elec_num; ++j) {
const double* dxj = &ee_distance_rescaled_gl[4*elec_num*(j+nw*elec_num)];
const double* xj = &ee_distance_rescaled [ elec_num*(j+nw*elec_num)];
double * restrict factor_ee_gl_0 = &(factor_ee_gl[nw*elec_num*4]);
double * restrict factor_ee_gl_1 = factor_ee_gl_0 + elec_num;
double * restrict factor_ee_gl_2 = factor_ee_gl_1 + elec_num;
double * restrict factor_ee_gl_3 = factor_ee_gl_2 + elec_num;
for (int i = 0; i < elec_num; ++i) {
if (j == i) continue;
@ -2837,29 +2849,27 @@ qmckl_exit_code qmckl_compute_jastrow_champ_factor_ee_gl_hpc(
(i < up_num && j < up_num ) || (i >= up_num && j >= up_num) ?
0.5 * b_vector[0] * invdenom2 : b_vector[0] * invdenom2;
double * restrict factor_ee_gl_0 = &(factor_ee_gl[nw*elec_num*4]);
double * restrict factor_ee_gl_1 = factor_ee_gl_0 + elec_num;
double * restrict factor_ee_gl_2 = factor_ee_gl_1 + elec_num;
double * restrict factor_ee_gl_3 = factor_ee_gl_2 + elec_num;
factor_ee_gl_0[i] += f*dx[0];
factor_ee_gl_1[i] += f*dx[1];
factor_ee_gl_2[i] += f*dx[2];
factor_ee_gl_3[i] += f*(dx[3] - 2.*b_vector[1]*grad_c2*invdenom);
factor_ee_gl_3[i] += f*dx[3];
factor_ee_gl_3[i] -= f*grad_c2*invdenom*2.0 * b_vector[1];
double kf=2.0;
double x1 = x;
x = 1.0;
double xk[bord_num+1];
xk[0] = 1.0;
for (int k=1 ; k<= bord_num ; ++k) {
xk[k] = xk[k-1]*x;
}
for (int k=2 ; k<= bord_num ; ++k) {
f = b_vector[k] * kf * x;
factor_ee_gl_0[i] += f*x1*dx[0];
factor_ee_gl_1[i] += f*x1*dx[1];
factor_ee_gl_2[i] += f*x1*dx[2];
factor_ee_gl_3[i] += f*(x1*dx[3] + (kf-1.)*grad_c2);
x *= x1;
kf += 1.;
const double f1 = b_vector[k] * kf[k] * xk[k-2];
const double f2 = f1*xk[1];
factor_ee_gl_0[i] += f2*dx[0];
factor_ee_gl_1[i] += f2*dx[1];
factor_ee_gl_2[i] += f2*dx[2];
factor_ee_gl_3[i] += f2*dx[3];
factor_ee_gl_3[i] += f1*kf[k-1]*grad_c2;
}
}
}
@ -7094,6 +7104,142 @@ assert(fabs( 0.005359281880312882 - een_rescaled_n_gl[0][2][1][0][5]) < 1.e-12
calculation of the three-body jastrow ~factor_een~ and its derivative
~factor_een_gl~.
**** Compute dim_c_vector
:PROPERTIES:
:Name: qmckl_compute_dim_c_vector
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
Computes the dimension of the vector of parameters.
#+begin_src python :exports results
def compute(cord_num):
dim_c_vector = 0
for p in range(2,cord_num+1):
for k in range(p-1, -1, -1):
if k != 0:
lmax = p - k
else:
lmax = p - k - 2
for l in range(lmax, -1, -1):
if ( ((p - k - l) & 1)==1): continue
dim_c_vector += 1
return dim_c_vector
return [ ("$N_{ord}$", "Number of parameters"), ("","") ] + \
[ (i, compute(i)) for i in range(2,11) ]
#+end_src
#+RESULTS:
| $N_{ord}$ | Number of parameters |
| | |
| 2 | 2 |
| 3 | 6 |
| 4 | 13 |
| 5 | 23 |
| 6 | 37 |
| 7 | 55 |
| 8 | 78 |
| 9 | 106 |
| 10 | 140 |
#+NAME: qmckl_factor_dim_c_vector_args
| Variable | Type | In/Out | Description |
|----------------+-----------------+--------+------------------------------------|
| ~context~ | ~qmckl_context~ | in | Global state |
| ~cord_num~ | ~int64_t~ | in | Order of polynomials |
| ~dim_c_vector~ | ~int64_t~ | out | Number of parameters per atom type |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_dim_c_vector_f( &
context, cord_num, dim_c_vector) &
result(info)
use qmckl
implicit none
integer(qmckl_context), intent(in) :: context
integer*8 , intent(in) :: cord_num
integer*8 , intent(out) :: dim_c_vector
double precision :: x
integer*8 :: i, a, k, l, p, lmax
info = QMCKL_SUCCESS
if (context == QMCKL_NULL_CONTEXT) then
info = QMCKL_INVALID_CONTEXT
return
endif
if (cord_num < 0) then
info = QMCKL_INVALID_ARG_2
return
endif
dim_c_vector = 0
do p = 2, cord_num
do k = p - 1, 0, -1
if (k .ne. 0) then
lmax = p - k
else
lmax = p - k - 2
endif
do l = lmax, 0, -1
if (iand(p - k - l, 1_8) == 1) cycle
dim_c_vector = dim_c_vector + 1
end do
end do
end do
end function qmckl_compute_dim_c_vector_f
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_compute_dim_c_vector (
const qmckl_context context,
const int64_t cord_num,
int64_t* const dim_c_vector){
int lmax;
if (context == QMCKL_NULL_CONTEXT) {
return QMCKL_INVALID_CONTEXT;
}
if (cord_num < 0) {
return QMCKL_INVALID_ARG_2;
}
,*dim_c_vector = 0;
for (int p=2; p <= cord_num; ++p){
for (int k=p-1; k >= 0; --k) {
if (k != 0) {
lmax = p - k;
} else {
lmax = p - k - 2;
}
for (int l = lmax; l >= 0; --l) {
if ( ((p - k - l) & 1)==1) continue;
,*dim_c_vector=*dim_c_vector+1;
}
}
}
return QMCKL_SUCCESS;
}
#+end_src
# #+CALL: generate_c_header(table=qmckl_factor_dim_c_vector_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_compute_dim_c_vector (
const qmckl_context context,
const int64_t cord_num,
int64_t* const dim_c_vector );
#+end_src
**** Get
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes
@ -7412,109 +7558,6 @@ qmckl_exit_code qmckl_provide_dtmp_c(qmckl_context context)
}
#+end_src
**** Compute dim_c_vector
:PROPERTIES:
:Name: qmckl_compute_dim_c_vector
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
#+NAME: qmckl_factor_dim_c_vector_args
| Variable | Type | In/Out | Description |
|-----------------+-----------------+--------+-----------------------------------|
| ~context~ | ~qmckl_context~ | in | Global state |
| ~cord_num~ | ~int64_t~ | in | Order of polynomials |
| ~dim_c_vector~ | ~int64_t~ | out | dimension of c_vector_full table |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_dim_c_vector_f( &
context, cord_num, dim_c_vector) &
result(info)
use qmckl
implicit none
integer(qmckl_context), intent(in) :: context
integer*8 , intent(in) :: cord_num
integer*8 , intent(out) :: dim_c_vector
double precision :: x
integer*8 :: i, a, k, l, p, lmax
info = QMCKL_SUCCESS
if (context == QMCKL_NULL_CONTEXT) then
info = QMCKL_INVALID_CONTEXT
return
endif
if (cord_num < 0) then
info = QMCKL_INVALID_ARG_2
return
endif
dim_c_vector = 0
do p = 2, cord_num
do k = p - 1, 0, -1
if (k .ne. 0) then
lmax = p - k
else
lmax = p - k - 2
endif
do l = lmax, 0, -1
if (iand(p - k - l, 1_8) == 1) cycle
dim_c_vector = dim_c_vector + 1
end do
end do
end do
end function qmckl_compute_dim_c_vector_f
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes
qmckl_exit_code qmckl_compute_dim_c_vector (
const qmckl_context context,
const int64_t cord_num,
int64_t* const dim_c_vector){
int lmax;
if (context == QMCKL_NULL_CONTEXT) {
return QMCKL_INVALID_CONTEXT;
}
if (cord_num < 0) {
return QMCKL_INVALID_ARG_2;
}
*dim_c_vector = 0;
for (int p=2; p <= cord_num; ++p){
for (int k=p-1; k >= 0; --k) {
if (k != 0) {
lmax = p - k;
} else {
lmax = p - k - 2;
}
for (int l = lmax; l >= 0; --l) {
if ( ((p - k - l) & 1)==1) continue;
*dim_c_vector=*dim_c_vector+1;
}
}
}
return QMCKL_SUCCESS;
}
#+end_src
# #+CALL: generate_c_header(table=qmckl_factor_dim_c_vector_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_compute_dim_c_vector (
const qmckl_context context,
const int64_t cord_num,
int64_t* const dim_c_vector );
#+end_src
**** Compute c_vector_full
:PROPERTIES:
:Name: qmckl_compute_c_vector_full
@ -7552,26 +7595,11 @@ integer function qmckl_compute_c_vector_full_doc_f( &
info = QMCKL_SUCCESS
if (context == QMCKL_NULL_CONTEXT) then
info = QMCKL_INVALID_CONTEXT
return
endif
if (nucl_num <= 0) then
info = QMCKL_INVALID_ARG_2
return
endif
if (type_nucl_num <= 0) then
info = QMCKL_INVALID_ARG_4
return
endif
if (dim_c_vector < 0) then
info = QMCKL_INVALID_ARG_5
return
endif
if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT
if (nucl_num <= 0) info = QMCKL_INVALID_ARG_2
if (dim_c_vector < 0) info = QMCKL_INVALID_ARG_3
if (type_nucl_num <= 0) info = QMCKL_INVALID_ARG_4
if (info /= QMCKL_SUCCESS) return
do a = 1, nucl_num
c_vector_full(a,1:dim_c_vector) = c_vector(type_nucl_vector(a)+1,1:dim_c_vector)
@ -7616,21 +7644,13 @@ qmckl_exit_code qmckl_compute_c_vector_full_hpc (
const double* c_vector,
double* const c_vector_full ) {
if (context == QMCKL_NULL_CONTEXT) {
return QMCKL_INVALID_CONTEXT;
}
if (nucl_num <= 0) {
return QMCKL_INVALID_ARG_2;
}
if (type_nucl_num <= 0) {
return QMCKL_INVALID_ARG_4;
}
if (dim_c_vector < 0) {
return QMCKL_INVALID_ARG_5;
}
if (context == QMCKL_NULL_CONTEXT) return QMCKL_INVALID_CONTEXT;
if (nucl_num <= 0) return QMCKL_INVALID_ARG_2;
if (dim_c_vector < 0) return QMCKL_INVALID_ARG_3;
if (type_nucl_num <= 0) return QMCKL_INVALID_ARG_4;
if (type_nucl_vector == NULL) return QMCKL_INVALID_ARG_5;
if (c_vector == NULL) return QMCKL_INVALID_ARG_6;
if (c_vector_full == NULL) return QMCKL_INVALID_ARG_7;
for (int i=0; i < dim_c_vector; ++i) {
for (int a=0; a < nucl_num; ++a){
@ -7639,7 +7659,7 @@ qmckl_exit_code qmckl_compute_c_vector_full_hpc (
}
return QMCKL_SUCCESS;
}
}
#+end_src
@ -7705,7 +7725,7 @@ qmckl_exit_code qmckl_compute_c_vector_full (
#+NAME: qmckl_factor_lkpm_combined_index_args
| Variable | Type | In/Out | Description |
|-----------------------+-----------------------------+--------+-------------------------------|
|-----------------------+----------------------------+--------+-------------------------------|
| ~context~ | ~qmckl_context~ | in | Global state |
| ~cord_num~ | ~int64_t~ | in | Order of polynomials |
| ~dim_c_vector~ | ~int64_t~ | in | dimension of cord full table |
@ -7726,21 +7746,10 @@ integer function qmckl_compute_lkpm_combined_index_f( &
info = QMCKL_SUCCESS
if (context == QMCKL_NULL_CONTEXT) then
info = QMCKL_INVALID_CONTEXT
return
endif
if (cord_num < 0) then
info = QMCKL_INVALID_ARG_2
return
endif
if (dim_c_vector < 0) then
info = QMCKL_INVALID_ARG_3
return
endif
if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT
if (cord_num < 0) info = QMCKL_INVALID_ARG_2
if (dim_c_vector < 0) info = QMCKL_INVALID_ARG_3
if (info /= QMCKL_SUCCESS) return
kk = 0
do p = 2, cord_num
@ -7774,20 +7783,10 @@ qmckl_exit_code qmckl_compute_lkpm_combined_index (
int kk, lmax, m;
if (context == QMCKL_NULL_CONTEXT) {
return QMCKL_INVALID_CONTEXT;
}
if (context == QMCKL_NULL_CONTEXT) return QMCKL_INVALID_CONTEXT;
if (cord_num < 0) return QMCKL_INVALID_ARG_2;
if (dim_c_vector < 0) return QMCKL_INVALID_ARG_3;
if (cord_num < 0) {
return QMCKL_INVALID_ARG_2;
}
if (dim_c_vector < 0) {
return QMCKL_INVALID_ARG_3;
}
/*
*/
kk = 0;
for (int p = 2; p <= cord_num; ++p) {
for (int k=(p-1); k >= 0; --k) {
@ -7901,25 +7900,13 @@ integer function qmckl_compute_tmp_c_doc_f( &
info = QMCKL_SUCCESS
if (context == QMCKL_NULL_CONTEXT) then
info = QMCKL_INVALID_CONTEXT
return
endif
if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT
if (cord_num < 0) info = QMCKL_INVALID_ARG_2
if (elec_num <= 0) info = QMCKL_INVALID_ARG_3
if (nucl_num <= 0) info = QMCKL_INVALID_ARG_4
if (walk_num <= 0) info = QMCKL_INVALID_ARG_5
if (info /= QMCKL_SUCCESS) return
if (cord_num < 0) then
info = QMCKL_INVALID_ARG_2
return
endif
if (elec_num <= 0) then
info = QMCKL_INVALID_ARG_3
return
endif
if (nucl_num <= 0) then
info = QMCKL_INVALID_ARG_4
return
endif
M = elec_num
N = nucl_num*(cord_num + 1)
@ -7993,25 +7980,11 @@ qmckl_exit_code qmckl_compute_tmp_c_hpc (
const double* een_rescaled_n,
double* const tmp_c ) {
if (context == QMCKL_NULL_CONTEXT) {
return QMCKL_INVALID_CONTEXT;
}
if (cord_num < 0) {
return QMCKL_INVALID_ARG_2;
}
if (elec_num <= 0) {
return QMCKL_INVALID_ARG_3;
}
if (nucl_num <= 0) {
return QMCKL_INVALID_ARG_4;
}
if (walk_num <= 0) {
return QMCKL_INVALID_ARG_5;
}
if (context == QMCKL_NULL_CONTEXT) return QMCKL_INVALID_CONTEXT;
if (cord_num < 0) return QMCKL_INVALID_ARG_2;
if (elec_num <= 0) return QMCKL_INVALID_ARG_3;
if (nucl_num <= 0) return QMCKL_INVALID_ARG_4;
if (walk_num <= 0) return QMCKL_INVALID_ARG_5;
qmckl_exit_code info = QMCKL_SUCCESS;
@ -8104,7 +8077,7 @@ qmckl_exit_code qmckl_compute_tmp_c_hpc (const qmckl_context context,
#+NAME: qmckl_factor_dtmp_c_args
| Variable | Type | In/Out | Description |
|--------------------------+------------------------------------------------------------------+--------+-----------------------------------------------|
|---------------------+------------------------------------------------------------------+--------+-----------------------------------------------|
| ~context~ | ~qmckl_context~ | in | Global state |
| ~cord_num~ | ~int64_t~ | in | Order of polynomials |
| ~elec_num~ | ~int64_t~ | in | Number of electrons |
@ -8169,33 +8142,20 @@ integer function qmckl_compute_dtmp_c_doc_f( &
double precision :: alpha, beta
integer*8 :: M, N, K, LDA, LDB, LDC
info = QMCKL_SUCCESS
if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT
if (cord_num < 0) info = QMCKL_INVALID_ARG_2
if (elec_num <= 0) info = QMCKL_INVALID_ARG_3
if (nucl_num <= 0) info = QMCKL_INVALID_ARG_4
if (walk_num <= 0) info = QMCKL_INVALID_ARG_5
if (info /= QMCKL_SUCCESS) return
TransA = 'N'
TransB = 'N'
alpha = 1.0d0
beta = 0.0d0
info = QMCKL_SUCCESS
if (context == QMCKL_NULL_CONTEXT) then
info = QMCKL_INVALID_CONTEXT
return
endif
if (cord_num < 0) then
info = QMCKL_INVALID_ARG_2
return
endif
if (elec_num <= 0) then
info = QMCKL_INVALID_ARG_3
return
endif
if (nucl_num <= 0) then
info = QMCKL_INVALID_ARG_4
return
endif
M = 4*elec_num
N = nucl_num*(cord_num + 1)
K = elec_num
@ -8269,25 +8229,11 @@ qmckl_compute_dtmp_c_hpc (const qmckl_context context,
double* const dtmp_c )
{
if (context == QMCKL_NULL_CONTEXT) {
return QMCKL_INVALID_CONTEXT;
}
if (cord_num < 0) {
return QMCKL_INVALID_ARG_2;
}
if (elec_num <= 0) {
return QMCKL_INVALID_ARG_3;
}
if (nucl_num <= 0) {
return QMCKL_INVALID_ARG_4;
}
if (walk_num <= 0) {
return QMCKL_INVALID_ARG_5;
}
if (context == QMCKL_NULL_CONTEXT) return QMCKL_INVALID_CONTEXT;
if (cord_num < 0) return QMCKL_INVALID_ARG_2;
if (elec_num <= 0) return QMCKL_INVALID_ARG_3;
if (nucl_num <= 0) return QMCKL_INVALID_ARG_4;
if (walk_num <= 0) return QMCKL_INVALID_ARG_5;
qmckl_exit_code info = QMCKL_SUCCESS;
@ -8622,30 +8568,12 @@ integer function qmckl_compute_jastrow_champ_factor_een_naive_f( &
info = QMCKL_SUCCESS
if (context == QMCKL_NULL_CONTEXT) then
info = QMCKL_INVALID_CONTEXT
return
endif
if (walk_num <= 0) then
info = QMCKL_INVALID_ARG_2
return
endif
if (elec_num <= 0) then
info = QMCKL_INVALID_ARG_3
return
endif
if (nucl_num <= 0) then
info = QMCKL_INVALID_ARG_4
return
endif
if (cord_num < 0) then
info = QMCKL_INVALID_ARG_5
return
endif
if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT
if (walk_num <= 0) info = QMCKL_INVALID_ARG_2
if (elec_num <= 0) info = QMCKL_INVALID_ARG_3
if (nucl_num <= 0) info = QMCKL_INVALID_ARG_4
if (cord_num < 0) info = QMCKL_INVALID_ARG_5
if (info /= QMCKL_SUCCESS) return
factor_een = 0.0d0
@ -8789,30 +8717,12 @@ integer function qmckl_compute_jastrow_champ_factor_een_doc_f( &
info = QMCKL_SUCCESS
if (context == QMCKL_NULL_CONTEXT) then
info = QMCKL_INVALID_CONTEXT
return
endif
if (walk_num <= 0) then
info = QMCKL_INVALID_ARG_2
return
endif
if (elec_num <= 0) then
info = QMCKL_INVALID_ARG_3
return
endif
if (nucl_num <= 0) then
info = QMCKL_INVALID_ARG_4
return
endif
if (cord_num < 0) then
info = QMCKL_INVALID_ARG_5
return
endif
if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT
if (walk_num <= 0) info = QMCKL_INVALID_ARG_2
if (elec_num <= 0) info = QMCKL_INVALID_ARG_3
if (nucl_num <= 0) info = QMCKL_INVALID_ARG_4
if (cord_num < 0) info = QMCKL_INVALID_ARG_5
if (info /= QMCKL_SUCCESS) return
factor_een = 0.0d0
@ -9038,6 +8948,22 @@ qmckl_get_jastrow_champ_factor_een_gl(qmckl_context context,
}
#+end_src
***** Fortran interface
#+begin_src f90 :tangle (eval fh_func) :comments org
interface
integer(qmckl_exit_code) function qmckl_get_jastrow_champ_factor_een_gl (context, &
factor_een_gl, size_max) bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (qmckl_context) , intent(in), value :: context
integer(c_int64_t), intent(in), value :: size_max
double precision, intent(out) :: factor_een_gl(size_max)
end function qmckl_get_jastrow_champ_factor_een_gl
end interface
#+end_src
#
**** Provide :noexport:
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_provide_jastrow_champ_factor_een_gl(qmckl_context context);
@ -9195,30 +9121,12 @@ integer function qmckl_compute_jastrow_champ_factor_een_gl_naive_f( &
info = QMCKL_SUCCESS
if (context == QMCKL_NULL_CONTEXT) then
info = QMCKL_INVALID_CONTEXT
return
endif
if (walk_num <= 0) then
info = QMCKL_INVALID_ARG_2
return
endif
if (elec_num <= 0) then
info = QMCKL_INVALID_ARG_3
return
endif
if (nucl_num <= 0) then
info = QMCKL_INVALID_ARG_4
return
endif
if (cord_num < 0) then
info = QMCKL_INVALID_ARG_5
return
endif
if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT
if (walk_num <= 0) info = QMCKL_INVALID_ARG_2
if (elec_num <= 0) info = QMCKL_INVALID_ARG_3
if (nucl_num <= 0) info = QMCKL_INVALID_ARG_4
if (cord_num < 0) info = QMCKL_INVALID_ARG_5
if (info /= QMCKL_SUCCESS) return
factor_een_gl = 0.0d0
@ -9237,10 +9145,8 @@ integer function qmckl_compute_jastrow_champ_factor_een_gl_naive_f( &
daccu = 0.0d0
daccu2 = 0.0d0
do i = 1, elec_num
accu = accu + een_rescaled_e(k, i, j, nw) * &
een_rescaled_n(m, a, i, nw)
accu2 = accu2 + een_rescaled_e(k, i, j, nw) * &
een_rescaled_n(m + l, a, i, nw)
accu = accu + een_rescaled_e(k, i, j, nw) * een_rescaled_n(m, a, i, nw)
accu2 = accu2 + een_rescaled_e(k, i, j, nw) * een_rescaled_n(m + l, a, i, nw)
daccu(1:4) = daccu(1:4) + een_rescaled_e_gl(k, j, 1:4, i, nw) * &
een_rescaled_n(m, a, i, nw)
daccu2(1:4) = daccu2(1:4) + een_rescaled_e_gl(k, j, 1:4, i, nw) * &
@ -9392,30 +9298,12 @@ integer function qmckl_compute_jastrow_champ_factor_een_gl_doc_f( &
info = QMCKL_SUCCESS
if (context == QMCKL_NULL_CONTEXT) then
info = QMCKL_INVALID_CONTEXT
return
endif
if (walk_num <= 0) then
info = QMCKL_INVALID_ARG_2
return
endif
if (elec_num <= 0) then
info = QMCKL_INVALID_ARG_3
return
endif
if (nucl_num <= 0) then
info = QMCKL_INVALID_ARG_4
return
endif
if (cord_num < 0) then
info = QMCKL_INVALID_ARG_5
return
endif
if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT
if (walk_num <= 0) info = QMCKL_INVALID_ARG_2
if (elec_num <= 0) info = QMCKL_INVALID_ARG_3
if (nucl_num <= 0) info = QMCKL_INVALID_ARG_4
if (cord_num < 0) info = QMCKL_INVALID_ARG_5
if (info /= QMCKL_SUCCESS) return
factor_een_gl = 0.0d0