From 1b6cf47f0dd81c155bff1b2e3fa04cb5b4f0f64e Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 15 Jun 2022 23:21:31 +0200 Subject: [PATCH] Fixed bug in HPC version of AOs --- org/qmckl_ao.org | 80 +++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 66 insertions(+), 14 deletions(-) diff --git a/org/qmckl_ao.org b/org/qmckl_ao.org index 43ca5bc..b53e4f7 100644 --- a/org/qmckl_ao.org +++ b/org/qmckl_ao.org @@ -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 ; ilao_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 ; iao_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,