From d257e28b92ac944aafdd371735f9312c1509cecf Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Thu, 11 Jan 2024 14:33:40 +0100
Subject: [PATCH] Fix bug in HPC laplacian AO

---
 org/qmckl_ao.org     | 30 ++++++++++++++++++------------
 org/qmckl_mo.org     |  1 +
 org/qmckl_trexio.org |  2 +-
 3 files changed, 20 insertions(+), 13 deletions(-)

diff --git a/org/qmckl_ao.org b/org/qmckl_ao.org
index e215e8f..aa59f38 100644
--- a/org/qmckl_ao.org
+++ b/org/qmckl_ao.org
@@ -114,6 +114,12 @@ int main() {
 
   #+end_src
 
+  #+NAME:fortran_cutoff
+  #+begin_src fortran
+ cutoff = qmckl_get_numprec_precision(context)
+ cutoff = dlog(2.d0**(cutoff-2))
+  #+end_src
+
 * Context
 
 ** Constant data
@@ -3605,7 +3611,7 @@ function qmckl_compute_ao_basis_primitive_gaussian_vgl &
 
   use qmckl_constants
 
-  use qmckl, only: qmckl_get_numprec_epsilon
+  use qmckl, only: qmckl_get_numprec_precision
   implicit none
   integer (qmckl_context), intent(in)  , value :: context
   integer (c_int64_t) , intent(in)  , value :: prim_num
@@ -3624,7 +3630,7 @@ function qmckl_compute_ao_basis_primitive_gaussian_vgl &
   info = QMCKL_SUCCESS
 
   ! Don't compute exponentials when the result will be almost zero.
-  cutoff = -dlog(qmckl_get_numprec_epsilon(context))
+  <<fortran_cutoff>>
 
   do inucl=1,nucl_num
      ! C is zero-based, so shift bounds by one
@@ -3908,7 +3914,7 @@ function qmckl_compute_ao_basis_shell_gaussian_vgl( &
 
   use, intrinsic :: iso_c_binding
   use qmckl_constants
-  use qmckl, only: qmckl_get_numprec_epsilon
+  use qmckl, only: qmckl_get_numprec_precision
 
   implicit none
   integer (qmckl_context), intent(in)  , value :: context
@@ -3936,7 +3942,7 @@ function qmckl_compute_ao_basis_shell_gaussian_vgl( &
   info = QMCKL_SUCCESS
 
   ! Don't compute exponentials when the result will be almost zero.
-  cutoff = -dlog(qmckl_get_numprec_epsilon(context))
+  <<fortran_cutoff>>
 
   do ipoint = 1, point_num
 
@@ -4358,7 +4364,7 @@ function test_qmckl_ao_power(context) bind(C)
   integer*8                     :: i,j
   double precision              :: epsilon
 
-  epsilon = qmckl_get_numprec_epsilon(context)
+  epsilon = qmckl_get_numprec_precision(context)
 
   n = 100;
   LDP = 10;
@@ -5247,7 +5253,7 @@ function test_qmckl_ao_polynomial_vgl(context) bind(C)
   double precision              :: epsilon
   integer(qmckl_exit_code) ::  test_qmckl_ao_polynomial_vgl
 
-  epsilon = qmckl_get_numprec_epsilon(context)
+  epsilon = qmckl_get_numprec_precision(context)
 
   X = (/ 1.1 , 2.2 ,  3.3 /)
   R = (/ 0.1 , 1.2 , -2.3 /)
@@ -5405,7 +5411,7 @@ function qmckl_compute_nucleus_range_gaussian(context, &
      bind(C) result(info)
   use, intrinsic :: iso_c_binding
   use qmckl_constants
-  use qmckl, only: qmckl_ao_polynomial_vgl, qmckl_get_numprec_epsilon
+  use qmckl, only: qmckl_ao_polynomial_vgl, qmckl_get_numprec_precision
   implicit none
   integer (qmckl_context), intent(in)  , value :: context
   integer (c_int64_t) , intent(in)  , value :: ao_num
@@ -5566,7 +5572,7 @@ function qmckl_compute_ao_value_doc(context, &
      bind(C) result(info)
   use, intrinsic :: iso_c_binding
   use qmckl_constants
-  use qmckl, only: qmckl_ao_polynomial_vgl, qmckl_get_numprec_epsilon
+  use qmckl, only: qmckl_ao_polynomial_vgl, qmckl_get_numprec_precision
   implicit none
   integer (qmckl_context), intent(in)  , value :: context
   integer (c_int64_t) , intent(in)  , value :: ao_num
@@ -5618,7 +5624,7 @@ function qmckl_compute_ao_value_doc(context, &
   info = QMCKL_SUCCESS
 
   ! Don't compute polynomials when the radial part is zero.
-  cutoff = -dlog(qmckl_get_numprec_epsilon(context))
+  <<fortran_cutoff>>
 
   do ipoint = 1, point_num
      e_coord(1) = coord(ipoint,1)
@@ -6354,7 +6360,7 @@ function qmckl_compute_ao_vgl_doc(context, &
      bind(C) result(info)
   use, intrinsic :: iso_c_binding
   use qmckl_constants
-  use qmckl, only : qmckl_ao_polynomial_vgl, qmckl_get_numprec_epsilon
+  use qmckl, only : qmckl_ao_polynomial_vgl, qmckl_get_numprec_precision
   implicit none
   integer (qmckl_context), intent(in)  , value :: context
   integer (c_int64_t) , intent(in)  , value :: ao_num
@@ -6405,7 +6411,7 @@ function qmckl_compute_ao_vgl_doc(context, &
   info = QMCKL_SUCCESS
 
   ! Don't compute polynomials when the radial part is zero.
-  cutoff = -dlog(qmckl_get_numprec_epsilon(context))
+  <<fortran_cutoff>>
 
   do ipoint = 1, point_num
      e_coord(1) = coord(ipoint,1)
@@ -6805,7 +6811,7 @@ IVDEP
               ao_vgl_2[0] = s2 * f[0];
               ao_vgl_3[0] = s3 * f[0];
               ao_vgl_4[0] = s4 * f[0];
-              ao_vgl_5[0] = s5;
+              ao_vgl_5[0] = s5 * f[0];
               break;
             case 3:
 IVDEP
diff --git a/org/qmckl_mo.org b/org/qmckl_mo.org
index 85f0607..6599d98 100644
--- a/org/qmckl_mo.org
+++ b/org/qmckl_mo.org
@@ -119,6 +119,7 @@ int main() {
 #include <string.h>
 #include <stdbool.h>
 #include <assert.h>
+#include <math.h>
 
 #include "qmckl.h"
 #include "qmckl_context_private_type.h"
diff --git a/org/qmckl_trexio.org b/org/qmckl_trexio.org
index 8ca99f1..528907d 100644
--- a/org/qmckl_trexio.org
+++ b/org/qmckl_trexio.org
@@ -368,7 +368,7 @@ qmckl_trexio_read_ao_X(qmckl_context context, trexio_t* const file)
 *** Number of atomic orbitals
 
    #+begin_src c :tangle (eval c)
-  int64_t ao_num = 0LL;
+  int64_t ao_num = 0;
 
   rcio = trexio_read_ao_num_64(file, &ao_num);
   if (rcio != TREXIO_SUCCESS) {