Function to calculate rescaled distances for elec-nucl. #15

org/qmckl_distance.org

@@ -807,6 +807,7 @@ assert(test_qmckl_dist(context) == QMCKL_SUCCESS);
    | int64_t       | ldb       | in  | Leading dimension of array ~B~                |
    | double        | C[n][ldc] | out | Array containing the $m \times n$ matrix $C$  |
    | int64_t       | ldc       | in  | Leading dimension of array ~C~                |
+   | double        | kappa     | in  | Factor for calculating rescaled distances     |
 
 *** Requirements
 
@@ -839,13 +840,14 @@ assert(test_qmckl_dist(context) == QMCKL_SUCCESS);
           const double* B,
           const int64_t ldb,
           double* const C,
-          const int64_t ldc );
+          const int64_t ldc,
+          const double  kappa);
     #+end_src
 
 *** Source
     #+begin_src f90 :tangle (eval f)
 integer function qmckl_distance_rescaled_f(context, transa, transb, m, n, &
-     A, LDA, B, LDB, C, LDC) &
+     A, LDA, B, LDB, C, LDC, kappa) &
      result(info)
   use qmckl
   implicit none
@@ -858,6 +860,7 @@ integer function qmckl_distance_rescaled_f(context, transa, transb, m, n, &
   real*8     , intent(in)  :: B(ldb,*)
   integer*8  , intent(in)  :: ldc
   real*8     , intent(out) :: C(ldc,*)
+  real*8     , intent(in)  :: kappa
 
   integer*8 :: i,j
   real*8    :: x, y, z
@@ -988,7 +991,7 @@ end function qmckl_distance_rescaled_f
    #+RESULTS:
    #+begin_src f90 :tangle (eval f) :comments org :exports none
    integer(c_int32_t) function qmckl_distance_rescaled &
-       (context, transa, transb, m, n, A, lda, B, ldb, C, ldc) &
+       (context, transa, transb, m, n, A, lda, B, ldb, C, ldc, kappa) &
        bind(C) result(info)
 
      use, intrinsic :: iso_c_binding
@@ -1005,10 +1008,11 @@ end function qmckl_distance_rescaled_f
      integer (c_int64_t) , intent(in)  , value :: ldb
      real    (c_double ) , intent(out)         :: C(ldc,n)
      integer (c_int64_t) , intent(in)  , value :: ldc
+     real    (c_double ) , intent(in)  , value :: kappa
 
      integer(c_int32_t), external :: qmckl_distance_rescaled_f
      info = qmckl_distance_rescaled_f &
-            (context, transa, transb, m, n, A, lda, B, ldb, C, ldc)
+            (context, transa, transb, m, n, A, lda, B, ldb, C, ldc, kappa)
 
    end function qmckl_distance_rescaled
    #+end_src
@@ -1019,7 +1023,7 @@ end function qmckl_distance_rescaled_f
    #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
    interface
      integer(c_int32_t) function qmckl_distance_rescaled &
-         (context, transa, transb, m, n, A, lda, B, ldb, C, ldc) &
+         (context, transa, transb, m, n, A, lda, B, ldb, C, ldc, kappa) &
          bind(C)
        use, intrinsic :: iso_c_binding
        import
@@ -1036,6 +1040,7 @@ end function qmckl_distance_rescaled_f
        integer (c_int64_t) , intent(in)  , value :: ldb
        real    (c_double ) , intent(out)         :: C(ldc,n)
        integer (c_int64_t) , intent(in)  , value :: ldc
+       real    (c_double ) , intent(in)  , value :: kappa
 
      end function qmckl_distance_rescaled
    end interface

org/qmckl_electron.org

@@ -1276,6 +1276,285 @@ assert(fabs(en_distance[1][0][1] - 3.1804527583077356) < 1.e-12);
 
      #+end_src
 
+** Electron-nucleus rescaled distances
+
+*** Get
+
+    #+begin_src c :comments org :tangle (eval h_func) :noweb yes
+qmckl_exit_code qmckl_get_electron_en_distance_rescaled(qmckl_context context, double* distance_rescaled);
+    #+end_src
+
+    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
+qmckl_exit_code qmckl_get_electron_en_distance_rescaled(qmckl_context context, double* distance_rescaled)
+{
+
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return QMCKL_NULL_CONTEXT;
+  }
+
+  qmckl_exit_code rc;
+
+  rc = qmckl_provide_en_distance_rescaled(context);
+  if (rc != QMCKL_SUCCESS) return rc;
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  assert (ctx != NULL);
+
+  size_t sze = ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num;
+  memcpy(distance_rescaled, ctx->electron.en_distance_rescaled, sze * sizeof(double));
+
+  return QMCKL_SUCCESS;
+}
+    #+end_src
+
+*** Provide                                                        :noexport:
+
+    #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
+qmckl_exit_code qmckl_provide_en_distance_rescaled(qmckl_context context);
+    #+end_src
+
+    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
+qmckl_exit_code qmckl_provide_en_distance_rescaled(qmckl_context context)
+{
+
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return QMCKL_NULL_CONTEXT;
+  }
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  assert (ctx != NULL);
+
+  if (!(ctx->nucleus.provided)) {
+    return QMCKL_NOT_PROVIDED;
+  }
+
+ /* Compute if necessary */
+  if (ctx->electron.coord_new_date > ctx->electron.en_distance_rescaled_date) {
+
+ /* Allocate array */
+    if (ctx->electron.en_distance_rescaled == NULL) {
+
+      qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
+      mem_info.size = ctx->electron.num * ctx->nucleus.num *
+        ctx->electron.walk_num * sizeof(double);
+      double* en_distance_rescaled = (double*) qmckl_malloc(context, mem_info);
+
+      if (en_distance_rescaled == NULL) {
+        return qmckl_failwith( context,
+                               QMCKL_ALLOCATION_FAILED,
+                               "qmckl_en_distance_rescaled",
+                               NULL);
+      }
+      ctx->electron.en_distance_rescaled = en_distance_rescaled;
+    }
+
+    qmckl_exit_code rc =
+      qmckl_compute_en_distance_rescaled(context,
+                                ctx->electron.num,
+                                ctx->nucleus.num,
+                                ctx->electron.kappa_en,
+                                ctx->electron.walk_num,
+                                ctx->electron.coord_new,
+                                ctx->nucleus.coord,
+                                ctx->electron.en_distance_rescaled);
+    if (rc != QMCKL_SUCCESS) {
+      return rc;
+    }
+
+    ctx->electron.en_distance_rescaled_date = ctx->date;
+  }
+
+  return QMCKL_SUCCESS;
+}
+    #+end_src
+
+*** Compute
+   :PROPERTIES:
+   :Name:     qmckl_compute_en_distance_rescaled
+   :CRetType: qmckl_exit_code
+   :FRetType: qmckl_exit_code
+   :END:
+
+    #+NAME: qmckl_en_distance_rescaled_args
+   | qmckl_context | context                                                 | in  | Global state                      |
+   | int64_t       | elec_num                                                | in  | Number of electrons               |
+   | int64_t       | nucl_num                                                | in  | Number of nuclei                  |
+   | double        | kappa_en                                                | in  | The factor for rescaled distances |
+   | int64_t       | walk_num                                                | in  | Number of walkers                 |
+   | double        | elec_coord[walk_num][3][elec_num]                       | in  | Electron coordinates              |
+   | double        | nucl_coord[3][elec_num]                                 | in  | Nuclear coordinates               |
+   | double        | en_distance_rescaled_date[walk_num][nucl_num][elec_num] | out | Electron-nucleus distances        |
+
+    #+begin_src f90 :comments org :tangle (eval f) :noweb yes
+integer function qmckl_compute_en_distance_rescaled_f(context, elec_num, nucl_num, kappa_en, walk_num, elec_coord, &
+     nucl_coord, en_distance_rescaled) &
+     result(info)
+  use qmckl
+  implicit none
+  integer(qmckl_context), intent(in)  :: context
+  integer*8             , intent(in)  :: elec_num
+  integer*8             , intent(in)  :: nucl_num
+  double precision      , intent(in)  :: kappa_en
+  integer*8             , intent(in)  :: walk_num
+  double precision      , intent(in)  :: elec_coord(elec_num,3,walk_num)
+  double precision      , intent(in)  :: nucl_coord(nucl_num,3)
+  double precision      , intent(out) :: en_distance_rescaled(elec_num,nucl_num,walk_num)
+
+  integer*8 :: k
+
+  info = QMCKL_SUCCESS
+
+  if (context == QMCKL_NULL_CONTEXT) then
+     info = QMCKL_INVALID_CONTEXT
+     return
+  endif
+
+  if (elec_num <= 0) then
+     info = QMCKL_INVALID_ARG_2
+     return
+  endif
+
+  if (nucl_num <= 0) then
+     info = QMCKL_INVALID_ARG_3
+     return
+  endif
+
+  ! TODO: comparison with 0
+  !if (kappa_en <= 0) then
+  !   info = QMCKL_INVALID_ARG_4
+  !   return
+  !endif
+
+  if (walk_num <= 0) then
+     info = QMCKL_INVALID_ARG_5
+     return
+  endif
+
+  do k=1,walk_num
+     info = qmckl_distance_rescaled(context, 'T', 'T', elec_num, nucl_num, &
+          elec_coord(1,1,k), elec_num, &
+          nucl_coord, nucl_num, &
+          en_distance_rescaled(1,1,k), elec_num, kappa_en)
+     if (info /= QMCKL_SUCCESS) then
+        exit
+     endif
+  end do
+
+end function qmckl_compute_en_distance_rescaled_f
+    #+end_src
+
+    #+begin_src c :tangle (eval h_private_func) :comments org :exports none
+qmckl_exit_code qmckl_compute_en_distance_rescaled (
+          const qmckl_context context,
+          const int64_t elec_num,
+          const int64_t nucl_num,
+          const double  kappa_en,
+          const int64_t walk_num,
+          const double* elec_coord,
+          const double* nucl_coord,
+          double* const en_distance_rescaled );
+    #+end_src
+
+    #+CALL: generate_c_interface(table=qmckl_en_distance_rescaled_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
+
+    #+RESULTS:
+    #+begin_src f90 :tangle (eval f) :comments org :exports none
+    integer(c_int32_t) function qmckl_compute_en_distance_rescaled &
+        (context, elec_num, nucl_num, kappa_en, walk_num, elec_coord, nucl_coord, en_distance_rescaled) &
+        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 :: nucl_num
+      real    (c_double ) , intent(in)  , value :: kappa_en
+      integer (c_int64_t) , intent(in)  , value :: walk_num
+      real    (c_double ) , intent(in)          :: elec_coord(elec_num,3,walk_num)
+      real    (c_double ) , intent(in)          :: nucl_coord(elec_num,3)
+      real    (c_double ) , intent(out)         :: en_distance_rescaled(elec_num,nucl_num,walk_num)
+
+      integer(c_int32_t), external :: qmckl_compute_en_distance_rescaled_f
+      info = qmckl_compute_en_distance_rescaled_f &
+             (context, elec_num, nucl_num, kappa_en, walk_num, elec_coord, nucl_coord, en_distance_rescaled)
+
+    end function qmckl_compute_en_distance_rescaled
+    #+end_src
+
+*** Test
+
+    #+begin_src python :results output :exports none
+import numpy as np
+
+elec_1_w1 = np.array( [ -2.26995253563, -5.15737533569, -2.22940072417 ])
+elec_2_w1 = np.array( [  3.51983380318, -1.08717381954, -1.19617708027 ])
+elec_1_w2 = np.array( [ -2.34410619736, -3.20016115904, -1.53496759012 ])
+elec_2_w2 = np.array( [  3.17996025085, -1.40260577202,  1.49473607540 ])
+nucl_1    = np.array( [ 1.096243353458458e+00, 8.907054016973815e-01, 7.777092280258892e-01 ] )
+nucl_2    = np.array( [ 1.168459237342663e+00, 1.125660720053393e+00, 2.833370314829343e+00 ] )
+
+print ( "[0][0][0] : ", np.linalg.norm(elec_1_w1-nucl_1) )
+print ( "[0][1][0] : ", np.linalg.norm(elec_1_w1-nucl_2) )
+print ( "[0][0][1] : ", np.linalg.norm(elec_2_w1-nucl_1) )
+print ( "[1][0][0] : ", np.linalg.norm(elec_1_w2-nucl_1) )
+print ( "[1][1][0] : ", np.linalg.norm(elec_1_w2-nucl_2) )
+print ( "[1][0][1] : ", np.linalg.norm(elec_2_w2-nucl_1) )
+    #+end_src
+
+    #+RESULTS:
+    : [0][0][0] :  7.546738741619978
+    : [0][1][0] :  8.77102435246984
+    : [0][0][1] :  3.698922010513608
+    : [1][0][0] :  5.824059436060509
+    : [1][1][0] :  7.080482110317645
+    : [1][0][1] :  3.1804527583077356
+
+     #+begin_src c :tangle (eval c_test)
+
+assert(qmckl_electron_provided(context));
+
+rc = qmckl_set_nucleus_num (context, nucl_num);
+assert(rc == QMCKL_SUCCESS);
+
+rc = qmckl_set_nucleus_kappa (context, nucl_kappa);
+assert(rc == QMCKL_SUCCESS);
+
+rc = qmckl_set_nucleus_charge (context, charge);
+assert (rc == QMCKL_SUCCESS);
+
+rc = qmckl_set_nucleus_coord (context, 'T', nucl_coord);
+assert (rc == QMCKL_SUCCESS);
+
+assert(qmckl_nucleus_provided(context));
+
+double en_distance_rescaled[walk_num][nucl_num][elec_num];
+
+rc = qmckl_get_electron_en_distance_rescaled(context, &(en_distance[0][0][0]));
+assert (rc == QMCKL_SUCCESS);
+
+// TODO: check exact values
+//// (e,n,w) in Fortran notation
+//// (1,1,1)
+//assert(fabs(en_distance[0][0][0] - 7.546738741619978) < 1.e-12);
+//
+//// (1,2,1)
+//assert(fabs(en_distance[0][1][0] - 8.77102435246984) < 1.e-12);
+//
+//// (2,1,1)
+//assert(fabs(en_distance[0][0][1] - 3.698922010513608) < 1.e-12);
+//
+//// (1,1,2)
+//assert(fabs(en_distance[1][0][0] - 5.824059436060509) < 1.e-12);
+//
+//// (1,2,2)
+//assert(fabs(en_distance[1][1][0] - 7.080482110317645) < 1.e-12);
+//
+//// (2,1,2)
+//assert(fabs(en_distance[1][0][1] - 3.1804527583077356) < 1.e-12);
+
+     #+end_src
+
 * End of files                                                     :noexport:
 
   #+begin_src c :tangle (eval h_private_type)

