From 1c30e3e207ebc2fafde56d174217a4a0a1188228 Mon Sep 17 00:00:00 2001 From: Emiel Slootman Date: Wed, 16 Apr 2025 23:43:47 +0200 Subject: [PATCH] test for ao_hessian: failing --- org/qmckl_ao.org | 122 ++++++++++++++++++++++++++++++++++++++++++- org/qmckl_forces.org | 2 +- 2 files changed, 122 insertions(+), 2 deletions(-) diff --git a/org/qmckl_ao.org b/org/qmckl_ao.org index 438828d..d0a8d89 100644 --- a/org/qmckl_ao.org +++ b/org/qmckl_ao.org @@ -8229,7 +8229,7 @@ function qmckl_compute_ao_hessian_doc(context, & e_coord(1) = coord(ipoint,1) e_coord(2) = coord(ipoint,2) e_coord(3) = coord(ipoint,3) - ao_hessian(:,:,ipoint,:,:) = 0.d0 + !ao_hessian(:,:,ipoint,:,:) = 0.d0 do inucl=1,nucl_num n_coord(1) = nucl_coord(inucl,1) n_coord(2) = nucl_coord(inucl,2) @@ -8318,6 +8318,126 @@ end function qmckl_compute_ao_hessian_doc double* const ao_hessian ); #+end_src +*** Test :noexport: + + #+begin_src c :tangle (eval c_test) :exports none +printf("AO hessian\n"); + +int64_t point_num = elec_num; + +double delta_x = 0.00001; +double coef[9] = { 1.0/280.0, -4.0/105.0, 1.0/5.0, -4.0/5.0, 0.0, 4.0/5.0, -1.0/5.0, 4.0/105.0, -1.0/280.0 }; + +rc = qmckl_set_nucleus_coord(context, 'T', &(nucl_coord[0]), 3*nucl_num); +assert(rc == QMCKL_SUCCESS); + +double * ao_hessian = (double*) malloc(ao_num*4*point_num*nucl_num*3 *sizeof(double)); +rc = qmckl_get_ao_basis_ao_hessian(context, &ao_hessian[0], 3*4*nucl_num*ao_num*point_num); +assert(rc == QMCKL_SUCCESS); + +double * finite_difference_ao_hessian= (double*) malloc(3*nucl_num* 3 * point_num * ao_num * sizeof(double)); + +double* nucleus_coord = (double*) malloc(3 * nucl_num * sizeof(double)); +if (nucleus_coord == NULL) { + return QMCKL_ALLOCATION_FAILED; +} + +rc = qmckl_get_nucleus_coord(context, 'N', nucleus_coord, 3*nucl_num); + +double* temp_coord = (double*) malloc(3 * nucl_num * sizeof(double)); +if (temp_coord == NULL) { + free(nucleus_coord); + return QMCKL_ALLOCATION_FAILED; +} + +double* ao_vgl= (double*) malloc(5 * point_num * ao_num * sizeof(double)); +if (ao_vgl == NULL) { + free(temp_coord); + free(nucleus_coord); + return QMCKL_ALLOCATION_FAILED; +} + + +// Copy original coordinates +for (int i = 0; i < 3 * nucl_num; i++) { + temp_coord[i] = nucleus_coord[i]; +} + + +for (int64_t a = 0; a < nucl_num; a++) { + for (int64_t k = 0; k < 3; k++) { + for (int64_t m = -4; m <= 4; m++) { + + // Apply finite difference displacement + temp_coord[k+a*3] = nucleus_coord[k+3*a] + (double) m * delta_x; + + // Update coordinates in the context + rc = qmckl_set_nucleus_coord(context, 'N', temp_coord, 3*nucl_num); + assert(rc == QMCKL_SUCCESS); + + rc = qmckl_context_touch(context); + assert(rc == QMCKL_SUCCESS); + + // Call the provided function + rc = qmckl_get_ao_basis_ao_vgl(context,&ao_vgl[0], 5*point_num*ao_num); + assert(rc == QMCKL_SUCCESS); + + // Accumulate derivative using finite-difference coefficients + for (int i = 0; i < point_num; i++) { + for (int n = 0; n < 3; n++){ + for (int j = 0; j < ao_num; j++) { + if (m == -4) { + finite_difference_ao_hessian[k*3*ao_num*point_num*nucl_num + a*3*ao_num*point_num + i*3*ao_num + n*ao_num + j] = 0.0; + } + finite_difference_ao_hessian[k*3*ao_num*point_num*nucl_num + a*3*ao_num*point_num + i*3*ao_num + n*ao_num + j] -= + coef[m + 4] * ao_vgl[i*ao_num*5 + (n+1)*ao_num + j]/delta_x; + } + } + } + } + temp_coord[k+a*3] = nucleus_coord[k+3*a]; + } +} + +// Reset coordinates in the context +rc = qmckl_set_nucleus_coord(context, 'N', temp_coord, 3*nucl_num); +assert(rc == QMCKL_SUCCESS); + +rc = qmckl_context_touch(context); +assert(rc == QMCKL_SUCCESS); + +free(nucleus_coord); +free(temp_coord); +free(ao_vgl); + + + +for (int a = 0; a < nucl_num; a++) { + for (int n = 0; n < 3; n++){ + for (int j = 0; j < point_num; j++){ + for (int m = 0; m < 3; m++){ + for (int i = 0; i < ao_num; i++){ + printf("a=%i n=%i j=%i m=%i i=%i\n", a, n, j, m, i); + printf("%.10f\n", ao_hessian[a*12*ao_num*point_num + n*3*ao_num*point_num + j*3*ao_num + m*ao_num + i]); + printf("%.10f\n", finite_difference_ao_hessian[n*3*ao_num*point_num*nucl_num + a*3*ao_num*point_num + j*3*ao_num + m*ao_num + i]); + assert(fabs(finite_difference_ao_hessian[n*3*ao_num*point_num*nucl_num + + a*3*ao_num*point_num + j*3*ao_num + + m*ao_num + i] - + ao_hessian[a*12*ao_num*point_num + n*3*ao_num*point_num + + j*3*ao_num + m*ao_num + i]) < 1.e-9); + } + } + } + } +} + +free(ao_hessian); +free(finite_difference_ao_hessian); + +printf("OK\n"); + + #+end_src + * End of files :noexport: diff --git a/org/qmckl_forces.org b/org/qmckl_forces.org index 47b99ae..3ae7c42 100644 --- a/org/qmckl_forces.org +++ b/org/qmckl_forces.org @@ -7006,7 +7006,7 @@ for (int a = 0; a < nucl_num; a++) { a*3*mo_num*point_num + j*3*mo_num + m*mo_num + i] - forces_mo_g[a*9*mo_num*point_num + n*3*mo_num*point_num + - j*3*mo_num + m*mo_num + i]) < 1.e-8); + j*3*mo_num + m*mo_num + i]) < 1.e-9); } } }