Jastrow specific test b2.h header added. #22

org/qmckl_jastrow.org

@@ -27,7 +27,7 @@ these factors along with their derivatives.
 #include "config.h"
 #endif
 
-#include "chbrclf.h"
+#include "b2.h"
 
 int main() {
   qmckl_context context;
@@ -76,10 +76,10 @@ int main() {
   | ~int64_t~  | ~aord_num~                                | in  | The number of a coeffecients                                      |
   | ~int64_t~  | ~bord_num~                                | in  | The number of b coeffecients                                      |
   | ~int64_t~  | ~cord_num~                                | in  | The number of c coeffecients                                      |
-  | ~uint64_t~ | ~type_nuc_num~                            | in  | Number of Nucleii types                                           |
+  | ~uint64_t~ | ~type_nucl_num~                            | in  | Number of Nucleii types                                           |
-  | ~double~   | ~aord_vector[aord_num + 1][type_nuc_num]~ | in  | Order of a polynomial coefficients                                |
+  | ~double~   | ~aord_vector[aord_num + 1][type_nucl_num]~ | in  | Order of a polynomial coefficients                                |
   | ~double~   | ~bord_vector[bord_num + 1]~               | in  | Order of b polynomial coefficients                                |
-  | ~double~   | ~cord_vector[cord_num][type_nuc_num]~     | in  | Order of c polynomial coefficients                                |
+  | ~double~   | ~cord_vector[cord_num][type_nucl_num]~     | in  | Order of c polynomial coefficients                                |
   | ~double~   | ~factor_ee~                               | out | Jastrow factor: electron-electron part                            |
   | ~double~   | ~factor_ee_date~                          | out | Jastrow factor: electron-electron part                            |
   | ~double~   | ~factor_en~                               | out | Jastrow factor: electron-nucleus  part                            |
@@ -96,15 +96,18 @@ int main() {
   | ~int64_t~ | ~dim_cord_vec~                                   | Number of unique C coefficients                 |
   | ~double~  | ~asymp_jasb[2]~                                 | Asymptotic component                            |
   | ~int64_t~ | ~asymp_jasb_date~                               | Asymptotic component                            |
-  | ~double~  | ~coord_vect_full[dim_cord_vec][nuc_num]~         | vector of non-zero coefficients                 |
+  | ~double~  | ~coord_vect_full[dim_cord_vec][nucl_num]~         | vector of non-zero coefficients                 |
   | ~int64_t~ | ~lkpm_of_cindex[4][dim_cord_vec]~                | Transform l,k,p, and m into consecutive indices |
-  | ~double~  | ~tmp_c[elec_num][nuc_num][ncord + 1][ncord]~     | vector of non-zero coefficients                 |
+  | ~double~  | ~tmp_c[elec_num][nucl_num][ncord + 1][ncord]~     | vector of non-zero coefficients                 |
-  | ~double~  | ~dtmp_c[elec_num][4][nuc_num][ncord + 1][ncord]~ | vector of non-zero coefficients                 |
+  | ~double~  | ~dtmp_c[elec_num][4][nucl_num][ncord + 1][ncord]~ | vector of non-zero coefficients                 |
 
   For H2O we have the following data:
 
-   #+BEGIN_EXAMPLE
+   #+NAME: jastrow_data
-type_nuc_num = 1
+   #+BEGIN_SRC python :results output
+import numpy as np
+
+type_nucl_num = 1
 aord_num     = 5
 bord_num     = 5
 cord_num     = 23
@@ -139,14 +142,38 @@ cord_vector_full = [
   2.614000000000000E-003, -1.477000000000000E-003, -1.137000000000000E-003,
  -4.010475000000000E-002,  6.106710000000000E-003 ],
 ]
-lkpm_of_cindex = 
+lkpm_of_cindex = [[1 , 1 , 2 , 0], 
-         [ 1, 1, 2, 0, 0, 0, 2, 1, 1, 2, 3, 0, 2, 1, 3, 0, 0, 1,
+          [0 , 0 , 2 , 1],
-           3, 1, 1, 0, 3, 1, 1, 3, 4, 0, 2, 2, 4, 0, 0, 2, 4, 1,
+          [1 , 2 , 3 , 0],
-           3, 1, 4, 0, 1, 1, 4, 1, 2, 0, 4, 1, 0, 0, 4, 2, 1, 4,
+          [2 , 1 , 3 , 0],
-           5, 0, 2, 3, 5, 0, 0, 3, 5, 1, 3, 2, 5, 0, 1, 2, 5, 1,
+          [0 , 1 , 3 , 1],
-           4, 1, 5, 0, 2, 1, 5, 1, 0, 1, 5, 2, 3, 0, 5, 1, 1, 0,
+          [1 , 0 , 3 , 1],
-           5, 2 ]
+          [1 , 3 , 4 , 0],
-   #+END_EXAMPLE
+          [2 , 2 , 4 , 0],
+          [0 , 2 , 4 , 1],
+          [3 , 1 , 4 , 0],
+          [1 , 1 , 4 , 1],
+          [2 , 0 , 4 , 1],
+          [0 , 0 , 4 , 2],
+          [1 , 4 , 5 , 0],
+          [2 , 3 , 5 , 0],
+          [0 , 3 , 5 , 1],
+          [3 , 2 , 5 , 0],
+          [1 , 2 , 5 , 1],
+          [4 , 1 , 5 , 0],
+          [2 , 1 , 5 , 1],
+          [0 , 1 , 5 , 2],
+          [3 , 0 , 5 , 1],
+          [1 , 0 , 5 , 2]]
+
+print(np.array(lkpm_of_cindex).reshape(4,dim_cord_vec))
+   #+END_SRC
+
+   #+RESULTS: jastrow_data
+   : [[1 1 2 0 0 0 2 1 1 2 3 0 2 1 3 0 0 1 3 1 1 0 3]
+   :  [1 1 3 4 0 2 2 4 0 0 2 4 1 3 1 4 0 1 1 4 1 2 0]
+   :  [4 1 0 0 4 2 1 4 5 0 2 3 5 0 0 3 5 1 3 2 5 0 1]
+   :  [2 5 1 4 1 5 0 2 1 5 1 0 1 5 2 3 0 5 1 1 0 5 2]]
 
 ** Data structure
 
@@ -156,7 +183,7 @@ typedef struct qmckl_jastrow_struct{
   int64_t  aord_num;
   int64_t  bord_num;
   int64_t  cord_num;
-  int64_t  type_nuc_num;
+  int64_t  type_nucl_num;
   int64_t  asymp_jasb_date;
   int64_t  tmp_c_date;
   int64_t  dtmp_c_date;
@@ -218,7 +245,7 @@ qmckl_exit_code qmckl_init_jastrow(qmckl_context context) {
 qmckl_exit_code  qmckl_get_jastrow_aord_num       (qmckl_context context, int64_t* const aord_num);
 qmckl_exit_code  qmckl_get_jastrow_bord_num       (qmckl_context context, int64_t* const bord_num);
 qmckl_exit_code  qmckl_get_jastrow_cord_num       (qmckl_context context, int64_t* const bord_num);
-qmckl_exit_code  qmckl_get_jastrow_type_nuc_num   (qmckl_context context, int64_t* const type_nuc_num);
+qmckl_exit_code  qmckl_get_jastrow_type_nucl_num   (qmckl_context context, int64_t* const type_nucl_num);
 qmckl_exit_code  qmckl_get_jastrow_aord_vector    (qmckl_context context, double * const aord_vector);
 qmckl_exit_code  qmckl_get_jastrow_bord_vector    (qmckl_context context, double * const bord_vector);
 qmckl_exit_code  qmckl_get_jastrow_cord_vector    (qmckl_context context, double * const cord_vector);
@@ -338,17 +365,17 @@ qmckl_exit_code qmckl_get_jastrow_cord_num (const qmckl_context context, int64_t
   return QMCKL_SUCCESS;
 }
 
-qmckl_exit_code qmckl_get_type_nuc_num (const qmckl_context context, int64_t* const type_nuc_num) {
+qmckl_exit_code qmckl_get_type_nucl_num (const qmckl_context context, int64_t* const type_nucl_num) {
 
   if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
     return (char) 0;
   }
   
-  if (type_nuc_num == NULL) {
+  if (type_nucl_num == NULL) {
     return qmckl_failwith( context,
                            QMCKL_INVALID_ARG_2,
-                           "qmckl_get_jastrow_type_nuc_num",
+                           "qmckl_get_jastrow_type_nucl_num",
-                           "type_nuc_num is a null pointer");
+                           "type_nucl_num is a null pointer");
   }
 
   qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
@@ -360,8 +387,8 @@ qmckl_exit_code qmckl_get_type_nuc_num (const qmckl_context context, int64_t* co
     return (char) 0;
   }
 
-  assert (ctx->jastrow.type_nuc_num > 0);
+  assert (ctx->jastrow.type_nucl_num > 0);
-  *type_nuc_num = ctx->jastrow.type_nuc_num;
+  *type_nucl_num = ctx->jastrow.type_nucl_num;
   return QMCKL_SUCCESS;
 }
 
@@ -455,7 +482,7 @@ qmckl_exit_code qmckl_get_cord_vector (const qmckl_context context, double * con
 
    #+begin_src c :comments org :tangle (eval h_func)
 qmckl_exit_code  qmckl_set_jastrow_ord_num       (qmckl_context context, const int64_t aord_num, const int64_t bord_num, const int64_t cord_num);
-qmckl_exit_code  qmckl_set_jastrow_type_nuc_num  (qmckl_context context, const int64_t type_nuc_num);
+qmckl_exit_code  qmckl_set_jastrow_type_nucl_num  (qmckl_context context, const int64_t type_nucl_num);
 qmckl_exit_code  qmckl_set_jastrow_aord_vector   (qmckl_context context, const double * aord_vector);
 qmckl_exit_code  qmckl_set_jastrow_bord_vector   (qmckl_context context, const double * bord_vector);
 qmckl_exit_code  qmckl_set_jastrow_cord_vector   (qmckl_context context, const double * cord_vector);
@@ -517,18 +544,18 @@ qmckl_exit_code qmckl_set_jastrow_ord_num(qmckl_context context, const int64_t a
   <<post2>>
 }
 
-qmckl_exit_code qmckl_set_jastrow_type_nuc_num(qmckl_context context, const int64_t type_nuc_num) {
+qmckl_exit_code qmckl_set_jastrow_type_nucl_num(qmckl_context context, const int64_t type_nucl_num) {
 <<pre2>>
 
-  if (type_nuc_num <= 0) {
+  if (type_nucl_num <= 0) {
     return qmckl_failwith( context,
                            QMCKL_INVALID_ARG_2,
-                           "qmckl_set_jastrow_type_nuc_num",
+                           "qmckl_set_jastrow_type_nucl_num",
-                           "type_nuc_num < 0");
+                           "type_nucl_num < 0");
   }
 
   int32_t mask = 1 << 2;
-  ctx->jastrow.type_nuc_num = type_nuc_num;
+  ctx->jastrow.type_nucl_num = type_nucl_num;
 
   <<post2>>
 }
@@ -539,7 +566,7 @@ qmckl_exit_code qmckl_set_jastrow_aord_vector(qmckl_context context, double cons
   int32_t mask = 1 << 3;
 
   int64_t aord_num;
-  int rc = qmckl_get_jastrow_aord_num(context, &aord_num);
+  qmckl_exit_code rc = qmckl_get_jastrow_aord_num(context, &aord_num);
   if (rc != QMCKL_SUCCESS) return rc;
  
   if (aord_num == 0) {
@@ -589,7 +616,7 @@ qmckl_exit_code qmckl_set_jastrow_bord_vector(qmckl_context context, double cons
   int32_t mask = 1 << 4;
 
   int64_t bord_num;
-  int rc = qmckl_get_jastrow_bord_num(context, &bord_num);
+  qmckl_exit_code rc = qmckl_get_jastrow_bord_num(context, &bord_num);
   if (rc != QMCKL_SUCCESS) return rc;
 
   if (bord_num == 0) {
@@ -639,7 +666,7 @@ qmckl_exit_code qmckl_set_jastrow_cord_vector(qmckl_context context, double cons
   int32_t mask = 1 << 5;
 
   int64_t cord_num;
-  int rc = qmckl_get_jastrow_cord_num(context, &cord_num);
+  qmckl_exit_code rc = qmckl_get_jastrow_cord_num(context, &cord_num);
   if (rc != QMCKL_SUCCESS) return rc;
 
   if (cord_num == 0) {
@@ -743,7 +770,7 @@ qmckl_exit_code qmckl_finalize_jastrow(qmckl_context context) {
   }
 
   /* Check for the nucleus data
-     1. nuc_num
+     1. nucl_num
      2. en_distances_rescaled
   */
   if (!(ctx->nucleus.provided)) {
@@ -764,6 +791,107 @@ qmckl_exit_code qmckl_finalize_jastrow(qmckl_context context) {
 }
    #+end_src
 
+** Test
+   #+begin_src c :tangle (eval c_test)
+/* Reference input data */
+int64_t walk_num      = b2_walk_num;
+int64_t elec_num      = b2_elec_num;
+int64_t elec_up_num   = b2_elec_up_num;
+int64_t elec_dn_num   = b2_elec_dn_num;
+double  rescale_factor_kappa_ee   = 1.0;
+double  rescale_factor_kappa_en   = 1.0;
+double  nucl_rescale_factor_kappa = 1.0;
+double* elec_coord    = &(b2_elec_coord[0][0][0]);
+
+int64_t  nucl_num      = b2_nucl_num;
+double*  charge        = b2_charge;
+double*  nucl_coord    = &(b2_nucl_coord[0][0]);
+
+/* --- */
+
+qmckl_exit_code rc;
+
+assert(!qmckl_electron_provided(context));
+
+int64_t n;
+rc = qmckl_get_electron_num (context, &n);
+assert(rc == QMCKL_NOT_PROVIDED);
+
+rc = qmckl_get_electron_up_num (context, &n);
+assert(rc == QMCKL_NOT_PROVIDED);
+
+rc = qmckl_get_electron_down_num (context, &n);
+assert(rc == QMCKL_NOT_PROVIDED);
+
+
+rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num);
+assert(rc == QMCKL_SUCCESS);
+assert(!qmckl_electron_provided(context));
+
+rc = qmckl_get_electron_up_num (context, &n);
+assert(rc == QMCKL_SUCCESS);
+assert(n == elec_up_num);
+
+rc = qmckl_get_electron_down_num (context, &n);
+assert(rc == QMCKL_SUCCESS);
+assert(n == elec_dn_num);
+
+rc = qmckl_get_electron_num (context, &n);
+assert(rc == QMCKL_SUCCESS);
+assert(n == elec_num);
+
+double k_ee = 0.;
+double k_en = 0.;
+rc = qmckl_get_electron_rescale_factor_ee (context, &k_ee);
+assert(rc == QMCKL_SUCCESS);
+assert(k_ee == 1.0);
+
+rc = qmckl_get_electron_rescale_factor_en (context, &k_en);
+assert(rc == QMCKL_SUCCESS);
+assert(k_en == 1.0);
+
+rc = qmckl_set_electron_rescale_factor_en(context, rescale_factor_kappa_en);
+assert(rc == QMCKL_SUCCESS);
+
+rc = qmckl_set_electron_rescale_factor_ee(context, rescale_factor_kappa_ee);
+assert(rc == QMCKL_SUCCESS);
+
+rc = qmckl_get_electron_rescale_factor_ee (context, &k_ee);
+assert(rc == QMCKL_SUCCESS);
+assert(k_ee == rescale_factor_kappa_ee);
+
+rc = qmckl_get_electron_rescale_factor_en (context, &k_en);
+assert(rc == QMCKL_SUCCESS);
+assert(k_en == rescale_factor_kappa_en);
+
+
+int64_t w;
+rc = qmckl_get_electron_walk_num (context, &w);
+assert(rc == QMCKL_NOT_PROVIDED);
+
+
+rc = qmckl_set_electron_walk_num (context, walk_num);
+assert(rc == QMCKL_SUCCESS);
+
+rc = qmckl_get_electron_walk_num (context, &w);
+assert(rc == QMCKL_SUCCESS);
+assert(w == walk_num);
+
+assert(qmckl_electron_provided(context));
+
+rc = qmckl_set_electron_coord (context, 'N', elec_coord);
+assert(rc == QMCKL_SUCCESS);
+
+double elec_coord2[walk_num*3*elec_num];
+
+rc = qmckl_get_electron_coord (context, 'N', elec_coord2);
+assert(rc == QMCKL_SUCCESS);
+for (int64_t i=0 ; i<3*elec_num ; ++i) {
+  assert( elec_coord[i] == elec_coord2[i] );
+ }
+
+    #+end_src
+
 * Computation
 
   The computed data is stored in the context so that it can be reused
@@ -963,6 +1091,59 @@ end function qmckl_compute_asymp_jasb_f
     #+end_src
 
 *** Test
+    #+begin_src python :results output :exports none :noweb yes
+import numpy as np
+
+<<jastrow_data>>
+
+kappa     = 1.0
+kappa_inv = 1.0
+
+asym_one = bord_vector[0] * kappa_inv / (1.0 + bord_vector[1]*kappa_inv)
+asymp_jasb = np.array([asym_one, 0.5 * asym_one])
+
+for i in range(2):
+  x = kappa_inv
+  for p in range(1,bord_num):
+    x = x * kappa_inv
+    asymp_jasb[i] += bord_vector[p + 1] * x
+
+print("asym_one         : ", asym_one)
+print("asymp_jasb       : ", asymp_jasb)
+    #+end_src
+
+    #+RESULTS:
+    : asym_one         :  0.43340325572525706
+    : asymp_jasb       :  [0.53237506 0.31567343]
+
+     #+begin_src c :tangle (eval c_test)
+assert(qmckl_electron_provided(context));
+
+
+double ee_distance[walk_num * elec_num * elec_num];
+rc = qmckl_get_electron_ee_distance(context, ee_distance);
+
+// (e1,e2,w)
+// (0,0,0) == 0.
+//assert(ee_distance[0] == 0.);
+//
+//// (1,0,0) == (0,1,0)
+//assert(ee_distance[1] == ee_distance[elec_num]);
+//
+//// value of (1,0,0)
+//assert(fabs(ee_distance[1]-7.152322512964209) < 1.e-12);
+//
+//// (0,0,1) == 0.
+//assert(ee_distance[elec_num*elec_num] == 0.);
+//
+//// (1,0,1) == (0,1,1)
+//assert(ee_distance[elec_num*elec_num+1] == ee_distance[elec_num*elec_num+elec_num]);
+//
+//// value of (1,0,1)
+//assert(fabs(ee_distance[elec_num*elec_num+1]-6.5517646321055665) < 1.e-12);
+
+     #+end_src
+
 
 * End of files                                                     :noexport:
 
@@ -972,7 +1153,7 @@ end function qmckl_compute_asymp_jasb_f
 
 *** Test
   #+begin_src c :tangle (eval c_test)
-    qmckl_exit_code rc = qmckl_context_destroy(context);
+    rc = qmckl_context_destroy(context);
     assert (rc == QMCKL_SUCCESS);
 
     return 0;

tests/b2.h

@@ -0,0 +1,106 @@
+#define b2_nucl_num ((int64_t) 2)
+
+double b2_charge[b2_nucl_num] = { 5., 5.};
+
+double b2_nucl_coord[3][b2_nucl_num] =
+{ {0.000000,  0.000000 },
+  {0.000000,  0.000000 },
+  {2.059801,  2.059801 } };
+
+#define b2_elec_up_num  ((int64_t) 5)
+#define b2_elec_dn_num  ((int64_t) 5)
+#define b2_elec_num     ((int64_t) 10)
+#define b2_walk_num     ((int64_t) 1)
+
+double b2_elec_coord[b2_walk_num][b2_elec_num][3] = { {
+    {-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 },   
+    { 0.766647499681200      , -0.293515395797937      ,   3.66454589201239 },   
+    {-0.127732483187947      , -0.138975497694196      ,  -8.669850480215846E-002},
+    {-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}}};
+
+/* Jastrow related */
+
+#define b2_type_nucl_num    ((int64_t) 1) 
+#define b2_aord_num         ((int64_t) 5) 
+#define b2_bord_num         ((int64_t) 5) 
+#define b2_cord_num         ((int64_t) 23) 
+#define b2_dim_cord_vec     ((int64_t) 23)
+
+double b2_aord_vector[b2_aord_num + 1][b2_type_nucl_num] = {  
+  { 0.      }, 
+  { 0.      },
+  {-0.380512},
+  {-0.157996},
+  {-0.031558},
+  { 0.021512}};
+
+double b2_bord_vector[b2_bord_num + 1] = {
+   0.5       ,
+   0.15366   ,
+   0.0672262 ,
+   0.02157   ,
+   0.0073096 , 
+   0.002866  };
+
+double b2_cord_vector[b2_cord_num][b2_type_nucl_num] = {  
+  { 5.717020e-01}, 
+  {-5.142530e-01},
+  {-5.130430e-01},
+  { 9.486000e-03},
+  {-4.205000e-03},
+  { 4.263258e-01},
+  { 8.288150e-02},
+  { 5.118600e-03},
+  {-2.997800e-03},
+  {-5.270400e-03},
+  {-7.500000e-05},
+  {-8.301650e-02},
+  { 1.454340e-02},
+  { 5.143510e-02},
+  { 9.250000e-04},
+  {-4.099100e-03},
+  { 4.327600e-03},
+  {-1.654470e-03},
+  { 2.614000e-03},
+  {-1.477000e-03},
+  {-1.137000e-03},
+  {-4.010475e-02},
+  { 6.106710e-03}};
+
+double b2_cord_vector_full[b2_dim_cord_vec][b2_nucl_num] = {
+   { 5.717020e-01,  5.717020e-01},
+   {-5.142530e-01, -5.142530e-01},
+   {-5.130430e-01, -5.130430e-01},
+   { 9.486000e-03,  9.486000e-03},
+   {-4.205000e-03, -4.205000e-03},
+   { 4.263258e-01,  4.263258e-01},
+   { 8.288150e-02,  8.288150e-02},
+   { 5.118600e-03,  5.118600e-03},
+   {-2.997800e-03, -2.997800e-03},
+   {-5.270400e-03, -5.270400e-03},
+   {-7.500000e-05, -7.500000e-05},
+   {-8.301650e-02, -8.301650e-02},
+   { 1.454340e-02,  1.454340e-02},
+   { 5.143510e-02,  5.143510e-02},
+   { 9.250000e-04,  9.250000e-04},
+   {-4.099100e-03, -4.099100e-03},
+   { 4.327600e-03,  4.327600e-03},
+   {-1.654470e-03, -1.654470e-03},
+   { 2.614000e-03,  2.614000e-03},
+   {-1.477000e-03, -1.477000e-03},
+   {-1.137000e-03, -1.137000e-03},
+   {-4.010475e-02, -4.010475e-02}, 
+   { 6.106710e-03,  6.106710e-03}};
+
+double b2_lkpm_of_cindex[4][b2_dim_cord_vec] = {
+  {1, 1, 2, 0, 0, 0, 2, 1, 1, 2, 3, 0, 2, 1, 3, 0, 0, 1, 3, 1, 1, 0, 3}, 
+  {1, 1, 3, 4, 0, 2, 2, 4, 0, 0, 2, 4, 1, 3, 1, 4, 0, 1, 1, 4, 1, 2, 0},
+  {4, 1, 0, 0, 4, 2, 1, 4, 5, 0, 2, 3, 5, 0, 0, 3, 5, 1, 3, 2, 5, 0, 1},
+  {2, 5, 1, 4, 1, 5, 0, 2, 1, 5, 1, 0, 1, 5, 2, 3, 0, 5, 1, 1, 0, 5, 2}};
+