Addd AO functions

org/qmckl_ao.org

@@ -106,6 +106,9 @@ int main() {
   | ~exponent~          | ~[prim_num]~  | Array of exponents                                                   |
   | ~coefficient~       | ~[prim_num]~  | Array of coefficients                                                |
   | ~prim_factor~       | ~[prim_num]~  | Normalization factors of the primtives                               |
+  | ~ao_num~            |               | Number of AOs                                                        |
+  | ~ao_factor~         | ~[ao_num]~    | Normalization factor of the AO                                       |
+  | ~ao_cartesian~      |               | If true, use polynomials. Otherwise, use spherical harmonics         |
 
   Computed data:
   
@@ -177,18 +180,19 @@ prim_factor = [ 1.0006253235944540e+01, 2.4169531573445120e+00, 7.96109248497664
 
    #+begin_src c :comments org :tangle (eval h_private_type)
 typedef struct qmckl_ao_basis_struct {
-  int32_t   uninitialized;
   int64_t   shell_num;
   int64_t   prim_num;
+  int64_t   ao_num;  
   int64_t * nucleus_index;
   int64_t * nucleus_shell_num;
   int32_t * shell_ang_mom;
   int64_t * shell_prim_num;
   int64_t * shell_prim_index;
   double  * shell_factor;
-  double  * exponent    ;
+  double  * exponent;
-  double  * coefficient ;
+  double  * coefficient;
-  double  * prim_factor ;
+  double  * prim_factor;
+  double  * ao_factor;
 
   int64_t * nucleus_prim_index;
   double  * coefficient_normalized;
@@ -198,7 +202,10 @@ typedef struct qmckl_ao_basis_struct {
   int64_t   primitive_vgl_date;
   double  * shell_vgl;
   int64_t   shell_vgl_date;
+
+  int32_t   uninitialized;
   bool      provided;
+  bool      ao_cartesian;
   char      type;
 } qmckl_ao_basis_struct;
    #+end_src
@@ -224,11 +231,10 @@ qmckl_exit_code qmckl_init_ao_basis(qmckl_context context) {
   qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
   assert (ctx != NULL);
 
-  ctx->ao_basis.uninitialized = (1 << 12) - 1;
+  ctx->ao_basis.uninitialized = (1 << 14) - 1;
 
  /* Default values */
-  /* ctx->ao_basis.
+  ctx->ao_basis.ao_cartesian = true;
-     ,*/
 
   return QMCKL_SUCCESS;
 }
@@ -240,6 +246,7 @@ qmckl_exit_code qmckl_init_ao_basis(qmckl_context context) {
 char      qmckl_get_ao_basis_type             (const qmckl_context context);
 int64_t   qmckl_get_ao_basis_shell_num        (const qmckl_context context);
 int64_t   qmckl_get_ao_basis_prim_num         (const qmckl_context context);
+int64_t   qmckl_get_ao_basis_ao_num           (const qmckl_context context);
 int64_t*  qmckl_get_ao_basis_nucleus_index    (const qmckl_context context);
 int32_t*  qmckl_get_ao_basis_shell_ang_mom    (const qmckl_context context);
 int64_t*  qmckl_get_ao_basis_shell_prim_num   (const qmckl_context context);
@@ -248,6 +255,7 @@ double*   qmckl_get_ao_basis_shell_factor     (const qmckl_context context);
 double*   qmckl_get_ao_basis_exponent         (const qmckl_context context);
 double*   qmckl_get_ao_basis_coefficient      (const qmckl_context context);
 double*   qmckl_get_ao_basis_prim_factor      (const qmckl_context context);
+double*   qmckl_get_ao_basis_ao_factor        (const qmckl_context context);
    #+end_src
 
    When all the data for the AOs have been provided, the following
@@ -494,6 +502,44 @@ double* qmckl_get_ao_basis_prim_factor (const qmckl_context context) {
 }
 
 
+int64_t qmckl_get_ao_basis_ao_num (const qmckl_context context) {
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return (int64_t) 0;
+  }
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  assert (ctx != NULL);
+
+  int32_t mask = 1 << 12;
+
+  if ( (ctx->ao_basis.uninitialized & mask) != 0) {
+    return (int64_t) 0;
+  }
+
+  assert (ctx->ao_basis.ao_num > (int64_t) 0);
+  return ctx->ao_basis.ao_num;
+}
+
+
+double* qmckl_get_ao_basis_ao_factor (const qmckl_context context) {
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return NULL;
+  }
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  assert (ctx != NULL);
+
+  int32_t mask = 1 << 13;
+
+  if ( (ctx->ao_basis.uninitialized & mask) != 0) {
+    return NULL;
+  }
+
+  assert (ctx->ao_basis.ao_factor != NULL);
+  return ctx->ao_basis.ao_factor;
+}
+
+
 bool qmckl_ao_basis_provided(const qmckl_context context) {
 
   if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
@@ -505,6 +551,9 @@ bool qmckl_ao_basis_provided(const qmckl_context context) {
 
   return ctx->ao_basis.provided;
 }
+
+
+
    #+end_src
 
 ** Initialization functions
@@ -516,6 +565,7 @@ bool qmckl_ao_basis_provided(const qmckl_context context) {
 qmckl_exit_code  qmckl_set_ao_basis_type             (qmckl_context context, const char      t);
 qmckl_exit_code  qmckl_set_ao_basis_shell_num        (qmckl_context context, const int64_t   shell_num);
 qmckl_exit_code  qmckl_set_ao_basis_prim_num         (qmckl_context context, const int64_t   prim_num);
+qmckl_exit_code  qmckl_set_ao_basis_ao_num           (qmckl_context context, const int64_t   ao_num);
 qmckl_exit_code  qmckl_set_ao_basis_nucleus_index    (qmckl_context context, const int64_t * nucleus_index);
 qmckl_exit_code  qmckl_set_ao_basis_nucleus_shell_num(qmckl_context context, const int64_t * nucleus_shell_num);
 qmckl_exit_code  qmckl_set_ao_basis_shell_ang_mom    (qmckl_context context, const int32_t * shell_ang_mom);
@@ -525,6 +575,8 @@ qmckl_exit_code  qmckl_set_ao_basis_shell_factor     (qmckl_context context, con
 qmckl_exit_code  qmckl_set_ao_basis_exponent         (qmckl_context context, const double  * exponent);
 qmckl_exit_code  qmckl_set_ao_basis_coefficient      (qmckl_context context, const double  * coefficient);
 qmckl_exit_code  qmckl_set_ao_basis_prim_factor      (qmckl_context context, const double  * prim_factor);
+qmckl_exit_code  qmckl_set_ao_basis_ao_factor        (qmckl_context context, const double  * ao_factor);
+qmckl_exit_code  qmckl_set_ao_basis_cartesian        (qmckl_context context, const bool      cartesian);
    #+end_src
 
    #+NAME:pre2
@@ -945,6 +997,7 @@ qmckl_exit_code  qmckl_set_ao_basis_coefficient(qmckl_context context, const dou
   <<post2>>
 }
 
+
 qmckl_exit_code  qmckl_set_ao_basis_prim_factor(qmckl_context context, const double* prim_factor) {
   <<pre2>>
 
@@ -986,6 +1039,83 @@ qmckl_exit_code  qmckl_set_ao_basis_prim_factor(qmckl_context context, const dou
   <<post2>>
 }
 
+
+qmckl_exit_code  qmckl_set_ao_basis_ao_num(qmckl_context context, const int64_t ao_num) {
+  <<pre2>>
+
+  if (ao_num <= 0) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_2,
+                           "qmckl_set_ao_basis_shell_num",
+                           "ao_num must be positive");
+  }
+
+  int64_t shell_num = qmckl_get_ao_basis_shell_num(context);
+
+  if (ao_num < shell_num) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_2,
+                           "qmckl_set_ao_basis_shell_num",
+                           "ao_num < shell_num");
+  }
+
+  int32_t mask = 1 << 12;
+  ctx->ao_basis.ao_num = ao_num;
+
+  <<post2>>
+}
+
+
+qmckl_exit_code  qmckl_set_ao_basis_ao_factor(qmckl_context context, const double* ao_factor) {
+  <<pre2>>
+
+  int32_t mask = 1 << 13;
+
+  const int64_t ao_num = qmckl_get_ao_basis_ao_num(context);
+  if (ao_num == 0L) {
+    return qmckl_failwith( context,
+                           QMCKL_FAILURE,
+                           "qmckl_set_ao_basis_ao_factor",
+                           "ao_num is not set");
+  }
+
+
+  if (ctx->ao_basis.ao_factor != NULL) {
+    qmckl_exit_code rc = qmckl_free(context, ctx->ao_basis.ao_factor);
+    if (rc != QMCKL_SUCCESS) {
+      return qmckl_failwith( context, rc,
+                             "qmckl_set_ao_basis_ao_factor",
+                             NULL);
+    }
+  }
+
+  qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
+  mem_info.size = ao_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_ao_basis_ao_factor",
+                           NULL);
+  }
+
+  memcpy(new_array, ao_factor, mem_info.size);
+
+  ctx->ao_basis.ao_factor = new_array;
+
+  <<post2>>
+}
+
+
+qmckl_exit_code qmckl_set_ao_basis_cartesian(qmckl_context context, const bool t) {
+  <<pre2>>
+
+  int32_t mask = 1;
+  ctx->ao_basis.ao_cartesian = t;
+
+  <<post2>>
+}
    #+end_src
 
  When the basis set is completely entered, other data structures are
