Added asymp_a

org/qmckl_jastrow.org

@@ -18,9 +18,10 @@
   $J_{\text{eN}}$ contains electron-nucleus  terms:
 
   \[
-  J_{\text{eN}}(\mathbf{r},\mathbf{R}) = \sum_{i=1}^{N_\text{elec}} \sum_{\alpha=1}^{N_\text{nucl}}
+  J_{\text{eN}}(\mathbf{r},\mathbf{R}) =
+  \sum_{\alpha=1}^{N_\text{nucl}} \sum_{i=1}^{N_\text{elec}} 
   \frac{a_{1,\alpha}\, g_\alpha(R_{i\alpha})}{1+a_{2,\alpha}\, g_\alpha(R_{i\alpha})} +
-  \sum_{p=2}^{N_\text{ord}^a} a_{p+1,\alpha}\, [g_\alpha(R_{i\alpha})]^p - J_{eN}^\infty
+  \sum_{p=2}^{N_\text{ord}^a} a_{p+1,\alpha}\, [g_\alpha(R_{i\alpha})]^p - J_{eN}^{\infty \alpha}
   \]
 
   $J_{\text{ee}}$ contains electron-electron terms:
@@ -166,6 +167,8 @@ int main() {
     |-------------------------------------+-------------------------------------------------------------------+---------------------------------------------------------------------------------------------------------|
     | ~dim_c_vector~                      | ~int64_t~                                                         | Number of unique C coefficients                                                                         |
     | ~dim_c_vector_date~                 | ~uint64_t~                                                        | Number of unique C coefficients                                                                         |
+    | ~asymp_jasa~                        | ~double[nucl_num]~                                                | Asymptotic component                                                                                    |
+    | ~asymp_jasa_date~                   | ~uint64_t~                                                        | Ladt modification of the asymptotic component                                                           |
     | ~asymp_jasb~                        | ~double[2]~                                                       | Asymptotic component (up- or down-spin)                                                                 |
     | ~asymp_jasb_date~                   | ~uint64_t~                                                        | Ladt modification of the asymptotic component                                                           |
     | ~c_vector_full~                     | ~double[dim_c_vector][nucl_num]~                                  | vector of non-zero coefficients                                                                         |
@@ -202,7 +205,7 @@ nucl_coord = np.array([ [0.000000,  0.000000 ],
    [0.000000,  0.000000 ],
    [0.000000,  2.059801 ] ])
 
-elec_coord =    [[[-0.250655104764153      ,  0.503070975550133      ,  -0.166554344502303],
+elec_coord = np.array(   [[[-0.250655104764153      ,  0.503070975550133      ,  -0.166554344502303],
      [-0.587812193472177      , -0.128751981129274      ,   0.187773606533075],
      [ 1.61335569047166       , -0.615556732874863      ,  -1.43165470979934 ],
      [-4.901239896295210E-003 , -1.120440036458986E-002 ,   1.99761909330422 ],
@@ -211,9 +214,9 @@ elec_coord =    [[[-0.250655104764153      ,  0.503070975550133      ,  -0.16655
      [-0.232271834949124      , -1.059321673434182E-002 ,  -0.504862241464867],
      [ 1.09360863531826       , -2.036103063808752E-003 ,  -2.702796910818986E-002],
      [-0.108090166832043      ,  0.189161729653261      ,   2.15398313919894],
-     [ 0.397978144318712      , -0.254277292595981      ,   2.54553335476344]]];
+     [ 0.397978144318712      , -0.254277292595981      ,   2.54553335476344]]])
 
-ee_distance_rescaled = [
+ee_distance_rescaled = np.array([
 [ 0.000000000000000, 0.000000000000000, 0.000000000000000,
   0.000000000000000, 0.000000000000000, 0.000000000000000,
   0.000000000000000, 0.000000000000000, 0.000000000000000,
@@ -253,7 +256,7 @@ ee_distance_rescaled = [
 [ 0.944400908070716, 0.922589018494961, 0.984615718580670,
   0.514328661540623, 0.692362267147064, 0.931894098453677,
   0.956034127544344, 0.931221472309472, 0.540903688625053,
-  0.000000000000000 ]]
+  0.000000000000000 ]])
 
 en_distance_rescaled = np.transpose(np.array([
 [ 0.443570948411811  ,    0.467602196999105    ,  0.893870160799932 ,
@@ -276,16 +279,16 @@ bord_num     = 5
 cord_num     = 5
 dim_c_vector= 23
 type_nucl_vector = [ 1, 1]
-a_vector = [
+a_vector = np.array([
 [0.000000000000000E+000],
 [0.000000000000000E+000],
 [-0.380512000000000E+000],
 [-0.157996000000000E+000],
 [-3.155800000000000E-002],
-[2.151200000000000E-002]]
+[2.151200000000000E-002]])
 
-b_vector = [ 0.500000000000000E-000,  0.153660000000000E-000,  6.722620000000000E-002,
+b_vector =np.array( [ 0.500000000000000E-000,  0.153660000000000E-000,  6.722620000000000E-002,
-  2.157000000000000E-002,  7.309600000000000E-003,  2.866000000000000E-003]
+  2.157000000000000E-002,  7.309600000000000E-003,  2.866000000000000E-003])
 c_vector = [ 0.571702000000000E-000, -0.514253000000000E-000, -0.513043000000000E-000,
   9.486000000000000E-003, -4.205000000000000E-003,  0.426325800000000E-000,
   8.288150000000000E-002,  5.118600000000000E-003, -2.997800000000000E-003,
@@ -349,6 +352,7 @@ typedef struct qmckl_jastrow_struct{
   int64_t   bord_num;
   int64_t   cord_num;
   int64_t   type_nucl_num;
+  uint64_t  asymp_jasa_date;
   uint64_t  asymp_jasb_date;
   uint64_t  tmp_c_date;
   uint64_t  dtmp_c_date;
@@ -364,6 +368,7 @@ typedef struct qmckl_jastrow_struct{
   double *  a_vector;
   double *  b_vector;
   double *  c_vector;
+  double *  asymp_jasa;
   double *  asymp_jasb;
   double *  factor_ee;
   double *  factor_en;
@@ -1696,6 +1701,22 @@ qmckl_get_jastrow_asymp_jasb(qmckl_context context,
 }
     #+end_src
 
+**** Fortran interface
+
+#+begin_src f90 :tangle (eval fh_func) :comments org
+interface
+   integer(qmckl_exit_code) function qmckl_get_jastrow_asymp_jasb(context, &
+        asymp_jasb, size_max) bind(C)
+     use, intrinsic :: iso_c_binding
+     import
+     implicit none
+     integer (qmckl_context) , intent(in), value :: context
+     integer(c_int64_t), intent(in), value       :: size_max
+     double precision, intent(out)               :: asymp_jasb(size_max)
+   end function qmckl_get_jastrow_asymp_jasb
+end interface
+#+end_src
+
 *** Provide                                                        :noexport:
     #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
 qmckl_exit_code qmckl_provide_jastrow_asymp_jasb(qmckl_context context);
@@ -2957,6 +2978,277 @@ assert(fabs(factor_ee_deriv_e[0][2][0]-0.073267755223968  ) < 1.e-12);
 assert(fabs(factor_ee_deriv_e[0][3][0]-1.5111672803213185 ) < 1.e-12);
      #+end_src
 
+** Asymptotic component for \(J_{eN}\)
+
+  Calculate the asymptotic component ~asymp_jasa~ to be substracted from the final
+  electron-nucleus jastrow factor \(J_{\text{eN}}\). The asymptotic component is calculated
+  via the ~a_vector~ and the electron-nucleus rescale factors ~rescale_factor_en~.
+
+  \[
+  J_{\text{en}}^{\infty \alpha} = \frac{a_1 \kappa_\alpha^{-1}}{1 + a_2  \kappa_\alpha^{-1}}
+  \]
+
+*** Get
+    #+begin_src c :comments org :tangle (eval h_func) :noweb yes
+qmckl_exit_code
+qmckl_get_jastrow_asymp_jasa(qmckl_context context,
+                             double* const asymp_jasa,
+                             const int64_t size_max);
+    #+end_src
+
+    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
+qmckl_exit_code
+qmckl_get_jastrow_asymp_jasa(qmckl_context context,
+                             double* const asymp_jasa,
+                             const int64_t size_max)
+{
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_CONTEXT,
+                           "qmckl_get_jastrow_asymp_jasa",
+                           NULL);
+  }
+
+  qmckl_exit_code rc;
+
+  rc = qmckl_provide_jastrow_asymp_jasa(context);
+  if (rc != QMCKL_SUCCESS) return rc;
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
+  assert (ctx != NULL);
+
+  int64_t sze = ctx->nucleus.num;
+  if (size_max < sze) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_3,
+                           "qmckl_get_jastrow_asymp_jasa",
+                           "Array too small. Expected nucleus.num");
+  }
+  memcpy(asymp_jasa, ctx->jastrow.asymp_jasa, sze * sizeof(double));
+
+  return QMCKL_SUCCESS;
+}
+    #+end_src
+
+**** Fortran interface
+
+#+begin_src f90 :tangle (eval fh_func) :comments org
+interface
+   integer(qmckl_exit_code) function qmckl_get_jastrow_asymp_jasa(context, &
+        asymp_jasa, size_max) bind(C)
+     use, intrinsic :: iso_c_binding
+     import
+     implicit none
+     integer (qmckl_context) , intent(in), value :: context
+     integer(c_int64_t), intent(in), value       :: size_max
+     double precision, intent(out)               :: asymp_jasa(size_max)
+   end function qmckl_get_jastrow_asymp_jasa
+end interface
+#+end_src
+
+*** Provide                                                        :noexport:
+    #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
+qmckl_exit_code qmckl_provide_jastrow_asymp_jasa(qmckl_context context);
+    #+end_src
+
+    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
+qmckl_exit_code qmckl_provide_jastrow_asymp_jasa(qmckl_context context)
+{
+
+  qmckl_exit_code rc;
+
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_CONTEXT,
+                           "qmckl_provide_jastrow_asymp_jasa",
+                           NULL);
+  }
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
+  assert (ctx != NULL);
+
+  if (!ctx->jastrow.provided) {
+    return qmckl_failwith( context,
+                           QMCKL_NOT_PROVIDED,
+                           "qmckl_provide_jastrow_asymp_jasa",
+                           NULL);
+  }
+
+ /* Compute if necessary */
+  if (ctx->date > ctx->jastrow.asymp_jasa_date) {
+
+ /* Allocate array */
+    if (ctx->jastrow.asymp_jasa == NULL) {
+
+      qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
+      mem_info.size = ctx->nucleus.num * sizeof(double);
+      double* asymp_jasa = (double*) qmckl_malloc(context, mem_info);
+
+      if (asymp_jasa == NULL) {
+        return qmckl_failwith( context,
+                               QMCKL_ALLOCATION_FAILED,
+                               "qmckl_asymp_jasa",
+                               NULL);
+      }
+      ctx->jastrow.asymp_jasa = asymp_jasa;
+    }
+
+    rc = qmckl_compute_jastrow_asymp_jasa(context,
+                                  ctx->jastrow.aord_num,
+                                  ctx->jastrow.a_vector,
+                                  ctx->jastrow.rescale_factor_en,
+                                  ctx->nucleus.num,
+                                  ctx->jastrow.asymp_jasa);
+    if (rc != QMCKL_SUCCESS) {
+      return rc;
+    }
+
+    ctx->jastrow.asymp_jasa_date = ctx->date;
+  }
+
+  return QMCKL_SUCCESS;
+}
+    #+end_src
+
+*** Compute
+   :PROPERTIES:
+   :Name:     qmckl_compute_jastrow_asymp_jasa
+   :CRetType: qmckl_exit_code
+   :FRetType: qmckl_exit_code
+   :END:
+
+    #+NAME: qmckl_asymp_jasa_args
+    | Variable            | Type                 | In/Out | Description                |
+    |---------------------+----------------------+--------+----------------------------|
+    | ~context~           | ~qmckl_context~      | in     | Global state               |
+    | ~aord_num~          | ~int64_t~            | in     | Order of the polynomial    |
+    | ~a_vector~          | ~double[aord_num+1]~ | in     | Values of a                |
+    | ~rescale_factor_en~ | ~double~             | in     | Electron nucleus distances |
+    | ~nucl_num~          | ~int64_t~            | in     | Number of nuclei           |
+    | ~asymp_jasa~        | ~double[nucl_num]~   | out    | Asymptotic value           |
+
+    #+begin_src f90 :comments org :tangle (eval f) :noweb yes
+integer function qmckl_compute_jastrow_asymp_jasa_f(context, aord_num, a_vector, &
+     rescale_factor_en, nucl_num, asymp_jasa) &
+     result(info)
+  use qmckl
+  implicit none
+  integer(qmckl_context), intent(in)  :: context
+  integer*8             , intent(in)  :: aord_num
+  double precision      , intent(in)  :: a_vector(aord_num + 1)
+  integer*8             , intent(in)  :: nucl_num
+  double precision      , intent(in)  :: rescale_factor_en(nucl_num)
+  double precision      , intent(out) :: asymp_jasa(nucl_num)
+
+  integer*8 :: i, j, p
+  double precision   :: kappa_inv, x, asym_one
+
+  
+  info = QMCKL_SUCCESS
+
+  if (context == QMCKL_NULL_CONTEXT) then
+     info = QMCKL_INVALID_CONTEXT
+     return
+  endif
+
+  if (aord_num < 0) then
+     info = QMCKL_INVALID_ARG_2
+     return
+  endif
+
+  do i=1,nucl_num
+
+     kappa_inv = 1.0d0 / rescale_factor_en(i)
+     
+     asymp_jasa(i) = a_vector(1) * kappa_inv / (1.0d0 + a_vector(2) * kappa_inv)
+     
+     x = kappa_inv
+     do p = 1, aord_num
+        x = x * kappa_inv
+        asymp_jasa(i) = asymp_jasa(i) + a_vector(p + 1) * x
+     end do
+     
+  end do
+
+end function qmckl_compute_jastrow_asymp_jasa_f
+    #+end_src
+
+#+begin_src c :comments org :tangle (eval c) :noweb yes
+qmckl_exit_code qmckl_compute_jastrow_asymp_jasa (
+	  const qmckl_context context,
+	  const int64_t aord_num,
+	  const double* a_vector,
+	  double* const rescale_factor_en,
+	  const int64_t nucl_num,
+	  double* const asymp_jasa ) {
+
+    if (context == QMCKL_NULL_CONTEXT){
+      return QMCKL_INVALID_CONTEXT;
+    }
+
+    if (aord_num < 0) {
+      return QMCKL_INVALID_ARG_2;
+    }
+
+    for (int i = 0 ; i <= nucl_num; ++i) {
+      const double kappa_inv = 1.0 / rescale_factor_en[i];
+      asymp_jasa[i] =  a_vector[0] * kappa_inv / (1.0 + a_vector[1] * kappa_inv);
+      
+      double x = kappa_inv;
+      for (int p = 1; p < aord_num; ++p){
+        x *= kappa_inv;
+        asymp_jasa[i] = asymp_jasa[i] + a_vector[p + 1] * x;
+      }
+    }
+
+  return QMCKL_SUCCESS;
+}
+#+end_src
+
+#   #+CALL: generate_c_header(table=qmckl_asymp_jasa_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
+
+    #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
+    qmckl_exit_code qmckl_compute_jastrow_asymp_jasa (
+          const qmckl_context context,
+          const int64_t aord_num,
+          const double* a_vector,
+          double* const rescale_factor_en,
+          const int64_t nucl_num,
+          double* const asymp_jasa );
+    #+end_src
+
+*** Test
+    #+name: asymp_jasa
+    #+begin_src python :results output :exports none :noweb yes
+import numpy as np
+
+<<jastrow_data>>
+
+
+asymp_jasa = a_vector[0] * kappa_inv / (1.0 + a_vector[1]*kappa_inv)
+x = kappa_inv
+for p in range(1,aord_num):
+  x = x * kappa_inv
+  asymp_jasa += a_vector[p + 1] * x
+print("asymp_jasa[i]    : ", asymp_jasa)
+
+    #+end_src
+
+    #+RESULTS: asymp_jasa
+    : asymp_jasa[i]    :  [-0.548554]
+
+     #+begin_src c :tangle (eval c_test)
+double asymp_jasa[2];
+rc = qmckl_get_jastrow_asymp_jasa(context, asymp_jasa, nucl_num);
+
+// calculate asymp_jasb
+printf("%e %e\n", asymp_jasa[0], -0.548554);
+printf("%e %e\n", asymp_jasa[1], -0.548554);
+assert(fabs(-0.548554 - asymp_jasa[0]) < 1.e-12);
+assert(fabs(-0.548554 - asymp_jasa[1]) < 1.e-12);
+
+     #+end_src
+
 ** Electron-nucleus component \(f_{en}\)
 
   Calculate the electron-electron jastrow component ~factor_en~ using the ~a_vector~