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

Fixed bug in HPC version of AOs

Browse Source
This commit is contained in:
Anthony Scemama 2022-06-15 23:21:31 +02:00
parent 28fe475425
commit 1b6cf47f0d
1 changed files with 66 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

80
org/qmckl_ao.org
View File

@ -3528,7 +3528,7 @@ integer function qmckl_compute_ao_basis_primitive_gaussian_vgl_f( &
info = QMCKL_SUCCESS
! Don't compute exponentials when the result will be almost zero.
cutoff = -dlog(1.d-15)
cutoff = 27.631021115928547 ! -dlog(1.d-12)
do inucl=1,nucl_num
! C is zero-based, so shift bounds by one
@ -3752,10 +3752,15 @@ print ( "[7][4][26] : %e"% lf(a,x,y))
(int64_t) 5*elec_num*walk_num*prim_num );
assert (rc == QMCKL_SUCCESS);
printf("prim_vgl[26][0][7] = %e\n",prim_vgl[26][0][7]);
assert( fabs(prim_vgl[26][0][7] - ( 1.0501570432064878E-003)) < 1.e-14 );
printf("prim_vgl[26][1][7] = %e\n",prim_vgl[26][1][7]);
assert( fabs(prim_vgl[26][1][7] - (-7.5014974095310560E-004)) < 1.e-14 );
printf("prim_vgl[26][2][7] = %e\n",prim_vgl[26][2][7]);
assert( fabs(prim_vgl[26][2][7] - (-3.8250692897610380E-003)) < 1.e-14 );
printf("prim_vgl[26][3][7] = %e\n",prim_vgl[26][3][7]);
assert( fabs(prim_vgl[26][3][7] - ( 3.4950559194080275E-003)) < 1.e-14 );
printf("prim_vgl[26][4][7] = %e\n",prim_vgl[26][4][7]);
assert( fabs(prim_vgl[26][4][7] - ( 2.0392163767356572E-002)) < 1.e-14 );
}
@ -3875,7 +3880,7 @@ integer function qmckl_compute_ao_basis_shell_gaussian_vgl_f( &
! Don't compute exponentials when the result will be almost zero.
! TODO : Use numerical precision here
cutoff = -dlog(1.d-12)
cutoff = 27.631021115928547 !-dlog(1.d-12)
do ipoint = 1, point_num
@ -5501,7 +5506,7 @@ integer function qmckl_compute_ao_value_doc_f(context, &
info = QMCKL_SUCCESS
! Don't compute polynomials when the radial part is zero.
cutoff = -dlog(1.d-12)
cutoff = 27.631021115928547 !-dlog(1.d-12)
do ipoint = 1, point_num
e_coord(1) = coord(ipoint,1)
@ -5622,7 +5627,7 @@ qmckl_compute_ao_value_hpc_gaussian (const qmckl_context context,
}
/* Don't compute polynomials when the radial part is zero. */
double cutoff = -log(1.e-12);
double cutoff = 27.631021115928547; // -log(1.e-12)
#ifdef HAVE_OPENMP
#pragma omp parallel
@ -6247,7 +6252,7 @@ integer function qmckl_compute_ao_vgl_doc_f(context, &
info = QMCKL_SUCCESS
! Don't compute polynomials when the radial part is zero.
cutoff = -dlog(1.d-12)
cutoff = 27.631021115928547 ! -dlog(1.d-12)
do ipoint = 1, point_num
e_coord(1) = coord(ipoint,1)
@ -6398,8 +6403,8 @@ qmckl_compute_ao_vgl_hpc_gaussian (
ao_index[shell_num] = ao_num+1;
}
/* Don't compute polynomials when the radial part is zero. */
double cutoff = -log(1.e-12);
/* Don't compute when the radial part is zero. */
double cutoff = 27.631021115928547; // -log(1.e-12)
#ifdef HAVE_OPENMP
#pragma omp parallel
@ -6637,11 +6642,6 @@ qmckl_compute_ao_vgl_hpc_gaussian (
for (int64_t ishell = ishell_start ; ishell < ishell_end ; ++ishell) {
const double s1 = ce_mat[ishell-ishell_start][0];
if (s1 == 0.0) continue;
const double s2 = ce_mat[ishell-ishell_start][1];
const double s3 = ce_mat[ishell-ishell_start][2];
const double s4 = ce_mat[ishell-ishell_start][3];
const double s5 = ce_mat[ishell-ishell_start][4];
const int64_t k = ao_index[ishell];
double* restrict const ao_vgl_1 = ao_vgl + ipoint*5*ao_num + k;
@ -6654,6 +6654,22 @@ qmckl_compute_ao_vgl_hpc_gaussian (
double* restrict const ao_vgl_4 = ao_vgl_1 + (ao_num<<1) + ao_num;
double* restrict const ao_vgl_5 = ao_vgl_1 + (ao_num<<2);
if (s1 == 0.0) {
for (int64_t il=0 ; il<n ; ++il) {
ao_vgl_1[il] = 0.0;
ao_vgl_2[il] = 0.0;
ao_vgl_3[il] = 0.0;
ao_vgl_4[il] = 0.0;
ao_vgl_5[il] = 0.0;
}
continue;
}
const double s2 = ce_mat[ishell-ishell_start][1];
const double s3 = ce_mat[ishell-ishell_start][2];
const double s4 = ce_mat[ishell-ishell_start][3];
const double s5 = ce_mat[ishell-ishell_start][4];
double* restrict poly_vgl_1 = NULL;
double* restrict poly_vgl_2 = NULL;
double* restrict poly_vgl_3 = NULL;
@ -6726,7 +6742,7 @@ qmckl_compute_ao_vgl_hpc_gaussian (
break;
default:
#ifdef HAVE_OPENMP
#pragma omp simd simdlen(8)
#pragma omp simd
#endif
for (int il=0 ; il<n ; ++il) {
ao_vgl_1[il] = poly_vgl_1[il] * s1 * f[il];
@ -6900,7 +6916,9 @@ qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context)
qmckl_exit_code rc;
/* Provide required data */
#ifndef HAVE_HPC
#ifdef HAVE_HPC
#else
rc = qmckl_provide_ao_basis_shell_vgl(context);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context, rc, "qmckl_provide_ao_basis_shell_vgl", NULL);
@ -6941,6 +6959,40 @@ qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context)
ctx->ao_basis.expo_per_nucleus,
ctx->ao_basis.coef_per_nucleus,
ctx->ao_basis.ao_vgl);
/* DEBUG
rc = qmckl_provide_ao_basis_shell_vgl(context);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context, rc, "qmckl_provide_ao_basis_shell_vgl", NULL);
}
int64_t K= ctx->ao_basis.ao_num * 5 * ctx->point.num;
double* check = malloc(K*sizeof(double));
rc = qmckl_compute_ao_vgl_doc(context,
ctx->ao_basis.ao_num,
ctx->ao_basis.shell_num,
ctx->point.num,
ctx->nucleus.num,
ctx->point.coord.data,
ctx->nucleus.coord.data,
ctx->ao_basis.nucleus_index,
ctx->ao_basis.nucleus_shell_num,
ctx->ao_basis.nucleus_range,
ctx->ao_basis.nucleus_max_ang_mom,
ctx->ao_basis.shell_ang_mom,
ctx->ao_basis.ao_factor,
ctx->ao_basis.shell_vgl,
check);
for (int64_t i=0 ; i<K ; ++i) {
if (fabs(check[i] - ctx->ao_basis.ao_vgl[i]) > 1.e-10) {
int a, b, c;
a = i/(ctx->ao_basis.ao_num*5);
b = (i-a*ctx->ao_basis.ao_num*5)/ctx->ao_basis.ao_num;
c = (i-a*ctx->ao_basis.ao_num*5 -b*ctx->ao_basis.ao_num);
printf("%d: %d, %d, %d, %e %e\n", i, a, b, c, check[i], ctx->ao_basis.ao_vgl[i]);
}
}
*/
/*
} else if (ctx->ao_basis.type == 'S') {
rc = qmck_compute_ao_vgl_hpc_slater(context,