@@ -1273,8 +1403,9 @@ assert(rc == QMCKL_SUCCESS);
 assert(qmckl_nucleus_provided(context));
 
 
-const int64_t    shell_num         = chbrclf_shell_num;
+const int64_t    shell_num         =  chbrclf_shell_num;
-const int64_t    prim_num          = chbrclf_prim_num;
+const int64_t    prim_num          =  chbrclf_prim_num;
+const int64_t    ao_num            =  chbrclf_ao_num;
 const int64_t *  nucleus_index     =  &(chbrclf_basis_nucleus_index[0]);
 const int64_t *  nucleus_shell_num =  &(chbrclf_basis_nucleus_shell_num[0]);
 const int32_t *  shell_ang_mom     =  &(chbrclf_basis_shell_ang_mom[0]);
@@ -1284,6 +1415,7 @@ const double  *  shell_factor      =  &(chbrclf_basis_shell_factor[0]);
 const double  *  exponent          =  &(chbrclf_basis_exponent[0]);
 const double  *  coefficient       =  &(chbrclf_basis_coefficient[0]);
 const double  *  prim_factor       =  &(chbrclf_basis_prim_factor[0]);
+const double  *  ao_factor         =  &(chbrclf_basis_ao_factor[0]);
 
 char    typ = 'G';
 
