From 9553059bbe831edc4b36b4d44ef8f1fcdf482239 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 8 May 2024 13:56:32 +0200 Subject: [PATCH] Fix wrong gradient at nodes of AOs --- org/qmckl_ao.org | 88 ++++++++++++++++++++++++++---------------------- org/qmckl_mo.org | 3 +- 2 files changed, 50 insertions(+), 41 deletions(-) diff --git a/org/qmckl_ao.org b/org/qmckl_ao.org index a50e560..29322ae 100644 --- a/org/qmckl_ao.org +++ b/org/qmckl_ao.org @@ -3943,6 +3943,7 @@ function qmckl_compute_ao_basis_shell_gaussian_vgl( & do inucl=1,nucl_num + ! The shift below avoids having an exact zero on a node of the orbital x = coord(ipoint,1) - nucl_coord(inucl,1) y = coord(ipoint,2) - nucl_coord(inucl,2) z = coord(ipoint,3) - nucl_coord(inucl,3) @@ -4569,8 +4570,9 @@ function qmckl_ao_polynomial_vgl_doc (context, & endif + ! The shift below is such that polynomials will not make the AO equal to zero at the nodes of the orbitals do i=1,3 - Y(i) = X(i) - R(i) + Y(i) = X(i) - R(i) + 1.d-20 end do if (lmax == 0) then @@ -4927,6 +4929,7 @@ qmckl_ao_polynomial_transp_vgl_hpc_inline (const qmckl_context context, if (ldl < 3) return QMCKL_INVALID_ARG_7; int32_t m; + const double shift=1.e-20; switch (lmax) { case 0: @@ -4964,7 +4967,7 @@ qmckl_ao_polynomial_transp_vgl_hpc_inline (const qmckl_context context, if (ldv == 4) { - const double tmp[20] = {1.0, X[0]-R[0], X[1]-R[1], X[2]-R[2], + const double tmp[20] = {1.0, X[0]-R[0]+shift, X[1]-R[1]+shift, X[2]-R[2]+shift, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, @@ -4988,9 +4991,9 @@ qmckl_ao_polynomial_transp_vgl_hpc_inline (const qmckl_context context, vgl5[k] = 0.0; } vgl1[0] = 1.0; - vgl1[1] = X[0]-R[0]; - vgl1[2] = X[1]-R[1]; - vgl1[3] = X[2]-R[2]; + vgl1[1] = X[0]-R[0]+shift; + vgl1[2] = X[1]-R[1]+shift; + vgl1[3] = X[2]-R[2]+shift; vgl2[1] = 1.0; vgl3[2] = 1.0; vgl4[3] = 1.0; @@ -5032,10 +5035,11 @@ qmckl_ao_polynomial_transp_vgl_hpc_inline (const qmckl_context context, if (ldv == 50) { const double tmp[50] = { - 1.0, Y[0], Y[1], Y[2], Y[0] * Y[0], - Y[0] * Y[1], Y[0] * Y[2], Y[1] * Y[1], Y[1] * Y[2], Y[2] * Y[2], - 0.0, 1.0, 0.0, 0.0, Y[0] + Y[0], - Y[1], Y[2], 0.0, 0.0, 0.0, + 1.0, Y[0]+shift, Y[1]+shift, Y[2]+shift, + Y[0] * Y[0]+shift, Y[0] * Y[1]+shift, Y[0] * Y[2]+shift, + Y[1] * Y[1]+shift, Y[1] * Y[2]+shift, Y[2] * Y[2]+shift, + 0.0, 1.0, 0.0, 0.0, + Y[0] + Y[0], Y[1], Y[2], 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, Y[0], 0.0, Y[1] + Y[1], Y[2], 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, @@ -5051,10 +5055,10 @@ qmckl_ao_polynomial_transp_vgl_hpc_inline (const qmckl_context context, double* restrict const vgl4 = VGL + 3*ldv; double* restrict const vgl5 = VGL + (ldv << 2); - vgl1[0] = 1.0 ; vgl1[1] = Y[0] ; vgl1[2] = Y[1]; - vgl1[3] = Y[2] ; vgl1[4] = Y[0]*Y[0]; vgl1[5] = Y[0]*Y[1]; - vgl1[6] = Y[0]*Y[2]; vgl1[7] = Y[1]*Y[1]; vgl1[8] = Y[1]*Y[2]; - vgl1[9] = Y[2]*Y[2]; + vgl1[0] = 1.0 ; vgl1[1] = Y[0]+shift ; vgl1[2] = Y[1]+shift; + vgl1[3] = Y[2]+shift ; vgl1[4] = Y[0]*Y[0]+shift; vgl1[5] = Y[0]*Y[1]+shift; + vgl1[6] = Y[0]*Y[2]+shift; vgl1[7] = Y[1]*Y[1]+shift; vgl1[8] = Y[1]*Y[2]+shift; + vgl1[9] = Y[2]*Y[2]+shift; vgl2[0] = 0.0 ; vgl2[1] = 1.0 ; vgl2[2] = 0.0 ; vgl2[3] = 0.0 ; vgl2[4] = Y[0]+Y[0]; vgl2[5] = Y[1]; @@ -5092,7 +5096,7 @@ qmckl_ao_polynomial_transp_vgl_hpc_inline (const qmckl_context context, const double Y[3] = { X[0]-R[0], X[1]-R[1], - X[2]-R[2] }; + X[2]-R[2]}; assert(size_max > lmax+3); double pows[3][size_max]; @@ -5130,9 +5134,9 @@ qmckl_ao_polynomial_transp_vgl_hpc_inline (const qmckl_context context, vgl5[k] = 0.0; } vgl1[0] = 1.0; - vgl1[1] = Y[0]; - vgl1[2] = Y[1]; - vgl1[3] = Y[2]; + vgl1[1] = Y[0]+shift; + vgl1[2] = Y[1]+shift; + vgl1[3] = Y[2]+shift; vgl2[1] = 1.0; vgl3[2] = 1.0; @@ -5159,7 +5163,7 @@ qmckl_ao_polynomial_transp_vgl_hpc_inline (const qmckl_context context, l[m][1] = b; l[m][2] = c; - vgl1[m] = xy * pows[2][c+2]; + vgl1[m] = xy * pows[2][c+2]+shift; xy *= dc; xz *= db; @@ -5366,7 +5370,6 @@ for (int32_t ldl=3 ; ldl<=5 ; ++ldl) { * Combining radial and polynomial parts - ** Determination of nucleus ranges The range of a nucleus is defined as the distance beyond which all @@ -5793,7 +5796,7 @@ qmckl_compute_ao_value_hpc_gaussian (const qmckl_context context, for (size_t k=0 ; k cutoff * nucleus_range[inucl]) { continue; } @@ -5836,9 +5842,9 @@ qmckl_compute_ao_value_hpc_gaussian (const qmckl_context context, case 2: poly_vgl[0] = 1.; - poly_vgl[1] = x; - poly_vgl[2] = y; - poly_vgl[3] = z; + poly_vgl[1] = x+shift; + poly_vgl[2] = y+shift; + poly_vgl[3] = z+shift; poly_vgl[4] = x2; poly_vgl[5] = x*y; poly_vgl[6] = x*z; @@ -6584,7 +6590,7 @@ qmckl_compute_ao_vgl_hpc_gaussian ( for (size_t k=0 ; k<5*ao_num ; ++k) { ptr[k] = 0.; } -*/ +,*/ memset(&ao_vgl[ipoint*ao_num*5], 0, 5*ao_num*sizeof(double)); @@ -6599,7 +6605,8 @@ qmckl_compute_ao_vgl_hpc_gaussian ( nucl_coord[inucl + nucl_num], nucl_coord[inucl + 2*nucl_num] }; - /* Test if the point is in the range of the nucleus */ + const double shift = 1.e-20; + const double x = e_coord[0] - n_coord[0]; const double y = e_coord[1] - n_coord[1]; const double z = e_coord[2] - n_coord[2]; @@ -6609,6 +6616,7 @@ qmckl_compute_ao_vgl_hpc_gaussian ( const double z2 = z*z; const double r2 = x2 + y2 + z2; + /* Test if the point is in the range of the nucleus */ if (r2 > cutoff * nucleus_range[inucl]) { continue; } @@ -6619,21 +6627,21 @@ qmckl_compute_ao_vgl_hpc_gaussian ( break; case 1: - poly_vgl_l1[0][1] = x; - poly_vgl_l1[0][2] = y; - poly_vgl_l1[0][3] = z; + poly_vgl_l1[0][1] = x+shift; + poly_vgl_l1[0][2] = y+shift; + poly_vgl_l1[0][3] = z+shift; break; case 2: - poly_vgl_l2[0][1] = x; - poly_vgl_l2[0][2] = y; - poly_vgl_l2[0][3] = z; - poly_vgl_l2[0][4] = x2; - poly_vgl_l2[0][5] = x*y; - poly_vgl_l2[0][6] = x*z; - poly_vgl_l2[0][7] = y2; - poly_vgl_l2[0][8] = y*z; - poly_vgl_l2[0][9] = z2; + poly_vgl_l2[0][1] = x+shift; + poly_vgl_l2[0][2] = y+shift; + poly_vgl_l2[0][3] = z+shift; + poly_vgl_l2[0][4] = x2+shift; + poly_vgl_l2[0][5] = x*y+shift; + poly_vgl_l2[0][6] = x*z+shift; + poly_vgl_l2[0][7] = y2+shift; + poly_vgl_l2[0][8] = y*z+shift; + poly_vgl_l2[0][9] = z2+shift; poly_vgl_l2[1][4] = x+x; poly_vgl_l2[1][5] = y; poly_vgl_l2[1][6] = z; @@ -6772,7 +6780,7 @@ qmckl_compute_ao_vgl_hpc_gaussian ( ao_vgl_4[il] = 0.0; ao_vgl_5[il] = 0.0; } -*/ +,*/ continue; } @@ -6876,7 +6884,7 @@ qmckl_compute_ao_vgl_hpc_gaussian ( ao_vgl_4[il] = 0.0; ao_vgl_5[il] = 0.0; } -*/ +,*/ } } } diff --git a/org/qmckl_mo.org b/org/qmckl_mo.org index fafef30..33f5432 100644 --- a/org/qmckl_mo.org +++ b/org/qmckl_mo.org @@ -1825,7 +1825,8 @@ integer function qmckl_compute_mo_basis_mo_vgl_doc_f(context, & end do end do -! info = qmckl_dgemm(context,'N', 'N', mo_num, point_num, ao_num, 1.d0, & +! print *, 'coucou' +! info = qmckl_dgemm(context,'N', 'N', mo_num, point_num*5_8, ao_num, 1.d0, & ! coefficient_t, int(size(coefficient_t,1),8), & ! ao_vgl, int(size(ao_vgl,1),8), 0.d0, & ! mo_vgl, int(size(mo_vgl,1),8))