Skip to content

Commit

Permalink
Fix wrong gradient at nodes of AOs
Browse files Browse the repository at this point in the history
  • Loading branch information
scemama committed May 8, 2024
1 parent 8c5ec87 commit 9553059
Show file tree
Hide file tree
Showing 2 changed files with 50 additions and 41 deletions.
88 changes: 48 additions & 40 deletions org/qmckl_ao.org
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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,
Expand All @@ -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;
Expand Down Expand Up @@ -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,
Expand All @@ -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];
Expand Down Expand Up @@ -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];
Expand Down Expand Up @@ -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;
Expand All @@ -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;
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -5793,7 +5796,7 @@ qmckl_compute_ao_value_hpc_gaussian (const qmckl_context context,
for (size_t k=0 ; k<ao_num ; ++k) {
ptr[k] = 0.;
}
*/
,*/

memset(&ao_value[ipoint*ao_num], 0, ao_num*sizeof(double));

Expand All @@ -5808,7 +5811,9 @@ qmckl_compute_ao_value_hpc_gaussian (const qmckl_context context,
nucl_coord[inucl + nucl_num],
nucl_coord[inucl + 2*nucl_num] };

/* Test if the point is in the range of the nucleus */
/* Shift to avoid haing exact zeros at the nodes of the AOs */
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];
Expand All @@ -5818,6 +5823,7 @@ qmckl_compute_ao_value_hpc_gaussian (const qmckl_context context,
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;
}
Expand All @@ -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;
Expand Down Expand Up @@ -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));

Expand All @@ -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];
Expand All @@ -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;
}
Expand All @@ -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;
Expand Down Expand Up @@ -6772,7 +6780,7 @@ qmckl_compute_ao_vgl_hpc_gaussian (
ao_vgl_4[il] = 0.0;
ao_vgl_5[il] = 0.0;
}
*/
,*/
continue;
}

Expand Down Expand Up @@ -6876,7 +6884,7 @@ qmckl_compute_ao_vgl_hpc_gaussian (
ao_vgl_4[il] = 0.0;
ao_vgl_5[il] = 0.0;
}
*/
,*/
}
}
}
Expand Down
3 changes: 2 additions & 1 deletion org/qmckl_mo.org
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand Down

0 comments on commit 9553059

Please sign in to comment.