@@ -1335,6 +1467,13 @@ assert(!qmckl_ao_basis_provided(context));
 
 rc = qmckl_set_ao_basis_prim_factor (context, prim_factor);
 assert(rc == QMCKL_SUCCESS);
+
+rc = qmckl_set_ao_basis_ao_num(context, ao_num);
+assert(rc == QMCKL_SUCCESS);
+
+rc = qmckl_set_ao_basis_ao_factor (context, ao_factor);
+assert(rc == QMCKL_SUCCESS);
+
 assert(qmckl_ao_basis_provided(context));
 
    #+end_src

org/qmckl_tests.org

@@ -38,7 +38,7 @@ Br        -1.218470   -0.187436   -0.028227
 #+END_example
 
  Nuclear coordinates are stored in atomic units in transposed format.
- 
+
 #+begin_src c :tangle ../tests/chbrclf.h
 #define chbrclf_nucl_num ((int64_t) 5)
 
@@ -526,15 +526,16 @@ F   1
 #+begin_src c :tangle ../tests/chbrclf.h
 #define chbrclf_shell_num 72
 #define chbrclf_prim_num 297
+#define chbrclf_ao_num 263
 
 int64_t chbrclf_basis_nucleus_index[chbrclf_nucl_num] = {0, 14, 23, 27, 53};
 
 int64_t chbrclf_basis_nucleus_shell_num[chbrclf_nucl_num] = {14, 9, 14, 16, 19};
 
 int32_t chbrclf_basis_shell_ang_mom[chbrclf_shell_num] =
