test for ao_hessian: failing

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:
 

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);
         }
       }
     }