org/qmckl_nucleus.org

@@ -780,6 +780,7 @@ qmckl_exit_code qmckl_provide_nn_distance_rescaled(qmckl_context context)
   qmckl_exit_code rc =
     qmckl_compute_nn_distance_rescaled(context,
                               ctx->nucleus.num,
+                              ctx->nucleus.kappa,
                               ctx->nucleus.coord,
                               ctx->nucleus.nn_distance_rescaled);
   if (rc != QMCKL_SUCCESS) {
@@ -801,12 +802,13 @@ qmckl_exit_code qmckl_provide_nn_distance_rescaled(qmckl_context context)
    | double        | nn_distance_rescaled[nucl_num][nucl_num] | out | Nucleus-nucleus rescaled distances (au) |
 
     #+begin_src f90 :comments org :tangle (eval f) :noweb yes
-integer function qmckl_compute_nn_distance_rescaled_f(context, nucl_num, coord, nn_distance_rescaled) &
+integer function qmckl_compute_nn_distance_rescaled_f(context, nucl_num, kappa, coord, nn_distance_rescaled) &
      result(info)
   use qmckl
   implicit none
   integer(qmckl_context), intent(in)  :: context
   integer*8             , intent(in)  :: nucl_num
+  double precision      , intent(in)  :: kappa
   double precision      , intent(in)  :: coord(nucl_num,3)
   double precision      , intent(out) :: nn_distance_rescaled(nucl_num,nucl_num)
 
@@ -827,7 +829,7 @@ integer function qmckl_compute_nn_distance_rescaled_f(context, nucl_num, coord,
   info = qmckl_distance_rescaled(context, 'T', 'T', nucl_num, nucl_num, &
           coord, nucl_num, &
           coord, nucl_num, &
-          nn_distance_rescaled, nucl_num)
+          nn_distance_rescaled, nucl_num, kappa)
 
 end function qmckl_compute_nn_distance_rescaled_f
     #+end_src
@@ -836,6 +838,7 @@ end function qmckl_compute_nn_distance_rescaled_f
 qmckl_exit_code qmckl_compute_nn_distance_rescaled (
           const qmckl_context context,
           const int64_t nucl_num,
+          const double  kappa,
           const double* coord,
           double* const nn_distance_rescaled );
     #+end_src
@@ -846,7 +849,7 @@ qmckl_exit_code qmckl_compute_nn_distance_rescaled (
     #+RESULTS:
     #+begin_src f90 :tangle (eval f) :comments org :exports none
     integer(c_int32_t) function qmckl_compute_nn_distance_rescaled &
-        (context, nucl_num, coord, nn_distance_rescaled) &
+        (context, nucl_num, kappa, coord, nn_distance_rescaled) &
         bind(C) result(info)
 
       use, intrinsic :: iso_c_binding
@@ -854,12 +857,13 @@ qmckl_exit_code qmckl_compute_nn_distance_rescaled (
 
       integer (c_int64_t) , intent(in)  , value :: context
       integer (c_int64_t) , intent(in)  , value :: nucl_num
+      real    (c_double ) , intent(in)  , value :: kappa
       real    (c_double ) , intent(in)          :: coord(nucl_num,3)
       real    (c_double ) , intent(out)         :: nn_distance_rescaled(nucl_num,nucl_num)
 
       integer(c_int32_t), external :: qmckl_compute_nn_distance_rescaled_f
       info = qmckl_compute_nn_distance_rescaled_f &
-             (context, nucl_num, coord, nn_distance_rescaled)
+             (context, nucl_num, kappa, coord, nn_distance_rescaled)
 
     end function qmckl_compute_nn_distance_rescaled
     #+end_src