mirror of
https://github.com/TREX-CoE/qmckl.git
synced 2024-08-16 18:38:28 +02:00
qmckl_compute_een_rescaled_e_hpc (c version) working
This commit is contained in:
parent
586eb92801
commit
b79a23897d
@ -3241,7 +3241,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_e(qmckl_context context)
|
|||||||
| ~een_rescaled_e~ | ~double[walk_num][0:cord_num][elec_num][elec_num]~ | out | Electron-electron rescaled distances |
|
| ~een_rescaled_e~ | ~double[walk_num][0:cord_num][elec_num][elec_num]~ | out | Electron-electron rescaled distances |
|
||||||
|
|
||||||
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
|
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
|
||||||
integer function qmckl_compute_een_rescaled_e_f( &
|
integer function qmckl_compute_een_rescaled_e_doc_f( &
|
||||||
context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, &
|
context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, &
|
||||||
ee_distance, een_rescaled_e) &
|
ee_distance, een_rescaled_e) &
|
||||||
result(info)
|
result(info)
|
||||||
@ -3260,7 +3260,6 @@ integer function qmckl_compute_een_rescaled_e_f( &
|
|||||||
|
|
||||||
allocate(een_rescaled_e_ij(elec_num * (elec_num - 1) / 2, cord_num + 1))
|
allocate(een_rescaled_e_ij(elec_num * (elec_num - 1) / 2, cord_num + 1))
|
||||||
|
|
||||||
|
|
||||||
info = QMCKL_SUCCESS
|
info = QMCKL_SUCCESS
|
||||||
|
|
||||||
if (context == QMCKL_NULL_CONTEXT) then
|
if (context == QMCKL_NULL_CONTEXT) then
|
||||||
@ -3289,6 +3288,7 @@ integer function qmckl_compute_een_rescaled_e_f( &
|
|||||||
een_rescaled_e_ij = 0.0d0
|
een_rescaled_e_ij = 0.0d0
|
||||||
een_rescaled_e_ij(:, 1) = 1.0d0
|
een_rescaled_e_ij(:, 1) = 1.0d0
|
||||||
|
|
||||||
|
|
||||||
k = 0
|
k = 0
|
||||||
do j = 1, elec_num
|
do j = 1, elec_num
|
||||||
do i = 1, j - 1
|
do i = 1, j - 1
|
||||||
@ -3297,6 +3297,7 @@ integer function qmckl_compute_een_rescaled_e_f( &
|
|||||||
end do
|
end do
|
||||||
end do
|
end do
|
||||||
|
|
||||||
|
|
||||||
do l = 2, cord_num
|
do l = 2, cord_num
|
||||||
do k = 1, elec_num * (elec_num - 1)/2
|
do k = 1, elec_num * (elec_num - 1)/2
|
||||||
een_rescaled_e_ij(k, l + 1) = een_rescaled_e_ij(k, l + 1 - 1) * een_rescaled_e_ij(k, 2)
|
een_rescaled_e_ij(k, l + 1) = een_rescaled_e_ij(k, l + 1 - 1) * een_rescaled_e_ij(k, 2)
|
||||||
@ -3305,6 +3306,7 @@ integer function qmckl_compute_een_rescaled_e_f( &
|
|||||||
|
|
||||||
! prepare the actual een table
|
! prepare the actual een table
|
||||||
een_rescaled_e(:, :, 0, nw) = 1.0d0
|
een_rescaled_e(:, :, 0, nw) = 1.0d0
|
||||||
|
|
||||||
do l = 1, cord_num
|
do l = 1, cord_num
|
||||||
k = 0
|
k = 0
|
||||||
do j = 1, elec_num
|
do j = 1, elec_num
|
||||||
@ -3325,28 +3327,14 @@ integer function qmckl_compute_een_rescaled_e_f( &
|
|||||||
|
|
||||||
end do
|
end do
|
||||||
|
|
||||||
end function qmckl_compute_een_rescaled_e_f
|
end function qmckl_compute_een_rescaled_e_doc_f
|
||||||
#+end_src
|
#+end_src
|
||||||
|
|
||||||
#+CALL: generate_c_header(table=qmckl_factor_een_rescaled_e_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
#+CALL: generate_c_interface(table=qmckl_factor_een_rescaled_e_args,rettyp=get_value("CRetType"),fname="qmckl_compute_een_rescaled_e_doc")
|
||||||
|
|
||||||
#+RESULTS:
|
|
||||||
#+begin_src c :tangle (eval h_func) :comments org
|
|
||||||
qmckl_exit_code qmckl_compute_een_rescaled_e (
|
|
||||||
const qmckl_context context,
|
|
||||||
const int64_t walk_num,
|
|
||||||
const int64_t elec_num,
|
|
||||||
const int64_t cord_num,
|
|
||||||
const double rescale_factor_kappa_ee,
|
|
||||||
const double* ee_distance,
|
|
||||||
double* const een_rescaled_e );
|
|
||||||
#+end_src
|
|
||||||
|
|
||||||
#+CALL: generate_c_interface(table=qmckl_factor_een_rescaled_e_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
||||||
|
|
||||||
#+RESULTS:
|
#+RESULTS:
|
||||||
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
||||||
integer(c_int32_t) function qmckl_compute_een_rescaled_e &
|
integer(c_int32_t) function qmckl_compute_een_rescaled_e_doc &
|
||||||
(context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, &
|
(context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, &
|
||||||
ee_distance, een_rescaled_e) &
|
ee_distance, een_rescaled_e) &
|
||||||
bind(C) result(info)
|
bind(C) result(info)
|
||||||
@ -3362,13 +3350,188 @@ end function qmckl_compute_een_rescaled_e_f
|
|||||||
real (c_double ) , intent(in) :: ee_distance(elec_num,elec_num,walk_num)
|
real (c_double ) , intent(in) :: ee_distance(elec_num,elec_num,walk_num)
|
||||||
real (c_double ) , intent(out) :: een_rescaled_e(elec_num,elec_num,0:cord_num,walk_num)
|
real (c_double ) , intent(out) :: een_rescaled_e(elec_num,elec_num,0:cord_num,walk_num)
|
||||||
|
|
||||||
integer(c_int32_t), external :: qmckl_compute_een_rescaled_e_f
|
integer(c_int32_t), external :: qmckl_compute_een_rescaled_e_doc_f
|
||||||
info = qmckl_compute_een_rescaled_e_f &
|
info = qmckl_compute_een_rescaled_e_doc_f &
|
||||||
(context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, ee_distance, een_rescaled_e)
|
(context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, ee_distance, een_rescaled_e)
|
||||||
|
|
||||||
end function qmckl_compute_een_rescaled_e
|
end function qmckl_compute_een_rescaled_e_doc
|
||||||
#+end_src
|
#+end_src
|
||||||
|
|
||||||
|
#+begin_src c :comments org :tangle (eval c) :noweb yes
|
||||||
|
qmckl_exit_code qmckl_compute_een_rescaled_e_hpc (
|
||||||
|
const qmckl_context context,
|
||||||
|
const int64_t walk_num,
|
||||||
|
const int64_t elec_num,
|
||||||
|
const int64_t cord_num,
|
||||||
|
const double rescale_factor_kappa_ee,
|
||||||
|
const double* ee_distance,
|
||||||
|
double* const een_rescaled_e ) {
|
||||||
|
|
||||||
|
double *een_rescaled_e_ij;
|
||||||
|
double x;
|
||||||
|
const int64_t elec_pairs = (elec_num * (elec_num - 1)) / 2;
|
||||||
|
const int64_t len_een_ij = elec_pairs * (cord_num + 1);
|
||||||
|
int64_t k;
|
||||||
|
|
||||||
|
// number of element for the een_rescaled_e_ij[N_e*(N_e-1)/2][cord+1]
|
||||||
|
// probably in C is better [cord+1, Ne*(Ne-1)/2]
|
||||||
|
//elec_pairs = (elec_num * (elec_num - 1)) / 2;
|
||||||
|
//len_een_ij = elec_pairs * (cord_num + 1);
|
||||||
|
een_rescaled_e_ij = (double *) malloc (len_een_ij * sizeof(double));
|
||||||
|
|
||||||
|
if (context == QMCKL_NULL_CONTEXT) {
|
||||||
|
return QMCKL_INVALID_CONTEXT;
|
||||||
|
}
|
||||||
|
|
||||||
|
if (walk_num <= 0) {
|
||||||
|
return QMCKL_INVALID_ARG_2;
|
||||||
|
}
|
||||||
|
|
||||||
|
if (elec_num <= 0) {
|
||||||
|
return QMCKL_INVALID_ARG_3;
|
||||||
|
}
|
||||||
|
|
||||||
|
if (cord_num <= 0) {
|
||||||
|
return QMCKL_INVALID_ARG_4;
|
||||||
|
}
|
||||||
|
|
||||||
|
// Prepare table of exponentiated distances raised to appropriate power
|
||||||
|
// init
|
||||||
|
|
||||||
|
for (int kk = 0; kk < walk_num*(cord_num+1)*elec_num*elec_num; ++kk) {
|
||||||
|
een_rescaled_e[kk]= 0.0;
|
||||||
|
}
|
||||||
|
|
||||||
|
/*
|
||||||
|
for (int nw = 0; nw < walk_num; ++nw) {
|
||||||
|
for (int l = 0; l < (cord_num + 1); ++l) {
|
||||||
|
for (int i = 0; i < elec_num; ++i) {
|
||||||
|
for (int j = 0; j < elec_num; ++j) {
|
||||||
|
een_rescaled_e[j + i*elec_num + l*elec_num*elec_num + nw*(cord_num+1)*elec_num*elec_num]= 0.0;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
*/
|
||||||
|
|
||||||
|
for (int nw = 0; nw < walk_num; ++nw) {
|
||||||
|
|
||||||
|
for (int kk = 0; kk < len_een_ij; ++kk) {
|
||||||
|
// this array initialized at 0 except een_rescaled_e_ij(:, 1) = 1.0d0
|
||||||
|
// and the arrangement of indices is [cord_num+1, ne*(ne-1)/2]
|
||||||
|
een_rescaled_e_ij[kk]= ( kk < (elec_pairs) ? 1.0 : 0.0 );
|
||||||
|
}
|
||||||
|
|
||||||
|
k = 0;
|
||||||
|
for (int i = 0; i < elec_num; ++i) {
|
||||||
|
for (int j = 0; j < i; ++j) {
|
||||||
|
// een_rescaled_e_ij(k, 2) = dexp(-rescale_factor_kappa_ee * ee_distance(i, j, nw));
|
||||||
|
een_rescaled_e_ij[k + elec_pairs] = exp(-rescale_factor_kappa_ee * \
|
||||||
|
ee_distance[j + i*elec_num + nw*(elec_num*elec_num)]);
|
||||||
|
k = k + 1;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
for (int l = 2; l < (cord_num+1); ++l) {
|
||||||
|
for (int k = 0; k < elec_pairs; ++k) {
|
||||||
|
// een_rescaled_e_ij(k, l + 1) = een_rescaled_e_ij(k, l + 1 - 1) * een_rescaled_e_ij(k, 2)
|
||||||
|
een_rescaled_e_ij[k+l*elec_pairs] = een_rescaled_e_ij[k + (l - 1)*elec_pairs] * \
|
||||||
|
een_rescaled_e_ij[k + elec_pairs];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
// prepare the actual een table
|
||||||
|
for (int i = 0; i < elec_num; ++i){
|
||||||
|
for (int j = 0; j < elec_num; ++j) {
|
||||||
|
een_rescaled_e[j + i*elec_num + 0 + nw*(cord_num+1)*elec_num*elec_num] = 1.0;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// Up to here it should work.
|
||||||
|
for ( int l = 1; l < (cord_num+1); ++l) {
|
||||||
|
k = 0;
|
||||||
|
for (int i = 0; i < elec_num; ++i) {
|
||||||
|
for (int j = 0; j < i; ++j) {
|
||||||
|
x = een_rescaled_e_ij[k + l*elec_pairs];
|
||||||
|
een_rescaled_e[j + i*elec_num + l*elec_num*elec_num + nw*elec_num*elec_num*(cord_num+1)] = x;
|
||||||
|
een_rescaled_e[i + j*elec_num + l*elec_num*elec_num + nw*elec_num*elec_num*(cord_num+1)] = x;
|
||||||
|
k = k + 1;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
for (int l = 0; l < (cord_num + 1); ++l) {
|
||||||
|
for (int j = 0; j < elec_num; ++j) {
|
||||||
|
een_rescaled_e[j + j*elec_num + l*elec_num*elec_num + nw*elec_num*elec_num*(cord_num+1)] = 0.0;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
free(een_rescaled_e_ij);
|
||||||
|
|
||||||
|
return QMCKL_SUCCESS;
|
||||||
|
}
|
||||||
|
#+end_src
|
||||||
|
|
||||||
|
#+CALL: generate_c_header(table=qmckl_factor_een_rescaled_e_args,rettyp=get_value("CRetType"),fname="qmckl_compute_een_rescaled_e_doc")
|
||||||
|
|
||||||
|
#+RESULTS:
|
||||||
|
#+begin_src c :tangle (eval h_func) :comments org
|
||||||
|
qmckl_exit_code qmckl_compute_een_rescaled_e (
|
||||||
|
const qmckl_context context,
|
||||||
|
const int64_t walk_num,
|
||||||
|
const int64_t elec_num,
|
||||||
|
const int64_t cord_num,
|
||||||
|
const double rescale_factor_kappa_ee,
|
||||||
|
const double* ee_distance,
|
||||||
|
double* const een_rescaled_e );
|
||||||
|
#+end_src
|
||||||
|
|
||||||
|
#+begin_src c :tangle (eval h_private_func) :comments org
|
||||||
|
qmckl_exit_code qmckl_compute_een_rescaled_e_doc (
|
||||||
|
const qmckl_context context,
|
||||||
|
const int64_t walk_num,
|
||||||
|
const int64_t elec_num,
|
||||||
|
const int64_t cord_num,
|
||||||
|
const double rescale_factor_kappa_ee,
|
||||||
|
const double* ee_distance,
|
||||||
|
double* const een_rescaled_e );
|
||||||
|
#+end_src
|
||||||
|
|
||||||
|
#+begin_src c :tangle (eval h_private_func) :comments org
|
||||||
|
qmckl_exit_code qmckl_compute_een_rescaled_e_hpc (
|
||||||
|
const qmckl_context context,
|
||||||
|
const int64_t walk_num,
|
||||||
|
const int64_t elec_num,
|
||||||
|
const int64_t cord_num,
|
||||||
|
const double rescale_factor_kappa_ee,
|
||||||
|
const double* ee_distance,
|
||||||
|
double* const een_rescaled_e );
|
||||||
|
#+end_src
|
||||||
|
|
||||||
|
#+begin_src c :comments org :tangle (eval c) :noweb yes
|
||||||
|
qmckl_exit_code qmckl_compute_een_rescaled_e (
|
||||||
|
const qmckl_context context,
|
||||||
|
const int64_t walk_num,
|
||||||
|
const int64_t elec_num,
|
||||||
|
const int64_t cord_num,
|
||||||
|
const double rescale_factor_kappa_ee,
|
||||||
|
const double* ee_distance,
|
||||||
|
double* const een_rescaled_e ) {
|
||||||
|
|
||||||
|
#ifdef HAVE_HPC
|
||||||
|
return qmckl_compute_een_rescaled_e_hpc(context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, ee_distance, een_rescaled_e);
|
||||||
|
#else
|
||||||
|
return qmckl_compute_een_rescaled_e_doc(context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, ee_distance, een_rescaled_e);
|
||||||
|
#endif
|
||||||
|
}
|
||||||
|
#+end_src
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
*** Test
|
*** Test
|
||||||
|
|
||||||
#+begin_src python :results output :exports none :noweb yes
|
#+begin_src python :results output :exports none :noweb yes
|
||||||
@ -3443,7 +3606,6 @@ assert(fabs(een_rescaled_e[0][1][0][4]-0.01754273169464735) < 1.e-12);
|
|||||||
assert(fabs(een_rescaled_e[0][2][1][3]-0.02214680362033448) < 1.e-12);
|
assert(fabs(een_rescaled_e[0][2][1][3]-0.02214680362033448) < 1.e-12);
|
||||||
assert(fabs(een_rescaled_e[0][2][1][4]-0.0005700154999202759) < 1.e-12);
|
assert(fabs(een_rescaled_e[0][2][1][4]-0.0005700154999202759) < 1.e-12);
|
||||||
assert(fabs(een_rescaled_e[0][2][1][5]-0.3424402276009091) < 1.e-12);
|
assert(fabs(een_rescaled_e[0][2][1][5]-0.3424402276009091) < 1.e-12);
|
||||||
|
|
||||||
#+end_src
|
#+end_src
|
||||||
|
|
||||||
** Electron-electron rescaled distances for each order and derivatives
|
** Electron-electron rescaled distances for each order and derivatives
|
||||||
@ -5916,6 +6078,7 @@ rc = qmckl_get_jastrow_dtmp_c(context, &(dtmp_c[0][0][0][0][0][0]));
|
|||||||
assert(fabs(tmp_c[0][0][1][0][0] - 2.7083473948352403) < 1e-12);
|
assert(fabs(tmp_c[0][0][1][0][0] - 2.7083473948352403) < 1e-12);
|
||||||
|
|
||||||
assert(fabs(dtmp_c[0][1][0][0][0][0] - 0.237440520852232) < 1e-12);
|
assert(fabs(dtmp_c[0][1][0][0][0][0] - 0.237440520852232) < 1e-12);
|
||||||
|
return QMCKL_SUCCESS;
|
||||||
#+end_src
|
#+end_src
|
||||||
|
|
||||||
** Electron-electron-nucleus Jastrow \(f_{een}\)
|
** Electron-electron-nucleus Jastrow \(f_{een}\)
|
||||||
|
Loading…
Reference in New Issue
Block a user