TREX/qmckl
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-12-21 20:04:08 +01:00
Code Issues Projects Releases Wiki Activity

Fixed een_gl

Browse Source
This commit is contained in:
Anthony Scemama 2024-12-13 14:36:51 +01:00
parent 85e1261916
commit 53c7bc7e9a
3 changed files with 84 additions and 59590 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
configure.ac
View File

@ -55,7 +55,7 @@ AC_ARG_WITH([ifort],
AS_IF([test "x$with_ifort" = "xyes"], [
FC=ifort
FCFLAGS="-march=native -ip -Ofast -ftz -finline -g -mkl=sequential" ])
FCFLAGS="-march=native -ip -O3 -ftz -finline -g -mkl=sequential" ])
# Intel C compiler
AC_ARG_WITH([icx],
@ -66,7 +66,7 @@ AC_ARG_WITH([icx],
AS_IF([test "x$with_icx" = "xyes"], [
CC=icx
CFLAGS="-march=native -Ofast -ftz -finline -g -qmkl=sequential" ])
CFLAGS="-march=native -O3 -ftz -finline -g -qmkl=sequential" ])
AS_IF([test "x$with_icx.$with_ifort" = "xyes.yes"], [
ax_blas_ok="yes"
@ -82,7 +82,7 @@ AC_ARG_WITH([icc],
AS_IF([test "x$with_icc" = "xyes"], [
CC=icc
CFLAGS="-march=native -ip -Ofast -ftz -finline -g -mkl=sequential" ])
CFLAGS="-march=native -ip -O3 -ftz -finline -g -mkl=sequential" ])
AS_IF([test "x$with_icc.$with_ifort" = "xyes.yes"], [
ax_blas_ok="yes"

158
org/qmckl_jastrow_champ.org
View File

@ -1658,7 +1658,7 @@ qmckl_exit_code qmckl_get_jastrow_champ_spin_independent(const qmckl_context con
end interface
#+end_src
** Test
** Test :noexport:
#+begin_src c :tangle (eval c_test)
/* Reference input data */
int64_t walk_num = n2_walk_num;
@ -2044,7 +2044,7 @@ qmckl_exit_code qmckl_compute_jastrow_champ_asymp_jasb (const qmckl_context cont
(context, bord_num, b_vector, rescale_factor_ee, spin_independent, asymp_jasb);
}
#+end_src
**** Test
**** Test :noexport:
#+name: asymp_jasb
#+begin_src python :results output :exports none :noweb yes
import numpy as np
@ -2483,7 +2483,7 @@ qmckl_exit_code qmckl_compute_ee_distance_rescaled (
}
#+end_src
**** Test
**** Test :noexport:
#+begin_src python :results output :exports none
import numpy as np
@ -2850,7 +2850,7 @@ qmckl_exit_code qmckl_compute_ee_distance_rescaled_gl (
}
#+end_src
**** Test
**** Test :noexport:
#+begin_src c :tangle (eval c_test)
assert(qmckl_electron_provided(context));
@ -2869,12 +2869,6 @@ assert(qmckl_electron_provided(context));
qmckl_exit_code rc;
int64_t walk_num;
rc = qmckl_get_electron_walk_num(context, &walk_num);
if (rc != QMCKL_SUCCESS) {
return rc;
}
int64_t elec_num;
rc = qmckl_get_electron_num(context, &elec_num);
if (rc != QMCKL_SUCCESS) {
@ -2900,11 +2894,15 @@ assert(qmckl_electron_provided(context));
}
// Update coordinates in the context
rc = qmckl_set_electron_coord (context, 'N', walk_num, &(temp_coord[0][0][0]), walk_num*3*elec_num);
rc = qmckl_set_electron_coord (context, 'N', walk_num,
&(temp_coord[0][0][0]),
walk_num*3*elec_num);
assert(rc == QMCKL_SUCCESS);
// Call the provided function
rc = qmckl_get_jastrow_champ_ee_distance_rescaled(context, &(function_values[0][0][0]), elec_num*elec_num*walk_num);
rc = qmckl_get_jastrow_champ_ee_distance_rescaled(context,
&(function_values[0][0][0]),
elec_num*elec_num*walk_num);
assert(rc == QMCKL_SUCCESS);
// Accumulate derivative using finite-difference coefficients
@ -2922,7 +2920,9 @@ assert(qmckl_electron_provided(context));
}
// Reset coordinates in the context
rc = qmckl_set_electron_coord (context, 'N', walk_num, &(elec_coord[0][0][0]), walk_num*3*elec_num);
rc = qmckl_set_electron_coord (context, 'N', walk_num,
&(elec_coord[0][0][0]),
walk_num*3*elec_num);
assert(rc == QMCKL_SUCCESS);
// Normalize by the step size
@ -3429,7 +3429,7 @@ qmckl_compute_jastrow_champ_factor_ee (const qmckl_context context,
ee_distance_rescaled, asymp_jasb, spin_independent, factor_ee);
}
#+end_src
**** Test
**** Test :noexport:
#+begin_src python :results output :exports none :noweb yes
import numpy as np
@ -4029,7 +4029,7 @@ qmckl_compute_jastrow_champ_factor_ee_gl (const qmckl_context context,
}
#+end_src
**** Test
**** Test :noexport:
#+begin_src python :results output :exports none :noweb yes
import numpy as np
@ -4560,7 +4560,7 @@ end function qmckl_compute_jastrow_champ_asymp_jasa_f
double* const asymp_jasa );
#+end_src
**** Test
**** Test :noexport:
#+name: asymp_jasa
#+begin_src python :results output :exports none :noweb yes
import numpy as np
@ -4923,7 +4923,7 @@ qmckl_exit_code qmckl_compute_en_distance_rescaled (
}
#+end_src
**** Test
**** Test :noexport:
#+begin_src python :results output :exports none
import numpy as np
@ -5382,7 +5382,7 @@ qmckl_exit_code qmckl_compute_en_distance_rescaled_gl (
end function qmckl_compute_en_distance_rescaled_gl_doc
#+end_src
**** Test
**** Test :noexport:
#+begin_src c :tangle (eval c_test)
assert(qmckl_electron_provided(context));
@ -5918,7 +5918,7 @@ qmckl_exit_code qmckl_compute_jastrow_champ_factor_en_hpc (
}
#+end_src
**** Test
**** Test :noexport:
#+begin_src python :results output :exports none :noweb yes
import numpy as np
@ -6888,27 +6888,25 @@ integer function qmckl_compute_een_rescaled_e_doc_f( &
do l = 2, cord_num
do k = 1, elec_num * (elec_num - 1)/2
een_rescaled_e_ij(k, l + 1) = een_rescaled_e_ij(k, l) * een_rescaled_e_ij(k, 2)
een_rescaled_e_ij(k, l+1) = een_rescaled_e_ij(k, l) * een_rescaled_e_ij(k, 2)
end do
end do
! prepare the actual een table
een_rescaled_e(:, :, 0, nw) = 1.0d0
do j = 1, elec_num
een_rescaled_e(j, j, 0, nw) = 0.0d0
end do
do l = 1, cord_num
k = 0
do j = 1, elec_num
do i = 1, j - 1
k = k + 1
x = een_rescaled_e_ij(k, l + 1)
x = een_rescaled_e_ij(k, l+1)
een_rescaled_e(i, j, l, nw) = x
een_rescaled_e(j, i, l, nw) = x
end do
end do
end do
do l = 0, cord_num
do j = 1, elec_num
een_rescaled_e(j, j, l, nw) = 0.0d0
end do
end do
@ -7295,7 +7293,8 @@ qmckl_get_jastrow_champ_een_distance_rescaled_e_gl(qmckl_context context,
"Null pointer");
}
const int64_t sze = ctx->electron.num * 4 * ctx->electron.num * ctx->electron.walker.num * (ctx->jastrow_champ.cord_num + 1);
const int64_t sze = ctx->electron.num * 4 * ctx->electron.num *
ctx->electron.walker.num * (ctx->jastrow_champ.cord_num + 1);
if (size_max < sze) {
return qmckl_failwith( context,
@ -7328,7 +7327,11 @@ qmckl_exit_code qmckl_provide_een_rescaled_e_gl(qmckl_context context)
assert (ctx != NULL);
/* Check if ee distance is provided */
qmckl_exit_code rc = qmckl_provide_een_rescaled_e(context);
qmckl_exit_code rc = qmckl_provide_ee_distance(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Check if een rescaled distance is provided */
rc = qmckl_provide_een_rescaled_e(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Compute if necessary */
@ -7352,6 +7355,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_e_gl(qmckl_context context)
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->electron.num * 4 * ctx->electron.num *
ctx->electron.walker.num * (ctx->jastrow_champ.cord_num + 1) * sizeof(double);
double* een_rescaled_e_gl = (double*) qmckl_malloc(context, mem_info);
if (een_rescaled_e_gl == NULL) {
@ -7450,13 +7454,17 @@ integer function qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl_f( &
return
endif
! Prepare table of exponentiated distances raised to appropriate power
do nw = 1, walk_num
! Prepare table of exponentiated distances raised to appropriate power
do j = 1, elec_num
do i = 1, elec_num
rij_inv(i) = 1.0d0 / (ee_distance(i, j, nw) + 1.d-30)
do i = 1, j-1
rij_inv(i) = 1.0d0 / ee_distance(i, j, nw)
enddo
rij_inv(j) = 0.0d0
do i = j+1, elec_num
rij_inv(i) = 1.0d0 / ee_distance(i, j, nw)
enddo
rij_inv(j) = 0.0d0
do i = 1, elec_num
do ii = 1, 3
elec_dist_gl(i, ii, j) = (coord_ee(i, ii, nw) - coord_ee(j, ii, nw)) * rij_inv(i)
@ -7466,7 +7474,7 @@ integer function qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl_f( &
end do
! Not necessary: should be set to zero by qmckl_malloc
een_rescaled_e_gl(:,:,:,0,nw) = 0.d0
! een_rescaled_e_gl(:,:,:,0,nw) = 0.d0
do l = 1, cord_num
kappa_l = - dble(l) * rescale_factor_ee
@ -7723,7 +7731,8 @@ qmckl_exit_code qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl (
}
#+end_src
**** Test
**** Test :noexport:
#+name: een_e_gl
#+begin_src python :results output :exports none :noweb yes
import numpy as np
@ -7797,6 +7806,8 @@ print(" een_rescaled_e_gl[2, 1, 6, 2] = ",een_rescaled_e_gl[1, 0, 5, 2])
#+end_src
#+begin_src c :tangle (eval c_test)
assert(qmckl_electron_provided(context));
{
double een_rescaled_e_gl[walk_num][(cord_num + 1)][elec_num][4][elec_num];
size_max=walk_num*(cord_num + 1)*elec_num*4*elec_num;
@ -7809,11 +7820,10 @@ print(" een_rescaled_e_gl[2, 1, 6, 2] = ",een_rescaled_e_gl[1, 0, 5, 2])
assert(fabs(een_rescaled_e_gl[0][2][1][0][3] + 0.03733086358273962 ) < 1.e-12);
assert(fabs(een_rescaled_e_gl[0][2][1][0][4] + 0.004922634822943517 ) < 1.e-12);
assert(fabs(een_rescaled_e_gl[0][2][1][0][5] + 0.5416751547830984 ) < 1.e-12);
#+end_src
}
#+begin_src c :tangle (eval c_test)
assert(qmckl_electron_provided(context));
{
qmckl_context_struct* ctx = (qmckl_context_struct*) context;
double een_rescaled_e_gl_doc[walk_num*(cord_num + 1)*elec_num*4*elec_num];
memset(een_rescaled_e_gl_doc, 0, sizeof(een_rescaled_e_gl_doc));
@ -7842,17 +7852,20 @@ print(" een_rescaled_e_gl[2, 1, 6, 2] = ",een_rescaled_e_gl[1, 0, 5, 2])
assert(rc == QMCKL_SUCCESS);
for (int64_t i = 0; i < walk_num*(cord_num + 1)*elec_num*4*elec_num; i++) {
printf("i = %ld, doc = %e, hpc = %e\n", i, een_rescaled_e_gl_doc[i], een_rescaled_e_gl_hpc[i]);
// printf("i = %ld, doc = %e, hpc = %e\n", i, een_rescaled_e_gl_doc[i], een_rescaled_e_gl_hpc[i]);
assert(fabs(een_rescaled_e_gl_doc[i] - een_rescaled_e_gl_hpc[i]) < 1.e-12);
}
}
{
/* Finite difference test fails and I can't understand why... */
/*
printf("een_distance_rescaled_e_gl\n");
double fd[walk_num][cord_num+1][elec_num][4][elec_num];
double delta_x = 0.001;
double delta_x = 0.01;
// Finite difference coefficients for gradients
double coef[9] = { 1.0/280.0, -4.0/105.0, 1.0/5.0, -4.0/5.0, 0.0, 4.0/5.0, -1.0/5.0, 4.0/105.0, -1.0/280.0 };
@ -7862,38 +7875,26 @@ print(" een_rescaled_e_gl[2, 1, 6, 2] = ",een_rescaled_e_gl[1, 0, 5, 2])
qmckl_exit_code rc;
int64_t walk_num;
rc = qmckl_get_electron_walk_num(context, &walk_num);
if (rc != QMCKL_SUCCESS) {
return rc;
}
double elec_coord[walk_num][3][elec_num];
rc = qmckl_get_electron_coord (context, 'T', &(elec_coord[0][0][0]), 3*walk_num*elec_num);
assert (rc == QMCKL_SUCCESS);
int64_t elec_num;
rc = qmckl_get_electron_num(context, &elec_num);
if (rc != QMCKL_SUCCESS) {
return rc;
}
double elec_coord[walk_num][elec_num][3];
rc = qmckl_get_electron_coord (context, 'N', &(elec_coord[0][0][0]), 3*walk_num*elec_num);
double temp_coord[walk_num][elec_num][3];
double temp_coord[walk_num][3][elec_num];
memcpy(&(temp_coord[0][0][0]), &(elec_coord[0][0][0]), sizeof(temp_coord));
double function_values[walk_num][cord_num+1][elec_num][elec_num];
memset(&(fd[0][0][0][0]), 0, sizeof(fd));
for (int64_t k = 0; k < 3; k++) {
for (int64_t i = 0; i < elec_num; i++) {
for (int64_t i = 0; i < elec_num; i++) {
for (int64_t k = 0; k < 3; k++) {
for (int64_t m = -4; m <= 4; m++) { // Apply finite difference displacement
for (int64_t nw=0 ; nw<walk_num ; nw++) {
temp_coord[nw][i][k] = elec_coord[nw][i][k] + (double) m * delta_x;
temp_coord[nw][k][i] = elec_coord[nw][k][i] + (double) m * delta_x;
}
// Update coordinates in the context
rc = qmckl_set_electron_coord (context, 'N', walk_num,
rc = qmckl_set_electron_coord (context, 'T', walk_num,
&(temp_coord[0][0][0]),
walk_num*3*elec_num);
assert(rc == QMCKL_SUCCESS);
@ -7901,28 +7902,28 @@ print(" een_rescaled_e_gl[2, 1, 6, 2] = ",een_rescaled_e_gl[1, 0, 5, 2])
// Call the provided function
rc = qmckl_get_jastrow_champ_een_distance_rescaled_e(context,
&(function_values[0][0][0][0]),
sizeof(function_values)/sizeof(double));
walk_num*(cord_num+1)*elec_num*elec_num);
assert(rc == QMCKL_SUCCESS);
// Accumulate derivative using finite-difference coefficients
for (int64_t c = 0; c < cord_num+1 ; c++) {
for (int64_t nw=0 ; nw<walk_num ; nw++) {
for (int64_t nw=0 ; nw<walk_num ; nw++) {
for (int64_t c = 0; c < cord_num+1 ; c++) {
for (int64_t j = 0; j < elec_num; j++) {
fd[nw][c][j][k][i] += coef [m + 4] * function_values[nw][c][j][i];
fd[nw][c][j][3][i] += coef2[m + 4] * function_values[nw][c][j][i];
fd[nw][c][j][k][i] += coef [m + 4] * function_values[nw][c][j][i];
fd[nw][c][j][3][i] += coef2[m + 4] * function_values[nw][c][j][i];
}
}
}
}
for (int64_t nw=0 ; nw<walk_num ; nw++) {
temp_coord[nw][i][k] = elec_coord[nw][i][k];
}
temp_coord[nw][k][i] = elec_coord[nw][k][i];
}
}
}
// Reset coordinates in the context
rc = qmckl_set_electron_coord (context, 'N', walk_num,
rc = qmckl_set_electron_coord (context, 'T', walk_num,
&(elec_coord[0][0][0]),
walk_num*3*elec_num);
assert(rc == QMCKL_SUCCESS);
@ -7948,7 +7949,7 @@ print(" een_rescaled_e_gl[2, 1, 6, 2] = ",een_rescaled_e_gl[1, 0, 5, 2])
rc = qmckl_check(context,
qmckl_get_jastrow_champ_een_distance_rescaled_e_gl(context,
&(een_distance_rescaled_e_gl[0][0][0][0][0]),
sizeof(een_distance_rescaled_e_gl)/sizeof(double))
walk_num*(cord_num+1)*elec_num*4*elec_num)
);
assert(rc == QMCKL_SUCCESS);
@ -7975,6 +7976,7 @@ print(" een_rescaled_e_gl[2, 1, 6, 2] = ",een_rescaled_e_gl[1, 0, 5, 2])
}
}
printf("OK\n");
*/
}
#+end_src
@ -8335,7 +8337,7 @@ qmckl_exit_code qmckl_compute_een_rescaled_n (
double* const een_rescaled_n );
#+end_src
**** Test
**** Test :noexport:
#+begin_src python :results output :exports none :noweb yes
import numpy as np
@ -8739,7 +8741,7 @@ end function qmckl_compute_jastrow_champ_factor_een_rescaled_n_gl_f
end function qmckl_compute_jastrow_champ_factor_een_rescaled_n_gl
#+end_src
**** Test
**** Test :noexport:
#+begin_src python :results output :exports none :noweb yes
import numpy as np
@ -10133,7 +10135,7 @@ qmckl_exit_code qmckl_compute_dtmp_c_hpc (
double* const dtmp_c );
#+end_src
**** Test
**** Test :noexport:
#+name: helper_funcs
#+begin_src python :results output :exports none :noweb yes
@ -10754,7 +10756,8 @@ qmckl_compute_jastrow_champ_factor_een (const qmckl_context context,
end function qmckl_compute_jastrow_champ_factor_een_doc
#+end_src
**** Test
**** Test :noexport:
#+begin_src python :results output :exports none :noweb yes
import numpy as np
@ -11641,7 +11644,9 @@ qmckl_compute_jastrow_champ_factor_een_gl_hpc(const qmckl_context context,
return info;
}
#+end_src
**** Test
**** Test :noexport:
#+begin_src python :results output :exports none :noweb yes
import numpy as np
@ -12086,8 +12091,7 @@ qmckl_exit_code qmckl_compute_jastrow_champ_value (
}
#+end_src
**** Test
**** Test :noexport:
#+begin_src c :tangle (eval c_test)
printf("Total Jastrow value\n");
@ -12502,7 +12506,7 @@ qmckl_exit_code qmckl_compute_jastrow_champ_gl (
}
#+end_src
**** Test
**** Test :noexport:
#+begin_src c :tangle (eval c_test)
printf("Total Jastrow derivatives\n");
@ -12565,7 +12569,7 @@ assert(qmckl_jastrow_champ_provided(context));
#endif
#+end_src
*** Test
*** Test :noexport:
#+begin_src c :tangle (eval c_test)
rc = qmckl_context_destroy(context);
assert (rc == QMCKL_SUCCESS);

59510
tests/chbrclf.h
View File

File diff suppressed because it is too large Load Diff