Added factor_een. #22

org/qmckl_jastrow.org

@@ -3641,7 +3641,7 @@ qmckl_exit_code qmckl_provide_lkpm_combined_index(qmckl_context context)
     if (ctx->jastrow.lkpm_combined_index == NULL) {
 
       qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
-      mem_info.size = 4 * ctx->jastrow.dim_cord_vect * sizeof(double);
+      mem_info.size = 4 * ctx->jastrow.dim_cord_vect * sizeof(int64_t);
       int64_t* lkpm_combined_index = (int64_t*) qmckl_malloc(context, mem_info);
 
       if (lkpm_combined_index == NULL) {
@@ -3904,7 +3904,7 @@ integer function qmckl_compute_lkpm_combined_index_f(context, cord_num, dim_cord
   integer(qmckl_context), intent(in)  :: context
   integer*8             , intent(in)  :: cord_num
   integer*8             , intent(in)  :: dim_cord_vect
-  double precision      , intent(out) :: lkpm_combined_index(dim_cord_vect,4)
+  integer*8             , intent(out) :: lkpm_combined_index(dim_cord_vect, 4)
   double precision                    :: x
   integer*8                           :: i, a, k, l, kk, p, lmax, m
 
@@ -3937,10 +3937,10 @@ integer function qmckl_compute_lkpm_combined_index_f(context, cord_num, dim_cord
         if (iand(p - k - l, 1) .eq. 1) cycle
         m = (p - k - l)/2
         kk = kk + 1
-        lkpm_combined_index(1, kk) = l
+        lkpm_combined_index(kk, 1) = l
-        lkpm_combined_index(2, kk) = k
+        lkpm_combined_index(kk, 2) = k
-        lkpm_combined_index(3, kk) = p
+        lkpm_combined_index(kk, 3) = p
-        lkpm_combined_index(4, kk) = m
+        lkpm_combined_index(kk, 4) = m
       end do
     end do
   end do
@@ -4045,6 +4045,389 @@ assert(qmckl_electron_provided(context));
      #+end_src
 
 
+** Electron-electron-nucleus Jastrow \(f_{een}\)
+
+  Calculate the electron-electron-nuclear three-body jastrow component ~factor_een~ 
+  using the above prepared tables.
+
+  TODO: write equations.
+  
+*** Get
+    #+begin_src c :comments org :tangle (eval h_func) :noweb yes
+qmckl_exit_code qmckl_get_jastrow_factor_een(qmckl_context context, double* const factor_een);
+    #+end_src
+
+    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
+qmckl_exit_code qmckl_get_jastrow_factor_een(qmckl_context context, double* const factor_een)
+{
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return QMCKL_NULL_CONTEXT;
+  }
+
+  qmckl_exit_code rc;
+
+  rc = qmckl_provide_factor_een(context);
+  if (rc != QMCKL_SUCCESS) return rc;
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  assert (ctx != NULL);
+
+  int64_t sze = ctx->electron.walk_num * ctx->electron.num;
+  memcpy(factor_een, ctx->jastrow.factor_een, 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_factor_een(qmckl_context context);
+    #+end_src
+    
+    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
+qmckl_exit_code qmckl_provide_factor_een(qmckl_context context)
+{
+
+  qmckl_exit_code rc;
+
+  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);
+
+ /* Check if en rescaled distance is provided */
+  rc = qmckl_provide_een_rescaled_e(context);
+  if(rc != QMCKL_SUCCESS) return rc;
+
+ /* Check if en rescaled distance derivatives is provided */
+  rc = qmckl_provide_een_rescaled_n(context);
+  if(rc != QMCKL_SUCCESS) return rc;
+
+ /* Check if en rescaled distance derivatives is provided */
+  rc = qmckl_provide_cord_vect_full(context);
+  if(rc != QMCKL_SUCCESS) return rc;
+
+ /* Check if en rescaled distance derivatives is provided */
+  rc = qmckl_provide_lkpm_combined_index(context);
+  if(rc != QMCKL_SUCCESS) return rc;
+
+ /* Compute if necessary */
+  if (ctx->date > ctx->jastrow.factor_een_date) {
+
+ /* Allocate array */
+    if (ctx->jastrow.factor_een == NULL) {
+
+      qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
+      mem_info.size = ctx->electron.walk_num * sizeof(double);
+      double* factor_een = (double*) qmckl_malloc(context, mem_info);
+
+      if (factor_een == NULL) {
+        return qmckl_failwith( context,
+                               QMCKL_ALLOCATION_FAILED,
+                               "qmckl_provide_factor_een",
+                               NULL);
+      }
+      ctx->jastrow.factor_een = factor_een;
+    }
+    
+    qmckl_exit_code rc =
+      qmckl_compute_factor_een(context,
+                              ctx->electron.walk_num,
+                              ctx->electron.num,
+                              ctx->nucleus.num,
+                              ctx->jastrow.cord_num,
+                              ctx->jastrow.dim_cord_vect,
+                              ctx->jastrow.cord_vect_full,
+                              ctx->jastrow.lkpm_combined_index,
+                              ctx->jastrow.een_rescaled_e,
+                              ctx->jastrow.een_rescaled_n,
+                              ctx->jastrow.factor_een);
+    if (rc != QMCKL_SUCCESS) {
+      return rc;
+    }
+
+    ctx->jastrow.factor_en_deriv_e_date = ctx->date;
+  }
+
+  return QMCKL_SUCCESS;
+}
+    #+end_src
+
+*** Compute
+   :PROPERTIES:
+   :Name:     qmckl_compute_factor_een
+   :CRetType: qmckl_exit_code
+   :FRetType: qmckl_exit_code
+   :END:
+
+    #+NAME: qmckl_factor_een_args
+   | qmckl_context | context                                                  | in  | Global state                         |
+   | int64_t       | walk_num                                                 | in  | Number of walkers                    |
+   | int64_t       | elec_num                                                 | in  | Number of electrons                  |
+   | int64_t       | nucl_num                                                 | in  | Number of nucleii                    |
+   | int64_t       | cord_num                                                 | in  | order of polynomials                 |
+   | int64_t       | dim_cord_vect                                            | in  | dimension of full coefficient vector |
+   | double        | cord_vect_full[dim_cord_vect][nucl_num]                  | in  | full coefficient vector              |
+   | int64_t       | lkpm_combined_index[4][dim_cord_vect]                    | in  | combined indices                     |
+   | double        | een_rescaled_e[walk_num][elec_num][elec_num][0:cord_num] | in  | Electron-nucleus rescaled            |
+   | double        | een_rescaled_n[walk_num][elec_num][nucl_num][0:cord_num] | in  | Electron-nucleus rescaled factor     |
+   | double        | factor_een[walk_num]                                     | out | Electron-nucleus jastrow             |
+
+    #+begin_src f90 :comments org :tangle (eval f) :noweb yes
+integer function qmckl_compute_factor_een_f(context, walk_num, elec_num, nucl_num, cord_num, dim_cord_vect,            &
+                                           cord_vect_full, lkpm_combined_index,                                       &
+                                           een_rescaled_e, een_rescaled_n, factor_een) &
+     result(info)
+  use qmckl
+  implicit none
+  integer(qmckl_context), intent(in)  :: context
+  integer*8             , intent(in)  :: walk_num, elec_num, cord_num, nucl_num, dim_cord_vect
+  integer*8             , intent(in)  :: lkpm_combined_index(4,dim_cord_vect)
+  double precision      , intent(in)  :: cord_vect_full(dim_cord_vect, nucl_num)
+  double precision      , intent(in)  :: een_rescaled_e(walk_num, elec_num, elec_num, 0:cord_num)
+  double precision      , intent(in)  :: een_rescaled_n(walk_num, elec_num, nucl_num, 0:cord_num)
+  double precision      , intent(out) :: factor_een(walk_num)
+
+  integer*8 :: i, a, j, l, k, p, m, n, nw
+  double precision :: accu, accu2, cn
+
+  info = QMCKL_SUCCESS
+
+  if (context == QMCKL_NULL_CONTEXT) then
+     info = QMCKL_INVALID_CONTEXT
+     return
+  endif
+
+  if (walk_num <= 0) then
+     info = QMCKL_INVALID_ARG_2
+     return
+  endif
+
+  if (elec_num <= 0) then
+     info = QMCKL_INVALID_ARG_3
+     return
+  endif
+
+  if (nucl_num <= 0) then
+     info = QMCKL_INVALID_ARG_4
+     return
+  endif
+
+  if (cord_num <= 0) then
+     info = QMCKL_INVALID_ARG_5
+     return
+  endif
+
+  factor_een = 0.0d0
+  
+  do nw =1, walk_num
+  do n = 1, dim_cord_vect
+    l = lkpm_combined_index(1, n)  
+    k = lkpm_combined_index(2, n)  
+    p = lkpm_combined_index(3, n)  
+    m = lkpm_combined_index(4, n)  
+
+    do a = 1, nucl_num
+      accu2 = 0.0d0
+      cn = cord_vect_full(n, a)
+      do j = 1, elec_num
+        accu = 0.0d0
+        do i = 1, elec_num
+          accu = accu + een_rescaled_e(nw, i, j, k) *       &
+                        een_rescaled_n(nw, i, a, m)
+        end do
+        accu2 = accu2 + accu * een_rescaled_n(nw, j, a, m + l)
+      end do
+      factor_een(nw) = factor_een(nw) + accu2 + cn
+    end do
+  end do
+  end do
+
+end function qmckl_compute_factor_een_f
+    #+end_src
+
+    #+CALL: generate_c_header(table=qmckl_factor_een_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
+
+    #+RESULTS:
+    #+begin_src c :tangle (eval h_func) :comments org
+    qmckl_exit_code qmckl_compute_factor_een (
+	  const qmckl_context context,
+	  const int64_t walk_num,
+	  const int64_t elec_num,
+	  const int64_t nucl_num,
+	  const int64_t cord_num,
+	  const int64_t dim_cord_vect,
+	  const double* cord_vect_full,
+	  const int64_t* lkpm_combined_index,
+	  const double* een_rescaled_e,
+	  const double* een_rescaled_n,
+	  double* const factor_een ); 
+    #+end_src
+
+
+    #+CALL: generate_c_interface(table=qmckl_factor_een_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_factor_een &
+	(context, &
+		 walk_num, &
+		 elec_num, &
+		 nucl_num, &
+		 cord_num, &
+		 dim_cord_vect, &
+		 cord_vect_full, &
+		 lkpm_combined_index, &
+		 een_rescaled_e, &
+		 een_rescaled_n, &
+		 factor_een) &
+	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 :: walk_num
+      integer (c_int64_t) , intent(in)  , value :: elec_num
+      integer (c_int64_t) , intent(in)  , value :: nucl_num
+      integer (c_int64_t) , intent(in)  , value :: cord_num
+      integer (c_int64_t) , intent(in)  , value :: dim_cord_vect
+      real    (c_double ) , intent(in)          :: cord_vect_full(nucl_num,dim_cord_vect)
+      integer (c_int64_t) , intent(in)          :: lkpm_combined_index(dim_cord_vect,4)
+      real    (c_double ) , intent(in)          :: een_rescaled_e(0:cord_num,elec_num,elec_num,walk_num)
+      real    (c_double ) , intent(in)          :: een_rescaled_n(0:cord_num,nucl_num,elec_num,walk_num)
+      real    (c_double ) , intent(out)         :: factor_een(walk_num)
+
+      integer(c_int32_t), external :: qmckl_compute_factor_een_f
+      info = qmckl_compute_factor_een_f &
+	     (context, &
+		 walk_num, &
+		 elec_num, &
+		 nucl_num, &
+		 cord_num, &
+		 dim_cord_vect, &
+		 cord_vect_full, &
+		 lkpm_combined_index, &
+		 een_rescaled_e, &
+		 een_rescaled_n, &
+		 factor_een)
+
+    end function qmckl_compute_factor_een
+    #+end_src
+
+*** Test
+    #+begin_src python :results output :exports none :noweb yes
+import numpy as np
+
+<<jastrow_data>>
+
+kappa = 1.0
+
+elec_coord = np.array(elec_coord)[0]
+nucl_coord = np.array(nucl_coord)
+elnuc_dist = np.zeros(shape=(elec_num, nucl_num),dtype=float)
+for i in range(elec_num):
+  for j in range(nucl_num):
+    elnuc_dist[i, j] = np.linalg.norm(elec_coord[i] - nucl_coord[:,j])
+
+elnuc_dist_deriv_e = np.zeros(shape=(4, elec_num, nucl_num),dtype=float)
+for a in range(nucl_num):
+  for i in range(elec_num):
+    rij_inv = 1.0 / elnuc_dist[i, a]
+    for ii in range(3):
+      elnuc_dist_deriv_e[ii, i, a] = (elec_coord[i][ii] - nucl_coord[ii][a]) * rij_inv
+    elnuc_dist_deriv_e[3, i, a] = 2.0 * rij_inv
+
+en_distance_rescaled_deriv_e = np.zeros(shape=(4,elec_num,nucl_num),dtype=float)
+for a in range(nucl_num):
+  for i in range(elec_num):
+    f = 1.0 - kappa * en_distance_rescaled[i][a]
+    for ii in range(4):
+      en_distance_rescaled_deriv_e[ii][i][a] = elnuc_dist_deriv_e[ii][i][a]
+    en_distance_rescaled_deriv_e[3][i][a] = en_distance_rescaled_deriv_e[3][i][a] + \
+                              (-kappa * en_distance_rescaled_deriv_e[0][i][a] * en_distance_rescaled_deriv_e[0][i][a]) + \
+                              (-kappa * en_distance_rescaled_deriv_e[1][i][a] * en_distance_rescaled_deriv_e[1][i][a]) + \
+                              (-kappa * en_distance_rescaled_deriv_e[2][i][a] * en_distance_rescaled_deriv_e[2][i][a])
+    for ii in range(4):
+      en_distance_rescaled_deriv_e[ii][i][a] = en_distance_rescaled_deriv_e[ii][i][a] * f
+
+third = 1.0 / 3.0
+factor_en_deriv_e = np.zeros(shape=(4,elec_num),dtype=float)
+dx = np.zeros(shape=(4),dtype=float)
+pow_ser_g = np.zeros(shape=(3),dtype=float)
+for a in range(nucl_num):
+   for i in range(elec_num):
+      x = en_distance_rescaled[i][a]
+      if abs(x) < 1e-18:
+        continue
+      pow_ser_g = np.zeros(shape=(3),dtype=float)
+      den         = 1.0 + aord_vector[1][type_nucl_vector[a]-1] * x
+      invden      = 1.0 / den
+      invden2     = invden * invden
+      invden3     = invden2 * invden
+      xinv        = 1.0 / (x + 1.0E-18)
+
+      for ii in range(4):
+        dx[ii] = en_distance_rescaled_deriv_e[ii][i][a]
+
+      lap1 = 0.0
+      lap2 = 0.0
+      lap3 = 0.0
+      for ii in range(3):
+        x = en_distance_rescaled[i][a]
+        if x < 1e-18:
+          continue
+        for p in range(2,aord_num+1):
+          y = p * aord_vector[(p-1) + 1][type_nucl_vector[a]-1] * x
+          pow_ser_g[ii] = pow_ser_g[ii] + y * dx[ii]
+          lap1 = lap1 + (p - 1) * y * xinv * dx[ii] * dx[ii]
+          lap2 = lap2 + y
+          x = x * en_distance_rescaled[i][a]
+        
+        lap3 = lap3 - 2.0 * aord_vector[1][type_nucl_vector[a]-1] * dx[ii] * dx[ii]
+
+        factor_en_deriv_e[ii][i] = factor_en_deriv_e[ii][i] + aord_vector[0][type_nucl_vector[a]-1]  * \
+                              dx[ii] * invden2 + pow_ser_g[ii]
+
+      ii = 3
+      lap2 = lap2 * dx[ii] * third
+      lap3 = lap3 + den * dx[ii]
+      lap3 = lap3 * (aord_vector[0][type_nucl_vector[a]-1] * invden3)
+      factor_en_deriv_e[ii][i] = factor_en_deriv_e[ii][i] + lap1 + lap2 + lap3
+
+print("factor_en_deriv_e[0][0]:",factor_en_deriv_e[0][0])
+print("factor_en_deriv_e[1][0]:",factor_en_deriv_e[1][0])
+print("factor_en_deriv_e[2][0]:",factor_en_deriv_e[2][0])
+print("factor_en_deriv_e[3][0]:",factor_en_deriv_e[3][0])
+
+
+    #+end_src
+
+    #+RESULTS:
+    : factor_en_deriv_e[0][0]: 0.11609919541763383
+    : factor_en_deriv_e[1][0]: -0.23301394780804574
+    : factor_en_deriv_e[2][0]: 0.17548337641865783
+    : factor_en_deriv_e[3][0]: -0.9667363412285741
+    
+
+     #+begin_src c :tangle (eval c_test)
+/* Check if Jastrow is properly initialized */
+assert(qmckl_jastrow_provided(context));
+
+// calculate factor_en_deriv_e
+double factor_een[walk_num];
+rc = qmckl_get_jastrow_factor_een(context, &(factor_een[0]));
+
+//// check factor_en_deriv_e
+//assert(fabs(factor_en_deriv_e[0][0][0]-0.11609919541763383) < 1.e-12);
+//assert(fabs(factor_en_deriv_e[0][1][0]+0.23301394780804574) < 1.e-12);
+//assert(fabs(factor_en_deriv_e[0][2][0]-0.17548337641865783) < 1.e-12);
+//assert(fabs(factor_en_deriv_e[0][3][0]+0.9667363412285741 ) < 1.e-12);
+
+     #+end_src
+     
 * End of files                                                     :noexport:
 
   #+begin_src c :tangle (eval h_private_type)