From c0131d5da48eede9cc066107890e5e57ed24a648 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 14 Mar 2023 19:13:49 +0100 Subject: [PATCH] Improved en_distance --- org/qmckl_blas.org | 30 +++-- org/qmckl_electron.org | 111 ++++++++-------- org/qmckl_mo.org | 296 +++++++++++++++++++++++++++++++++-------- org/qmckl_trexio.org | 2 +- 4 files changed, 317 insertions(+), 122 deletions(-) diff --git a/org/qmckl_blas.org b/org/qmckl_blas.org index 5e566d2..ee878fe 100644 --- a/org/qmckl_blas.org +++ b/org/qmckl_blas.org @@ -248,7 +248,7 @@ qmckl_matrix_free( qmckl_context context, | Variable | Type | Description | |----------+-----------------------------------+-----------------------------| - | ~order~ | ~int64_t~ | Order of the tensor | + | ~order~ | ~int32_t~ | Order of the tensor | | ~size~ | ~int64_t[QMCKL_TENSOR_ORDER_MAX]~ | Dimension of each component | | ~data~ | ~double*~ | Elements | @@ -260,7 +260,8 @@ qmckl_matrix_free( qmckl_context context, typedef struct qmckl_tensor { double* restrict data; - int64_t order; + int32_t order; + int32_t __space__; int64_t size[QMCKL_TENSOR_ORDER_MAX]; } qmckl_tensor; #+end_src @@ -269,7 +270,7 @@ typedef struct qmckl_tensor { #+begin_src c :comments org :tangle (eval h_private_func) qmckl_tensor qmckl_tensor_alloc( qmckl_context context, - const int64_t order, + const int32_t order, const int64_t* size); #+end_src @@ -279,7 +280,7 @@ qmckl_tensor_alloc( qmckl_context context, #+begin_src c :comments org :tangle (eval c) :exports none qmckl_tensor qmckl_tensor_alloc( qmckl_context context, - const int64_t order, + const int32_t order, const int64_t* size) { /* Should always be true by contruction */ @@ -288,10 +289,11 @@ qmckl_tensor_alloc( qmckl_context context, assert (size != NULL); qmckl_tensor result; + memset(&result, 0, sizeof(qmckl_tensor)); result.order = order; int64_t prod_size = (int64_t) 1; - for (int64_t i=0 ; i (int64_t) 0); result.size[i] = size[i]; prod_size *= size[i]; @@ -383,7 +385,7 @@ qmckl_matrix_of_vector(const qmckl_vector vector, #+begin_src c :comments org :tangle (eval h_private_func) qmckl_tensor qmckl_tensor_of_vector(const qmckl_vector vector, - const int64_t order, + const int32_t order, const int64_t* size); #+end_src @@ -392,13 +394,13 @@ qmckl_tensor_of_vector(const qmckl_vector vector, #+begin_src c :comments org :tangle (eval c) :exports none qmckl_tensor qmckl_tensor_of_vector(const qmckl_vector vector, - const int64_t order, + const int32_t order, const int64_t* size) { qmckl_tensor result; int64_t prod_size = 1; - for (int64_t i=0 ; ielectron.walker.point.date > ctx->electron.ee_distance_date) { + if (ctx->point.date > ctx->electron.ee_distance_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { free(ctx->electron.ee_distance); @@ -965,7 +967,7 @@ qmckl_exit_code qmckl_provide_ee_potential(qmckl_context context) if (rc != QMCKL_SUCCESS) return rc; /* Compute if necessary */ - if (ctx->electron.walker.point.date > ctx->electron.ee_potential_date) { + if (ctx->point.date > ctx->electron.ee_potential_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { free(ctx->electron.ee_potential); @@ -1144,7 +1146,7 @@ qmckl_exit_code qmckl_get_electron_en_distance(qmckl_context context, double* di qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); - size_t sze = ctx->electron.num * ctx->nucleus.num * ctx->electron.walker.num; + size_t sze = ctx->point.num * ctx->nucleus.num; memcpy(distance, ctx->electron.en_distance, sze * sizeof(double)); return QMCKL_SUCCESS; @@ -1176,19 +1178,28 @@ qmckl_exit_code qmckl_provide_en_distance(qmckl_context context) } /* Compute if necessary */ - if (ctx->electron.walker.point.date > ctx->electron.en_distance_date) { + if (ctx->point.date > ctx->electron.en_distance_date) { - if (ctx->electron.walker.num > ctx->electron.walker_old.num) { - free(ctx->electron.en_distance); - ctx->electron.en_distance = NULL; + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->point.num * ctx->nucleus.num * sizeof(double); + + if (ctx->electron.en_distance != NULL) { + qmckl_memory_info_struct mem_info_test = qmckl_memory_info_struct_zero; + qmckl_exit_code rc = qmckl_get_malloc_info(context, ctx->electron.en_distance, &mem_info_test); + + /* if rc != QMCKL_SUCCESS, we are maybe in an _inplace function because the + memory was not allocated with qmckl_malloc */ + + if ((rc == QMCKL_SUCCESS) && (mem_info_test.size != mem_info.size)) { + rc = qmckl_free(context, ctx->electron.en_distance); + assert (rc == QMCKL_SUCCESS); + ctx->electron.en_distance = NULL; + } } /* Allocate array */ if (ctx->electron.en_distance == NULL) { - qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; - mem_info.size = ctx->electron.num * ctx->nucleus.num * - ctx->electron.walker.num * sizeof(double); double* en_distance = (double*) qmckl_malloc(context, mem_info); if (en_distance == NULL) { @@ -1202,10 +1213,9 @@ qmckl_exit_code qmckl_provide_en_distance(qmckl_context context) qmckl_exit_code rc = qmckl_compute_en_distance(context, - ctx->electron.num, + ctx->point.num, ctx->nucleus.num, - ctx->electron.walker.num, - ctx->electron.walker.point.coord.data, + ctx->point.coord.data, ctx->nucleus.coord.data, ctx->electron.en_distance); if (rc != QMCKL_SUCCESS) { @@ -1227,28 +1237,26 @@ qmckl_exit_code qmckl_provide_en_distance(qmckl_context context) :END: #+NAME: qmckl_en_distance_args - | Variable | Type | In/Out | Description | - |---------------+----------------------------------------+--------+----------------------------| - | ~context~ | ~qmckl_context~ | in | Global state | - | ~elec_num~ | ~int64_t~ | in | Number of electrons | - | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | - | ~walk_num~ | ~int64_t~ | in | Number of walkers | - | ~elec_coord~ | ~double[3][walk_num][elec_num]~ | in | Electron coordinates | - | ~nucl_coord~ | ~double[3][elec_num]~ | in | Nuclear coordinates | - | ~en_distance~ | ~double[walk_num][elec_num][nucl_num]~ | out | Electron-nucleus distances | + | Variable | Type | In/Out | Description | + |---------------+-------------------------------+--------+----------------------------| + | ~context~ | ~qmckl_context~ | in | Global state | + | ~point_num~ | ~int64_t~ | in | Number of points | + | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | + | ~elec_coord~ | ~double[3][point_num]~ | in | Electron coordinates | + | ~nucl_coord~ | ~double[3][nucl_num]~ | in | Nuclear coordinates | + | ~en_distance~ | ~double[point_num][nucl_num]~ | out | Electron-nucleus distances | #+begin_src f90 :comments org :tangle (eval f) :noweb yes -integer function qmckl_compute_en_distance_f(context, elec_num, nucl_num, walk_num, elec_coord, nucl_coord, en_distance) & +integer function qmckl_compute_en_distance_f(context, point_num, nucl_num, elec_coord, nucl_coord, en_distance) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context - integer*8 , intent(in) :: elec_num + integer*8 , intent(in) :: point_num integer*8 , intent(in) :: nucl_num - integer*8 , intent(in) :: walk_num - double precision , intent(in) :: elec_coord(elec_num,walk_num,3) + double precision , intent(in) :: elec_coord(point_num,3) double precision , intent(in) :: nucl_coord(nucl_num,3) - double precision , intent(out) :: en_distance(nucl_num,elec_num,walk_num) + double precision , intent(out) :: en_distance(nucl_num,point_num) integer*8 :: k @@ -1259,7 +1267,7 @@ integer function qmckl_compute_en_distance_f(context, elec_num, nucl_num, walk_n return endif - if (elec_num <= 0) then + if (point_num <= 0) then info = QMCKL_INVALID_ARG_2 return endif @@ -1269,14 +1277,9 @@ integer function qmckl_compute_en_distance_f(context, elec_num, nucl_num, walk_n return endif - if (walk_num <= 0) then - info = QMCKL_INVALID_ARG_4 - return - endif - - info = qmckl_distance(context, 'T', 'T', nucl_num, elec_num * walk_num, & + info = qmckl_distance(context, 'T', 'T', nucl_num, point_num, & nucl_coord, nucl_num, & - elec_coord, elec_num * walk_num, & + elec_coord, point_num, & en_distance, nucl_num) end function qmckl_compute_en_distance_f @@ -1285,9 +1288,8 @@ end function qmckl_compute_en_distance_f #+begin_src c :tangle (eval h_private_func) :comments org :exports none qmckl_exit_code qmckl_compute_en_distance ( const qmckl_context context, - const int64_t elec_num, + const int64_t point_num, const int64_t nucl_num, - const int64_t walk_num, const double* elec_coord, const double* nucl_coord, double* const en_distance ); @@ -1298,23 +1300,22 @@ qmckl_exit_code qmckl_compute_en_distance ( #+RESULTS: #+begin_src f90 :tangle (eval f) :comments org :exports none integer(c_int32_t) function qmckl_compute_en_distance & - (context, elec_num, nucl_num, walk_num, elec_coord, nucl_coord, en_distance) & + (context, point_num, nucl_num, elec_coord, nucl_coord, en_distance) & bind(C) result(info) use, intrinsic :: iso_c_binding implicit none integer (c_int64_t) , intent(in) , value :: context - integer (c_int64_t) , intent(in) , value :: elec_num + integer (c_int64_t) , intent(in) , value :: point_num integer (c_int64_t) , intent(in) , value :: nucl_num - integer (c_int64_t) , intent(in) , value :: walk_num - real (c_double ) , intent(in) :: elec_coord(elec_num,walk_num,3) - real (c_double ) , intent(in) :: nucl_coord(elec_num,3) - real (c_double ) , intent(out) :: en_distance(nucl_num,elec_num,walk_num) + real (c_double ) , intent(in) :: elec_coord(point_num,3) + real (c_double ) , intent(in) :: nucl_coord(nucl_num,3) + real (c_double ) , intent(out) :: en_distance(nucl_num,point_num) integer(c_int32_t), external :: qmckl_compute_en_distance_f info = qmckl_compute_en_distance_f & - (context, elec_num, nucl_num, walk_num, elec_coord, nucl_coord, en_distance) + (context, point_num, nucl_num, elec_coord, nucl_coord, en_distance) end function qmckl_compute_en_distance #+end_src @@ -1451,7 +1452,7 @@ qmckl_exit_code qmckl_provide_en_potential(qmckl_context context) if (rc != QMCKL_SUCCESS) return rc; /* Compute if necessary */ - if (ctx->electron.walker.point.date > ctx->electron.en_potential_date) { + if (ctx->point.date > ctx->electron.en_potential_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { free(ctx->electron.en_potential); diff --git a/org/qmckl_mo.org b/org/qmckl_mo.org index e2df46a..bf5b1ab 100644 --- a/org/qmckl_mo.org +++ b/org/qmckl_mo.org @@ -157,7 +157,7 @@ typedef struct qmckl_mo_basis_struct { double * restrict mo_vgl; double * restrict mo_value; - double * restrict r_cusp_param; + qmckl_tensor r_cusp_param; uint64_t mo_vgl_date; uint64_t mo_value_date; @@ -188,6 +188,8 @@ qmckl_exit_code qmckl_init_mo_basis(qmckl_context context) { qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); + ctx->mo_basis.r_cusp = NULL; + ctx->mo_basis.uninitialized = (1 << 2) - 1; return QMCKL_SUCCESS; @@ -252,7 +254,9 @@ qmckl_exit_code qmckl_set_mo_basis_mo_num(qmckl_context context, const int64_t m return QMCKL_SUCCESS; } + #+end_src + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_set_mo_basis_coefficient(qmckl_context context, const double* coefficient, const int64_t size_max) @@ -292,55 +296,6 @@ qmckl_exit_code qmckl_set_mo_basis_coefficient(qmckl_context context, <> } - - -qmckl_exit_code -qmckl_set_mo_basis_r_cusp(qmckl_context context, - const double* r_cusp, - const int64_t size_max) -{ - - int32_t mask = 0; - - <
>
-
-  if (r_cusp == NULL) {
-      return qmckl_failwith( context, QMCKL_INVALID_ARG_2,
-                             "qmckl_set_mo_basis_r_cusp",
-                             "r_cusp: Null pointer");
-  }
-        
-  if (size_max < ctx->nucleus.num) {
-    return qmckl_failwith( context, QMCKL_INVALID_ARG_3,
-                             "qmckl_set_mo_basis_r_cusp",
-                             "Array too small");
-  }
-  
-
-  if (ctx->mo_basis.r_cusp != NULL) {
-    qmckl_exit_code rc = qmckl_free(context, ctx->mo_basis.r_cusp);
-    if (rc != QMCKL_SUCCESS) {
-      return rc;
-    }
-  }
-
-  qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
-  mem_info.size = ctx->nucleus.num * sizeof(double);
-  double* new_array = (double*) qmckl_malloc(context, mem_info);
-  if (new_array == NULL) {
-    return qmckl_failwith( context,
-                           QMCKL_ALLOCATION_FAILED,
-                           "qmckl_set_mo_basis_r_cusp",
-                           NULL);
-  }
-
-  memcpy(new_array, r_cusp, mem_info.size);
-
-  ctx->mo_basis.r_cusp = new_array;
-  return QMCKL_SUCCESS;
-
-}
-
    #+end_src
 
  When the basis set is completely entered, other data structures are
@@ -406,10 +361,246 @@ qmckl_exit_code qmckl_finalize_mo_basis(qmckl_context context) {
     ctx->mo_basis.mo_value = NULL;
     ctx->mo_basis.mo_value_date = 0;
   }
+
   return qmckl_context_touch(context);
 }
    #+end_src
 