-  {0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 3, 3, 0, 0, 0, 0, 1, 1, 1, 2, 2, 0, 
+  {0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 3, 3, 0, 0, 0, 0, 1, 1, 1, 2, 2, 0,
    0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 3, 3, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1,
-   2, 2, 2, 3, 3, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3}; 
+   2, 2, 2, 3, 3, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3};
 
 int64_t chbrclf_basis_shell_prim_num[chbrclf_shell_num] =
   {10, 10, 1, 1, 1, 5, 1, 1, 1, 1, 1, 1, 1, 1, 5, 1, 1, 1, 1, 1, 1, 1, 1, 10,
@@ -553,7 +554,38 @@ double chbrclf_basis_shell_factor[chbrclf_shell_num] =
   {1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
    1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
    1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
-   1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.}; 
+   1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.};
+
+double chbrclf_basis_ao_factor[chbrclf_ao_num] =
+  {1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
+   1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
+   1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
+   1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
+   1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
+   1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
+   1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
+   1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
+   1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
+   1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
+   1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
+   1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
+   1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
+   1., 1., 1.};
+
+int64_t chbrclf_basis_ao_shell[chbrclf_ao_num] =
+  {0, 1, 2, 3, 4, 5, 5, 5, 6, 6, 6, 7, 7, 7, 8, 8, 8, 9, 9, 9, 9, 9, 9, 10, 10, 10,
+   10, 10, 10, 11, 11, 11, 11, 11, 11, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 13,
+   13, 13, 13, 13, 13, 13, 13, 13, 13, 14, 15, 16, 17, 18, 18, 18, 19, 19, 19, 20,
+   20, 20, 21, 21, 21, 21, 21, 21, 22, 22, 22, 22, 22, 22, 23, 24, 25, 26, 27, 28,
+   28, 28, 29, 29, 29, 30, 30, 30, 31, 31, 31, 32, 32, 32, 32, 32, 32, 33, 33, 33,
+   33, 33, 33, 34, 34, 34, 34, 34, 34, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 36,
+   36, 36, 36, 36, 36, 36, 36, 36, 36, 37, 38, 39, 40, 41, 42, 43, 43, 43, 44, 44,
+   44, 45, 45, 45, 46, 46, 46, 47, 47, 47, 48, 48, 48, 48, 48, 48, 49, 49, 49, 49,
+   49, 49, 50, 50, 50, 50, 50, 50, 51, 51, 51, 51, 51, 51, 51, 51, 51, 51, 52, 52,
+   52, 52, 52, 52, 52, 52, 52, 52, 53, 54, 55, 56, 57, 58, 59, 60, 60, 60, 61, 61,
+   61, 62, 62, 62, 63, 63, 63, 64, 64, 64, 65, 65, 65, 66, 66, 66, 66, 66, 66, 67,
+   67, 67, 67, 67, 67, 68, 68, 68, 68, 68, 68, 69, 69, 69, 69, 69, 69, 70, 70, 70,
+   70, 70, 70, 70, 70, 70, 70, 71, 71, 71, 71, 71, 71, 71, 71, 71, 71};
 
 double chbrclf_basis_exponent[chbrclf_prim_num] =
   {8.2360000000000000e+03, 1.2350000000000000e+03, 2.8080000000000001e+02,
@@ -1017,7 +1049,7 @@ double chbrclf_elec_coord[chbrclf_walk_num][chbrclf_elec_num][3] = { {
 
 * N2
 
-  This test is mainly for the Jastrow factor and was supplied by 
+  This test is mainly for the Jastrow factor and was supplied by
   Ramon Panades Baruetta. The coordinates and Jastrow coefficients
   have been taken from his fork of IRPJast. The core electrons are
   treated by pseudopotentials thus excluded from the actual calculation.
@@ -1037,7 +1069,7 @@ N2
 #+END_example
 
  Nuclear coordinates are stored in atomic units in transposed format.
- 
+
 #+begin_src c :tangle ../tests/n2.h
 #define n2_nucl_num ((int64_t) 2)
 
@@ -1061,11 +1093,11 @@ double n2_nucl_coord[3][n2_nucl_num] =
 #define n2_walk_num     ((int64_t) 1)
 
 double n2_elec_coord[n2_walk_num][n2_elec_num][3] = { {
-    {-0.250655104764153      ,  0.503070975550133      ,  -0.166554344502303},      
+    {-0.250655104764153      ,  0.503070975550133      ,  -0.166554344502303},
-    {-0.587812193472177      , -0.128751981129274      ,   0.187773606533075},    
+    {-0.587812193472177      , -0.128751981129274      ,   0.187773606533075},
-    { 1.61335569047166       , -0.615556732874863      ,  -1.43165470979934 },   
+    { 1.61335569047166       , -0.615556732874863      ,  -1.43165470979934 },
-    {-4.901239896295210E-003 , -1.120440036458986E-002 ,   1.99761909330422 },   
+    {-4.901239896295210E-003 , -1.120440036458986E-002 ,   1.99761909330422 },
-    { 0.766647499681200      , -0.293515395797937      ,   3.66454589201239 },   
+    { 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},
@@ -1075,25 +1107,25 @@ double n2_elec_coord[n2_walk_num][n2_elec_num][3] = { {
 #+end_src
 
 ** Jastrow related data
- 
+
-  This test is mainly for the Jastrow factor and was supplied by 
+  This test is mainly for the Jastrow factor and was supplied by
   Ramon Panades Baruetta.
 
 #+begin_src c :tangle ../tests/n2.h
 /* Jastrow related */
 
-#define n2_type_nucl_num    ((int64_t) 1) 
+#define n2_type_nucl_num    ((int64_t) 1)
-#define n2_aord_num         ((int64_t) 5) 
+#define n2_aord_num         ((int64_t) 5)
-#define n2_bord_num         ((int64_t) 5) 
+#define n2_bord_num         ((int64_t) 5)
-#define n2_cord_num         ((int64_t) 23) 
+#define n2_cord_num         ((int64_t) 23)
 #define n2_dim_cord_vec     ((int64_t) 23)
 
 int64_t n2_type_nucl_vector[n2_nucl_num] = {
   1,
   1};
 
-double n2_aord_vector[n2_aord_num + 1][n2_type_nucl_num] = {  
+double n2_aord_vector[n2_aord_num + 1][n2_type_nucl_num] = {
-  { 0.      }, 
+  { 0.      },
   { 0.      },
   {-0.380512},
   {-0.157996},
@@ -1105,11 +1137,11 @@ double n2_bord_vector[n2_bord_num + 1] = {
    0.15366   ,
    0.0672262 ,
    0.02157   ,
-   0.0073096 , 
+   0.0073096 ,
    0.002866  };
 
-double n2_cord_vector[n2_cord_num][n2_type_nucl_num] = {  
+double n2_cord_vector[n2_cord_num][n2_type_nucl_num] = {
-  { 5.717020e-01}, 
+  { 5.717020e-01},
   {-5.142530e-01},
   {-5.130430e-01},
   { 9.486000e-03},
@@ -1155,11 +1187,11 @@ double n2_cord_vector_full[n2_dim_cord_vec][n2_nucl_num] = {
    { 2.614000e-03,  2.614000e-03},
    {-1.477000e-03, -1.477000e-03},
    {-1.137000e-03, -1.137000e-03},
-   {-4.010475e-02, -4.010475e-02}, 
+   {-4.010475e-02, -4.010475e-02},
    { 6.106710e-03,  6.106710e-03}};
 
 double n2_lkpm_of_cindex[4][n2_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, 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}};