TREX/qmckl
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-08-24 22:21:46 +02:00
Code Issues Projects Releases Wiki Activity

Fixed bug in 3-body Jastrow. Tested in CHAMP->correct

Browse Source
This commit is contained in:
Anthony Scemama 2023-09-12 16:41:26 +02:00
parent fc7d8b2af5
commit 4133a6ba0e
1 changed files with 177 additions and 61 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

238
org/qmckl_jastrow_champ.org
View File

@ -42,7 +42,7 @@
\sum_{p=2}^{N_{\text{ord}}}
\sum_{k=0}^{p-1}
\sum_{l=0}^{p-k-2\delta_{k,0}}
c_{lkp\alpha} \left[ g_\text{ee}({r}_{ij}) \right]^k
c_{lkp\alpha} \left[ g_\text{e}({r}_{ij}) \right]^k
\left[ \left[ g_\alpha({R}_{i\alpha}) \right]^l + \left[ g_\alpha({R}_{j\alpha}) \right]^l \right]
\left[ g_\alpha({R}_{i\,\alpha}) \, g_\alpha({R}_{j\alpha}) \right]^{(p-k-l)/2}
\]
@ -60,7 +60,7 @@
$J_{\text{ee}}$ and $J_{\text{eN}}$ have an asymptotic value of zero.
The eN and eeN parameters are the same of all identical nuclei.
The types of nuclei use zero-based indexing.
Warning: The types of nuclei use zero-based indexing.
* Headers :noexport:
#+begin_src elisp :noexport :results none
@ -581,12 +581,12 @@ qmckl_set_jastrow_champ_type_nucl_vector(qmckl_context context,
for (int i=0 ; i<nucl_num ; ++i) {
if (type_nucl_vector[i] < 0) {
return qmckl_failwith( context, QMCKL_INVALID_ARG_2,
"qmckl_set_type_nucl_vector",
"qmckl_set_jastrow_champ_type_nucl_vector",
"Inconsistent values of type_nucl_vector (<0)" );
}
if (type_nucl_vector[i] >= type_nucl_num) {
return qmckl_failwith( context, QMCKL_INVALID_ARG_2,
"qmckl_set_type_nucl_vector",
"qmckl_set_jastrow_champ_type_nucl_vector",
"Inconsistent values of type_nucl_vector (>=nucl_num). Values should use 0-based indexing as in C." );
}
}
@ -595,7 +595,7 @@ qmckl_set_jastrow_champ_type_nucl_vector(qmckl_context context,
qmckl_exit_code rc = qmckl_free(context, ctx->jastrow_champ.type_nucl_vector);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context, rc,
"qmckl_set_type_nucl_vector",
"qmckl_set_jastrow_champ_type_nucl_vector",
"Unable to free ctx->jastrow_champ.type_nucl_vector");
}
}
@ -676,7 +676,7 @@ qmckl_set_jastrow_champ_a_vector(qmckl_context context,
if(new_array == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_set_jastrow_champ_coefficient",
"qmckl_set_jastrow_champ_a_vector",
NULL);
}
@ -736,7 +736,7 @@ qmckl_set_jastrow_champ_b_vector(qmckl_context context,
if(new_array == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_set_jastrow_champ_coefficient",
"qmckl_set_jastrow_champ_b_vector",
NULL);
}
@ -807,7 +807,7 @@ qmckl_set_jastrow_champ_c_vector(qmckl_context context,
if(new_array == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_set_jastrow_champ_coefficient",
"qmckl_set_jastrow_champ_c_vector",
NULL);
}
@ -2008,7 +2008,7 @@ rc = qmckl_check(context,
);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_check(context,
qmckl_set_jastrow_champ_c_vector(context, c_vector,dim_c_vector*type_nucl_num)
qmckl_set_jastrow_champ_c_vector(context, c_vector, dim_c_vector*type_nucl_num)
);
assert(rc == QMCKL_SUCCESS);
@ -7589,15 +7589,15 @@ qmckl_exit_code qmckl_provide_dtmp_c(qmckl_context context)
:END:
#+NAME: qmckl_factor_c_vector_full_args
| Variable | Type | In/Out | Description |
|--------------------+----------------------------------------+--------+------------------------------|
| ~context~ | ~qmckl_context~ | in | Global state |
| ~nucl_num~ | ~int64_t~ | in | Number of atoms |
| ~dim_c_vector~ | ~int64_t~ | in | dimension of cord full table |
| ~type_nucl_num~ | ~int64_t~ | in | dimension of cord full table |
| ~type_nucl_vector~ | ~int64_t[nucl_num]~ | in | dimension of cord full table |
| Variable | Type | In/Out | Description |
|--------------------+---------------------------------------+--------+------------------------------|
| ~context~ | ~qmckl_context~ | in | Global state |
| ~nucl_num~ | ~int64_t~ | in | Number of atoms |
| ~dim_c_vector~ | ~int64_t~ | in | dimension of cord full table |
| ~type_nucl_num~ | ~int64_t~ | in | dimension of cord full table |
| ~type_nucl_vector~ | ~int64_t[nucl_num]~ | in | dimension of cord full table |
| ~c_vector~ | ~double[dim_c_vector][type_nucl_num]~ | in | dimension of cord full table |
| ~c_vector_full~ | ~double[dim_c_vector][nucl_num]~ | out | Full list of coefficients |
| ~c_vector_full~ | ~double[nucl_num][dim_c_vector]~ | out | Full list of coefficients |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_c_vector_full_doc_f( &
@ -7611,8 +7611,8 @@ integer function qmckl_compute_c_vector_full_doc_f( &
integer*8 , intent(in) :: dim_c_vector
integer*8 , intent(in) :: type_nucl_num
integer*8 , intent(in) :: type_nucl_vector(nucl_num)
double precision , intent(in) :: c_vector(type_nucl_num, dim_c_vector)
double precision , intent(out) :: c_vector_full(nucl_num,dim_c_vector)
double precision , intent(in) :: c_vector(dim_c_vector, type_nucl_num)
double precision , intent(out) :: c_vector_full(nucl_num, dim_c_vector)
double precision :: x
integer*8 :: i, a, k, l, nw
@ -7625,7 +7625,7 @@ integer function qmckl_compute_c_vector_full_doc_f( &
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)
c_vector_full(a,1:dim_c_vector) = c_vector(1:dim_c_vector, type_nucl_vector(a)+1)
end do
end function qmckl_compute_c_vector_full_doc_f
@ -7677,7 +7677,7 @@ qmckl_exit_code qmckl_compute_c_vector_full_hpc (
for (int i=0; i < dim_c_vector; ++i) {
for (int a=0; a < nucl_num; ++a){
c_vector_full[a + i*nucl_num] = c_vector[(type_nucl_vector[a])+i*type_nucl_num];
c_vector_full[a + i*nucl_num] = c_vector[i + type_nucl_vector[a]*dim_c_vector];
}
}
@ -7746,7 +7746,7 @@ qmckl_exit_code qmckl_compute_c_vector_full (
:FRetType: qmckl_exit_code
:END:
#+NAME: qmckl_factor_lkpm_combined_index_args
#+NAME: lkpm_combined_index_args
| Variable | Type | In/Out | Description |
|-----------------------+----------------------------+--------+-------------------------------|
| ~context~ | ~qmckl_context~ | in | Global state |
@ -7755,7 +7755,7 @@ qmckl_exit_code qmckl_compute_c_vector_full (
| ~lkpm_combined_index~ | ~int64_t[4][dim_c_vector]~ | out | Full list of combined indices |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_lkpm_combined_index_f( &
integer function qmckl_compute_lkpm_combined_index_doc_f( &
context, cord_num, dim_c_vector, lkpm_combined_index) &
result(info)
use qmckl
@ -7777,13 +7777,13 @@ integer function qmckl_compute_lkpm_combined_index_f( &
kk = 0
do p = 2, cord_num
do k = p - 1, 0, -1
if (k .ne. 0) then
if (k /= 0) then
lmax = p - k
else
lmax = p - k - 2
end if
do l = lmax, 0, -1
if (iand(p - k - l, 1_8) .eq. 1) cycle
if (iand(p - k - l, 1_8) .eq. 1_8) cycle
m = (p - k - l)/2
kk = kk + 1
lkpm_combined_index(kk, 1) = l
@ -7794,11 +7794,11 @@ integer function qmckl_compute_lkpm_combined_index_f( &
end do
end do
end function qmckl_compute_lkpm_combined_index_f
end function qmckl_compute_lkpm_combined_index_doc_f
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes
qmckl_exit_code qmckl_compute_lkpm_combined_index (
qmckl_exit_code qmckl_compute_lkpm_combined_index_hpc (
const qmckl_context context,
const int64_t cord_num,
const int64_t dim_c_vector,
@ -7834,8 +7834,85 @@ qmckl_exit_code qmckl_compute_lkpm_combined_index (
}
#+end_src
# #+CALL: generate_c_header(table=qmckl_factor_lkpm_combined_index_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
#+CALL: generate_c_interface(table=lkpm_combined_index_args,rettyp=get_value("CRetType"),fname="qmckl_compute_lkpm_combined_index_doc")
#+RESULTS:
#+begin_src f90 :tangle (eval f) :comments org :exports none
integer(c_int32_t) function qmckl_compute_lkpm_combined_index_doc &
(context, cord_num, dim_c_vector, lkpm_combined_index) &
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 :: cord_num
integer (c_int64_t) , intent(in) , value :: dim_c_vector
integer (c_int64_t) , intent(out) :: lkpm_combined_index(dim_c_vector,4)
integer(c_int32_t), external :: qmckl_compute_lkpm_combined_index_doc_f
info = qmckl_compute_lkpm_combined_index_doc_f &
(context, cord_num, dim_c_vector, lkpm_combined_index)
end function qmckl_compute_lkpm_combined_index_doc
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes
qmckl_exit_code qmckl_compute_lkpm_combined_index (
const qmckl_context context,
const int64_t cord_num,
const int64_t dim_c_vector,
int64_t* const lkpm_combined_index ) {
#ifdef HAVE_HPC
return qmckl_compute_lkpm_combined_index_hpc(context, cord_num, dim_c_vector, lkpm_combined_index);
#else
return qmckl_compute_lkpm_combined_index_doc(context, cord_num, dim_c_vector, lkpm_combined_index);
#endif
}
#+end_src
#+CALL: generate_c_header(table=lkpm_combined_index_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_lkpm_combined_index (
const qmckl_context context,
const int64_t cord_num,
const int64_t dim_c_vector,
int64_t* const lkpm_combined_index );
#+end_src
#+CALL: generate_c_header(table=lkpm_combined_index_args,rettyp=get_value("CRetType"),fname="qmckl_compute_lkpm_combined_index_doc")
#+RESULTS:
#+begin_src c :tangle (eval h_func) :comments org
qmckl_exit_code qmckl_compute_lkpm_combined_index_doc (
const qmckl_context context,
const int64_t cord_num,
const int64_t dim_c_vector,
int64_t* const lkpm_combined_index );
#+end_src
#+CALL: generate_c_header(table=lkpm_combined_index_args,rettyp=get_value("CRetType"),fname="qmckl_compute_lkpm_combined_index_hpc")
#+RESULTS:
#+begin_src c :tangle (eval h_func) :comments org
qmckl_exit_code qmckl_compute_lkpm_combined_index_hpc (
const qmckl_context context,
const int64_t cord_num,
const int64_t dim_c_vector,
int64_t* const lkpm_combined_index );
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_compute_lkpm_combined_index (
const qmckl_context context,
const int64_t cord_num,
const int64_t dim_c_vector,
int64_t* const lkpm_combined_index );
#+end_src
#+CALL: generate_c_header(table=lkpm_combined_index_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_lkpm_combined_index (
const qmckl_context context,
@ -8526,6 +8603,20 @@ qmckl_exit_code qmckl_provide_jastrow_champ_factor_een(qmckl_context context)
ctx->jastrow_champ.factor_een = factor_een;
}
/*
rc = qmckl_compute_jastrow_champ_factor_een_naive(context,
ctx->electron.walker.num,
ctx->electron.num,
ctx->nucleus.num,
ctx->jastrow_champ.cord_num,
ctx->jastrow_champ.dim_c_vector,
ctx->jastrow_champ.c_vector_full,
ctx->jastrow_champ.lkpm_combined_index,
ctx->jastrow_champ.een_rescaled_e,
ctx->jastrow_champ.een_rescaled_n,
ctx->jastrow_champ.factor_een);
*/
rc = qmckl_compute_jastrow_champ_factor_een(context,
ctx->electron.walker.num,
ctx->electron.num,
@ -8566,8 +8657,8 @@ qmckl_exit_code qmckl_provide_jastrow_champ_factor_een(qmckl_context context)
| ~dim_c_vector~ | ~int64_t~ | in | dimension of full coefficient vector |
| ~c_vector_full~ | ~double[dim_c_vector][nucl_num]~ | in | full coefficient vector |
| ~lkpm_combined_index~ | ~int64_t[4][dim_c_vector]~ | in | combined indices |
| ~een_rescaled_e~ | ~double[walk_num][elec_num][elec_num][0:cord_num]~ | in | Electron-nucleus rescaled |
| ~een_rescaled_n~ | ~double[walk_num][elec_num][nucl_num][0:cord_num]~ | in | Electron-nucleus rescaled factor |
| ~een_rescaled_e~ | ~double[walk_num][0:cord_num][elec_num][elec_num]~ | in | Electron-nucleus rescaled |
| ~een_rescaled_n~ | ~double[walk_num][0:cord_num][nucl_num][elec_num]~ | in | Electron-nucleus rescaled factor |
| ~factor_een~ | ~double[walk_num]~ | out | Electron-nucleus jastrow |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
@ -8582,8 +8673,8 @@ integer function qmckl_compute_jastrow_champ_factor_een_naive_f( &
integer*8 , intent(in) :: walk_num, elec_num, cord_num, nucl_num, dim_c_vector
integer*8 , intent(in) :: lkpm_combined_index(dim_c_vector,4)
double precision , intent(in) :: c_vector_full(nucl_num, dim_c_vector)
double precision , intent(in) :: een_rescaled_e(0:cord_num, elec_num, elec_num, walk_num)
double precision , intent(in) :: een_rescaled_n(0:cord_num, nucl_num, elec_num, walk_num)
double precision , intent(in) :: een_rescaled_e(elec_num, elec_num, 0:cord_num, walk_num)
double precision , intent(in) :: een_rescaled_n(elec_num, nucl_num, 0:cord_num, walk_num)
double precision , intent(out) :: factor_een(walk_num)
integer*8 :: i, a, j, l, k, p, m, n, nw
@ -8597,34 +8688,58 @@ integer function qmckl_compute_jastrow_champ_factor_een_naive_f( &
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
! do nw =1, walk_num
! factor_een(nw) = 0.0d0
! do n = 1, dim_c_vector
! l = lkpm_combined_index(n, 1)
! k = lkpm_combined_index(n, 2)
! p = lkpm_combined_index(n, 3)
! m = lkpm_combined_index(n, 4)
! do a = 1, nucl_num
! accu2 = 0.0d0
! cn = c_vector_full(a, n)
! do j = 1, elec_num
! accu = 0.0d0
! do i = 1, elec_num
! accu = accu + een_rescaled_e(i,j,k,nw) * &
! een_rescaled_n(i,a,m,nw)
! end do
! accu2 = accu2 + accu * een_rescaled_n(j,a,m + l,nw)
! end do
! factor_een(nw) = factor_een(nw) + accu2 * cn
! end do
! end do
! end do
do nw =1, walk_num
do n = 1, dim_c_vector
l = lkpm_combined_index(n, 1)
k = lkpm_combined_index(n, 2)
p = lkpm_combined_index(n, 3)
m = lkpm_combined_index(n, 4)
factor_een(nw) = 0.d0
do n = 1, dim_c_vector
l = lkpm_combined_index(n, 1)
k = lkpm_combined_index(n, 2)
p = lkpm_combined_index(n, 3)
m = lkpm_combined_index(n, 4)
do a = 1, nucl_num
accu2 = 0.0d0
cn = c_vector_full(a, n)
do j = 1, elec_num
accu = 0.0d0
do i = 1, elec_num
accu = accu + een_rescaled_e(k,i,j,nw) * &
een_rescaled_n(m,a,i,nw)
!if(nw .eq. 1) then
! print *,l,k,p,m,j,i,een_rescaled_e(k,i,j,nw), een_rescaled_n(m,a,i,nw), accu
!endif
do a = 1, nucl_num
accu2 = 0.0d0
cn = c_vector_full(a, n)
print *, a, l, k, p, cn
do j = 1, elec_num
accu = 0.0d0
do i = 1, j-1
accu = accu + een_rescaled_e(i,j,k,nw) * &
(een_rescaled_n(i,a,l,nw) + een_rescaled_n(j,a,l,nw)) * &
(een_rescaled_n(i,a,m,nw) * een_rescaled_n(j,a,m,nw))
end do
accu2 = accu2 + accu
end do
factor_een(nw) = factor_een(nw) + accu2 * cn
end do
accu2 = accu2 + accu * een_rescaled_n(m + l,a,j,nw)
!print *, l,m,nw,accu, accu2, een_rescaled_n(m + l, a, j, nw), cn, factor_een(nw)
end do
factor_een(nw) = factor_een(nw) + accu2 * cn
end do
end do
end do
end do
end function qmckl_compute_jastrow_champ_factor_een_naive_f
@ -8786,11 +8901,12 @@ qmckl_compute_jastrow_champ_factor_een_doc (const qmckl_context context,
const int64_t dim_c_vector,
const double* c_vector_full,
const int64_t* lkpm_combined_index,
const double* een_rescaled_e,
const double* tmp_c,
const double* een_rescaled_n,
double* const factor_een );
qmckl_exit_code
qmckl_compute_jastrow_champ_factor_een (const qmckl_context context,
const int64_t walk_num,
const int64_t elec_num,
@ -8799,7 +8915,7 @@ qmckl_compute_jastrow_champ_factor_een (const qmckl_context context,
const int64_t dim_c_vector,
const double* c_vector_full,
const int64_t* lkpm_combined_index,
const double* een_rescaled_e,
const double* tmp_c,
const double* een_rescaled_n,
double* const factor_een );
#+end_src
@ -8814,17 +8930,17 @@ qmckl_compute_jastrow_champ_factor_een (const qmckl_context context,
const int64_t dim_c_vector,
const double* c_vector_full,
const int64_t* lkpm_combined_index,
const double* een_rescaled_e,
const double* tmp_c,
const double* een_rescaled_n,
double* const factor_een )
{
#ifdef HAVE_HPC
return qmckl_compute_jastrow_champ_factor_een_doc (context, walk_num, elec_num, nucl_num, cord_num, dim_c_vector,
c_vector_full, lkpm_combined_index, een_rescaled_e, een_rescaled_n,
c_vector_full, lkpm_combined_index, tmp_c, een_rescaled_n,
factor_een );
#else
return qmckl_compute_jastrow_champ_factor_een_doc (context, walk_num, elec_num, nucl_num, cord_num, dim_c_vector,
c_vector_full, lkpm_combined_index, een_rescaled_e, een_rescaled_n,
c_vector_full, lkpm_combined_index, tmp_c, een_rescaled_n,
factor_een );
#endif
}