+** Cusp adjsutment functions
+
+   To activate the cusp adjustment, the user must enter the radius of
+   the fitting function for each atom.
+   
+   This function requires the computation of the value and gradients
+   of the $s$ AOs at the distance equal to the radius, and the values
+   of the non-$s$ AOs at the center.
+   
+   #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
+qmckl_exit_code
+qmckl_set_mo_basis_r_cusp(qmckl_context context,
+                          const double* r_cusp,
+                          const int64_t size_max)
+{
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return QMCKL_NULL_CONTEXT;
+  }
+
+  qmckl_exit_code rc;
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
+  
+  if (r_cusp == NULL) {
+      return qmckl_failwith( context, QMCKL_INVALID_ARG_2,
+                             "qmckl_set_mo_basis_r_cusp",
+                             "r_cusp: Null pointer");
+  }
+        
+  if (size_max < ctx->nucleus.num) {
+    return qmckl_failwith( context, QMCKL_INVALID_ARG_3,
+                             "qmckl_set_mo_basis_r_cusp",
+                             "Array too small");
+  }
+  
+
+  // Nullify r_cusp
+  if (ctx->mo_basis.r_cusp != NULL) {
+    rc = qmckl_free(context, ctx->mo_basis.r_cusp);
+    if (rc != QMCKL_SUCCESS) {
+      return rc;
+    }
+    ctx->mo_basis.r_cusp = NULL;
+  }
+
+
+  // Save old points
+  int64_t old_point_num = ctx->point.num;
+  double* old_coord = NULL;
+  if (old_point_num > 0) {
+    qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
+    mem_info.size = old_point_num * 3 * sizeof(double);
+    old_coord = (double*) qmckl_malloc(context, mem_info);
+    
+    rc = qmckl_get_point(context, 'T', old_coord,  (old_point_num * 3));
+    if (rc != QMCKL_SUCCESS) {
+      return qmckl_failwith( context,
+                             QMCKL_FAILURE,
+                             "qmckl_set_mo_basis_r_cusp",
+                             "Unable to get old coordinates");
+    }
+  }
+
+  double* coord = (double*) malloc(ctx->nucleus.num * 3 * sizeof(double));
+
+  // Evaluate MO vgl at r_cusp
+  qmckl_tensor mo_vgl_at_r_cusp;
+  {
+    int64_t sze[3] = { ctx->mo_basis.mo_num, 5, ctx->nucleus.num };
+    mo_vgl_at_r_cusp = qmckl_tensor_alloc(context, 3, &(sze[0]));
+    
+    rc = qmckl_double_of_matrix(context, ctx->nucleus.coord, coord, ctx->nucleus.num * 3);
+    if (rc != QMCKL_SUCCESS) return rc;
+
+    int64_t i=0;
+    for (int64_t inucl=0 ; inuclnucleus.num ; ++inucl) {
+      coord[i] += r_cusp[inucl];
+      i+=3;
+    }
+      
+    rc = qmckl_set_point(context, 'T', ctx->nucleus.num, coord, (ctx->nucleus.num * 3));
+    if (rc != QMCKL_SUCCESS) {
+      return qmckl_failwith( context,
+                             QMCKL_FAILURE,
+                             "qmckl_set_mo_basis_r_cusp",
+                             "Unable to set coordinates at r_cusp");
+    }
+    
+    rc = qmckl_get_mo_basis_mo_vgl(context,
+                                   &(qmckl_ten3(mo_vgl_at_r_cusp,0,0,0)), 
+                                   ctx->mo_basis.mo_num * 5 * ctx->nucleus.num);
+    if (rc != QMCKL_SUCCESS) return rc;
+    
+  }
+
+  // Set r_cusp
+  {
+    assert (ctx->mo_basis.r_cusp == NULL);
+    qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
+    mem_info.size = ctx->nucleus.num * sizeof(double);
+    ctx->mo_basis.r_cusp = (double*) qmckl_malloc(context, mem_info);
+    if (ctx->mo_basis.r_cusp == NULL) {
+      return qmckl_failwith( context,
+                             QMCKL_ALLOCATION_FAILED,
+                             "qmckl_set_mo_basis_r_cusp",
+                             NULL);
+    }
+    memcpy(ctx->mo_basis.r_cusp, r_cusp, mem_info.size);
+  }
+
+
+  
+  // Allocate cusp parameters and set them to zero
+  {
+    if (ctx->mo_basis.r_cusp_param.size[0] != 0) {
+      rc = qmckl_tensor_free(context, &(ctx->mo_basis.r_cusp_param));
+      if (rc != QMCKL_SUCCESS) return rc;
+    }
+    
+    int64_t sze[3] = { 4, ctx->mo_basis.mo_num, ctx->nucleus.num };
+    ctx->mo_basis.r_cusp_param = qmckl_tensor_alloc(context, 3, &(sze[0]));
+    ctx->mo_basis.r_cusp_param = qmckl_tensor_set(ctx->mo_basis.r_cusp_param, 0.);
+  }
+
+  
+  // Evaluate MO value at nucleus without s components
+  qmckl_matrix mo_value_at_nucl_no_s;
+  {
+    mo_value_at_nucl_no_s = qmckl_matrix_alloc(context, ctx->mo_basis.mo_num, ctx->nucleus.num);
+    
+    rc = qmckl_double_of_matrix(context, ctx->nucleus.coord, coord, ctx->nucleus.num * 3);
+    if (rc != QMCKL_SUCCESS) return rc;
+
+    rc = qmckl_set_point(context, 'T', ctx->nucleus.num, coord, (ctx->nucleus.num * 3));
+    if (rc != QMCKL_SUCCESS) {
+      return qmckl_failwith( context,
+                             QMCKL_FAILURE,
+                             "qmckl_set_mo_basis_r_cusp",
+                             "Unable to set coordinates at the nuclei");
+    }
+    
+    rc = qmckl_get_mo_basis_mo_value(context,
+                                     &(qmckl_mat(mo_value_at_nucl_no_s,0,0)), 
+                                     ctx->mo_basis.mo_num * ctx->nucleus.num);
+    if (rc != QMCKL_SUCCESS) return rc;
+  }
+
+  
+  // Evaluate MO vgl at r_cusp without s components
+  qmckl_tensor mo_vgl_at_r_cusp_no_s;
+  {
+    int64_t sze[3] = { ctx->mo_basis.mo_num, 5, ctx->nucleus.num };
+    mo_vgl_at_r_cusp_no_s = qmckl_tensor_alloc(context, 3, &(sze[0]));
+    
+    rc = qmckl_double_of_matrix(context, ctx->nucleus.coord, coord, ctx->nucleus.num * 3);
+    if (rc != QMCKL_SUCCESS) return rc;
+
+    int64_t i=0;
+    for (int64_t inucl=0 ; inuclnucleus.num ; ++inucl) {
+      coord[i] += r_cusp[inucl];
+      i+=3;
+    }
+      
+    rc = qmckl_set_point(context, 'T', ctx->nucleus.num, coord, (ctx->nucleus.num * 3));
+    if (rc != QMCKL_SUCCESS) {
+      return qmckl_failwith( context,
+                             QMCKL_FAILURE,
+                             "qmckl_set_mo_basis_r_cusp",
+                             "Unable to set coordinates at r_cusp");
+    }
+    
+    rc = qmckl_get_mo_basis_mo_vgl(context,
+                                   &(qmckl_ten3(mo_vgl_at_r_cusp_no_s,0,0,0)), 
+                                   ctx->mo_basis.mo_num * 5 * ctx->nucleus.num);
+    if (rc != QMCKL_SUCCESS) return rc;
+    
+  }
+
+  // Compute parameters
+  {
+     for (int64_t inucl=0 ; inucl < ctx->nucleus.num ; ++inucl) {
+       const double Z = qmckl_vec(ctx->nucleus.charge,inucl);
+       if (Z < 0.1) continue;  // Avoid dummy atoms
+       const double R = r_cusp[inucl];
+       for (int64_t i=0 ; imo_basis.mo_num ; ++i) {
+         const double phi        = qmckl_ten3(mo_vgl_at_r_cusp,i,0,inucl) - qmckl_ten3(mo_vgl_at_r_cusp_no_s,i,0,inucl);
+         const double grad_phi   = qmckl_ten3(mo_vgl_at_r_cusp,i,1,inucl) - qmckl_ten3(mo_vgl_at_r_cusp_no_s,i,1,inucl);
+         const double lap_phi    = qmckl_ten3(mo_vgl_at_r_cusp,i,4,inucl) - qmckl_ten3(mo_vgl_at_r_cusp_no_s,i,4,inucl);
+         const double eta        = qmckl_mat(mo_value_at_nucl_no_s,i,inucl);
+
+         qmckl_ten3(ctx->mo_basis.r_cusp_param,0,i,inucl) =
+           -(R*(2.*eta*Z-6.*grad_phi)+lap_phi*R*R+6.*phi)/(2.*R*Z-6.);
+
+         qmckl_ten3(ctx->mo_basis.r_cusp_param,1,i,inucl) =
+           (lap_phi*R*R*Z-6.*grad_phi*R*Z+6.*phi*Z+6.*eta*Z)/(2.*R*Z-6.);
+
+         qmckl_ten3(ctx->mo_basis.r_cusp_param,2,i,inucl) =
+           -(R*(-5.*grad_phi*Z-1.5*lap_phi)+lap_phi*R*R*Z+3.*phi*Z+3.*eta*Z+6.*grad_phi)/(R*R*Z-3.*R);
+
+         qmckl_ten3(ctx->mo_basis.r_cusp_param,3,i,inucl) =
+           (R*(-2.*grad_phi*Z-lap_phi)+0.5*lap_phi*R*R*Z+phi*Z+eta*Z+3.*grad_phi)/(R*R*R*Z-3.*R*R);
+
+printf("%ld %ld %f %f %f %f\n",i, inucl, 
+                qmckl_ten3(ctx->mo_basis.r_cusp_param,0,i,inucl),
+                qmckl_ten3(ctx->mo_basis.r_cusp_param,1,i,inucl),
+                qmckl_ten3(ctx->mo_basis.r_cusp_param,2,i,inucl),
+                qmckl_ten3(ctx->mo_basis.r_cusp_param,3,i,inucl));
+       }
+     }
+  }
+  
+  free(coord);
+  qmckl_matrix_free(context, &mo_value_at_nucl_no_s);
+  qmckl_tensor_free(context, &mo_vgl_at_r_cusp_no_s);
+  qmckl_tensor_free(context, &mo_vgl_at_r_cusp);
+  
+  // Restore old points
+  if (old_point_num > 0) {
+    rc = qmckl_set_point(context, 'T', old_point_num, old_coord, (old_point_num * 3));
+    if (rc != QMCKL_SUCCESS) {
+      return qmckl_failwith( context,
+                             QMCKL_FAILURE,
+                             "qmckl_set_mo_basis_r_cusp",
+                             "Unable to set old coordinates");
+    }
+    rc = qmckl_free(context, old_coord);
+    if (rc != QMCKL_SUCCESS) return rc;
+    old_coord = NULL;
+  }
+
+  return QMCKL_SUCCESS;
+}
+
+   #+end_src
+
 ** Access functions
 
    #+begin_src c :comments org :tangle (eval h_func) :exports none
@@ -908,10 +1099,11 @@ qmckl_exit_code qmckl_provide_mo_basis_mo_value(qmckl_context context)
       } else {
         rc = qmckl_provide_en_distance(context);
         if (rc != QMCKL_SUCCESS) {
+          return rc;
           return qmckl_failwith( context,
                                  QMCKL_NOT_PROVIDED,
-                                 "qmckl_electron_en_distance",
-                                 NULL);
+                                 "qmckl_provide_mo_basis_mo_value",
+                                 "en_distance");
         }
 
         rc = qmckl_compute_mo_basis_mo_value_cusp(context,
diff --git a/org/qmckl_trexio.org b/org/qmckl_trexio.org
index 0deb13e..e73d424 100644
--- a/org/qmckl_trexio.org
+++ b/org/qmckl_trexio.org
@@ -267,7 +267,7 @@ qmckl_trexio_read_nucleus_X(qmckl_context context, trexio_t* const file)
 
 
    #+begin_src c :tangle (eval c)
-
+  assert ( qmckl_nucleus_provided(context) );
   return QMCKL_SUCCESS;
 }
 #endif