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

Fix wrong gradient at nodes of AOs

Browse Source
This commit is contained in:
Anthony Scemama 2024-05-08 13:56:32 +02:00
parent 8c5ec872ed
commit 9553059bbe
2 changed files with 50 additions and 41 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

88
org/qmckl_ao.org
View File

@ -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<ao_num ; ++k) {
ptr[k] = 0.;
}
*/
,*/
memset(&ao_value[ipoint*ao_num], 0, ao_num*sizeof(double));
@ -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];
@ -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;
}
@ -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;
}
*/
,*/
}
}
}

3
org/qmckl_mo.org
View File

@ -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))