Fixed bug in mo_vgl where ao_vgl was not given as dependency. #41

org/qmckl_context.org

@@ -244,6 +244,12 @@ qmckl_context qmckl_context_create() {
     
     rc = qmckl_init_ao_basis(context);
     assert (rc == QMCKL_SUCCESS);
+    
+    rc = qmckl_init_mo_basis(context);
+    assert (rc == QMCKL_SUCCESS);
+    
+    rc = qmckl_init_determinant(context);
+    assert (rc == QMCKL_SUCCESS);
   }
 
  /* Allocate qmckl_memory_struct */

org/qmckl_determinant.org

@@ -171,6 +171,26 @@ typedef struct qmckl_determinant_struct {
    Some values are initialized by default, and are not concerned by
    this mechanism.
 
+   #+begin_src c :comments org :tangle (eval h_private_func) 
+qmckl_exit_code qmckl_init_determinant(qmckl_context context);
+   #+end_src
+   
+   #+begin_src c :comments org :tangle (eval c)
+qmckl_exit_code qmckl_init_determinant(qmckl_context context) {
+
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return false;
+  }
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  assert (ctx != NULL);
+
+  ctx->det.uninitialized = (1 << 6) - 1;
+
+  return QMCKL_SUCCESS;
+}
+   #+end_src
+   
 ** Access functions
    
    #+begin_src c :comments org :tangle (eval h_private_func) :exports none
@@ -189,6 +209,20 @@ int64_t*   qmckl_get_determinant_mo_index_beta   (const qmckl_context context);
 bool      qmckl_determinant_provided           (const qmckl_context context);
    #+end_src
 
+      #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
+bool qmckl_determinant_provided(const qmckl_context context) {
+
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return false;
+  }
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  assert (ctx != NULL);
+
+  return ctx->det.provided;
+}
+      #+end_src
+
    #+NAME:post
    #+begin_src c :exports none
 if ( (ctx->det.uninitialized & mask) != 0) {
@@ -337,8 +371,8 @@ qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
 ctx->det.uninitialized &= ~mask;
 ctx->det.provided = (ctx->det.uninitialized == 0);
 if (ctx->det.provided) {
-  //qmckl_exit_code rc_ = qmckl_finalize_determinant(context);
+  qmckl_exit_code rc_ = qmckl_finalize_determinant(context);
-  //if (rc_ != QMCKL_SUCCESS) return rc_;
+  if (rc_ != QMCKL_SUCCESS) return rc_;
  }
 
 return QMCKL_SUCCESS;
@@ -412,7 +446,7 @@ qmckl_exit_code qmckl_set_determinant_det_num_beta(qmckl_context context, const
 qmckl_exit_code  qmckl_set_determinant_mo_index_alpha(qmckl_context context, const int64_t* mo_index_alpha) {
   <<pre2>>
 
-  int32_t mask = 1 << 5;
+  int32_t mask = 1 << 4;
 
   if (ctx->det.mo_index_alpha != NULL) {
     qmckl_exit_code rc = qmckl_free(context, ctx->det.mo_index_alpha);
@@ -444,7 +478,7 @@ qmckl_exit_code  qmckl_set_determinant_mo_index_alpha(qmckl_context context, con
 qmckl_exit_code  qmckl_set_determinant_mo_index_beta(qmckl_context context, const int64_t* mo_index_beta) {
   <<pre2>>
 
-  int32_t mask = 1 << 6;
+  int32_t mask = 1 << 5;
 
   if (ctx->det.mo_index_beta != NULL) {
     qmckl_exit_code rc = qmckl_free(context, ctx->det.mo_index_beta);
@@ -479,11 +513,50 @@ qmckl_exit_code  qmckl_set_determinant_mo_index_beta(qmckl_context context, cons
  computed to accelerate the calculations.
 
    #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
-qmckl_exit_code qmckl_finalize_basis(qmckl_context context);
+qmckl_exit_code qmckl_finalize_determinant(qmckl_context context);
    #+end_src
 
    #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
 qmckl_exit_code qmckl_finalize_determinant(qmckl_context context) {
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_CONTEXT,
+                           "qmckl_finalize_determinant",
+                           NULL);
+  }
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  assert (ctx != NULL);
+
+  qmckl_exit_code rc;
+  rc = qmckl_provide_det_vgl_alpha(context);
+  if (rc != QMCKL_SUCCESS) {
+    return qmckl_failwith( context,
+                           QMCKL_NOT_PROVIDED,
+                           "qmckl_finalize_determinant",
+                            NULL);
+  }
+  rc = qmckl_provide_det_vgl_beta(context);
+  if (rc != QMCKL_SUCCESS) {
+    return qmckl_failwith( context,
+                           QMCKL_NOT_PROVIDED,
+                           "qmckl_finalize_determinant",
+                            NULL);
+  }
+  rc = qmckl_provide_det_inv_matrix_alpha(context);
+  if (rc != QMCKL_SUCCESS) {
+    return qmckl_failwith( context,
+                           QMCKL_NOT_PROVIDED,
+                           "qmckl_finalize_determinant",
+                            NULL);
+  }
+  rc = qmckl_provide_det_inv_matrix_beta(context);
+  if (rc != QMCKL_SUCCESS) {
+    return qmckl_failwith( context,
+                           QMCKL_NOT_PROVIDED,
+                           "qmckl_finalize_determinant",
+                            NULL);
+  }
 }
    #+end_src
 
@@ -608,7 +681,7 @@ qmckl_exit_code qmckl_provide_det_vgl_alpha(qmckl_context context) {
   if (!ctx->det.provided) {
     return qmckl_failwith( context,
                            QMCKL_NOT_PROVIDED,
-                           "qmckl_mo_basis",
+                           "qmckl_determinant",
                            NULL);
   }
 
@@ -959,19 +1032,19 @@ integer function qmckl_compute_det_vgl_beta_f(context, &
       do imo = 1, beta_num
         mo_id = mo_index_beta(imo, iwalk, idet)
         ! Value
-        det_vgl_beta(imo, ielec, 1, iwalk, idet) = mo_vgl(mo_id, alpha_num + ielec, 1)
+        det_vgl_beta(imo, ielec, 1, iwalk, idet) = mo_vgl(mo_id, ielec, 1)
 
         ! Grad_x
-        det_vgl_beta(imo, ielec, 2, iwalk, idet) = mo_vgl(mo_id, alpha_num + ielec, 2)
+        det_vgl_beta(imo, ielec, 2, iwalk, idet) = mo_vgl(mo_id, ielec, 2)
 
         ! Grad_y
-        det_vgl_beta(imo, ielec, 3, iwalk, idet) = mo_vgl(mo_id, alpha_num + ielec, 3)
+        det_vgl_beta(imo, ielec, 3, iwalk, idet) = mo_vgl(mo_id, ielec, 3)
 
         ! Grad_z
-        det_vgl_beta(imo, ielec, 4, iwalk, idet) = mo_vgl(mo_id, alpha_num + ielec, 4)
+        det_vgl_beta(imo, ielec, 4, iwalk, idet) = mo_vgl(mo_id, ielec, 4)
 
         ! Lap
-        det_vgl_beta(imo, ielec, 5, iwalk, idet) = mo_vgl(mo_id, alpha_num + ielec, 5)
+        det_vgl_beta(imo, ielec, 5, iwalk, idet) = mo_vgl(mo_id, ielec, 5)
       end do
     end do
   end do
@@ -1280,7 +1353,7 @@ qmckl_exit_code qmckl_get_det_inv_matrix_beta(qmckl_context context, double * co
   rc = qmckl_provide_mo_vgl(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  rc = qmckl_provide_det_vgl_alpha(context);
+  rc = qmckl_provide_det_vgl_beta(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
   rc = qmckl_provide_det_inv_matrix_beta(context);
@@ -1496,7 +1569,7 @@ qmckl_exit_code qmckl_provide_det_inv_matrix_beta(qmckl_context context) {
 
       qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
       mem_info.size = ctx->det.walk_num * ctx->det.det_num_beta * 
-                      ctx->electron.up_num * ctx->electron.up_num * sizeof(double);
+                      ctx->electron.down_num * ctx->electron.down_num * sizeof(double);
       double* det_adj_matrix_beta = (double*) qmckl_malloc(context, mem_info);
 
       if (det_adj_matrix_beta == NULL) {

org/qmckl_local_energy.org

@@ -412,6 +412,16 @@ integer function qmckl_compute_kinetic_energy_f(context, walk_num, &
      return
   endif
 
+  print *,"Size of mo_vgl=",size(mo_vgl)
+  print *,"Size of mo_index_alpha=",size(mo_index_alpha)
+  print *,"Size of mo_index_beta=",size(mo_index_beta)
+  print *,"Size of det_adj_matrix_alpha=",size(det_adj_matrix_alpha)
+  print *,"Size of det_adj_matrix_beta=",size(det_adj_matrix_beta)
+  print *,"num_alpha = ",alpha_num
+  print *,"num_beta = ",beta_num
+  print *,"det_num_alpha = ",det_num_alpha
+  print *,"det_num_beta = ",det_num_beta
+  print *,"walk_num = ",walk_num
   e_kin = 0.0d0
   do idet = 1, det_num_alpha
   do iwalk = 1, walk_num
@@ -427,8 +437,9 @@ integer function qmckl_compute_kinetic_energy_f(context, walk_num, &
     do imo = 1, beta_num
     do ielec = 1, beta_num
       mo_id = mo_index_beta(ielec, iwalk, idet)
-      e_kin(iwalk) = e_kin(iwalk) - 0.5d0 * det_adj_matrix_beta(imo, ielec, iwalk, idet) * &
+      print *,mo_id, imo, ielec, iwalk, idet
-                                    mo_vgl(mo_id, ielec, 5)
+      e_kin(iwalk) = e_kin(iwalk) - 0.5d0 * det_adj_matrix_beta(imo, ielec, iwalk, idet) 
+      !                              mo_vgl(mo_id, ielec, 5)
     end do
     end do
   end do
@@ -1447,7 +1458,7 @@ qmckl_exit_code qmckl_provide_drift_vector(qmckl_context context) {
    | ~int64_t~       | ~mo_index_alpha[det_num_alpha][walk_num][alpha_num]~                  | in  | MO indices for electrons                  |
    | ~int64_t~       | ~mo_index_beta[det_num_beta][walk_num][beta_num]~                     | in  | MO indices for electrons                  |
    | ~int64_t~       | ~mo_num~                                                              | in  | Number of MOs                             |
-   | ~double~        | ~mo_vgl[5][elec_num][mo_num]~                               | in  | Value, gradients and Laplacian of the MOs |
+   | ~double~        | ~mo_vgl[5][elec_num][mo_num]~                                         | in  | Value, gradients and Laplacian of the MOs |
    | ~double~        | ~det_adj_matrix_alpha[det_num_alpha][walk_num][alpha_num][alpha_num]~ | in  | Value, gradients and Laplacian of the Det |
    | ~double~        | ~det_adj_matrix_beta[det_num_beta][walk_num][beta_num][beta_num]~     | in  | Value, gradients and Laplacian of the Det |
    | ~double~        | ~r_drift[walk_num][3]~                                                | out | Kinetic energy                            |

org/qmckl_mo.org

@@ -118,6 +118,26 @@ typedef struct qmckl_mo_basis_struct {
    Some values are initialized by default, and are not concerned by
    this mechanism.
 
+   #+begin_src c :comments org :tangle (eval h_private_func) 
+qmckl_exit_code qmckl_init_mo_basis(qmckl_context context);
+   #+end_src
+   
+   #+begin_src c :comments org :tangle (eval c)
+qmckl_exit_code qmckl_init_mo_basis(qmckl_context context) {
+
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return false;
+  }
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  assert (ctx != NULL);
+
+  ctx->mo_basis.uninitialized = (1 << 2) - 1;
+
+  return QMCKL_SUCCESS;
+}
+   #+end_src
+   
 ** Access functions
 
    #+begin_src c :comments org :tangle (eval h_func) :exports none
@@ -259,10 +279,9 @@ qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
 ctx->mo_basis.uninitialized &= ~mask;
 ctx->mo_basis.provided = (ctx->mo_basis.uninitialized == 0);
 if (ctx->mo_basis.provided) {
-  qmckl_exit_code rc_ = qmckl_finalize_basis(context);
+  qmckl_exit_code rc_ = qmckl_finalize_mo_basis(context);
   if (rc_ != QMCKL_SUCCESS) return rc_;
- }
+}
-
 return QMCKL_SUCCESS;
    #+end_src
 
@@ -319,6 +338,34 @@ qmckl_exit_code  qmckl_set_mo_basis_coefficient(qmckl_context context, const dou
  When the basis set is completely entered, other data structures are
  computed to accelerate the calculations.
 
+   #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
+qmckl_exit_code qmckl_finalize_mo_basis(qmckl_context context);
+   #+end_src
+
+   #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
+qmckl_exit_code qmckl_finalize_mo_basis(qmckl_context context) {
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_CONTEXT,
+                           "qmckl_finalize_mo_basis",
+                           NULL);
+  }
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  assert (ctx != NULL);
+
+  qmckl_exit_code rc;
+  rc = qmckl_provide_mo_vgl(context);
+  if (rc != QMCKL_SUCCESS) {
+    return qmckl_failwith( context,
+                           QMCKL_NOT_PROVIDED,
+                           "qmckl_finalize_mo_basis",
+                            NULL);
+  }
+  return rc;
+}
+   #+end_src
+
 * Computation
 
 ** Computation of MOs
@@ -378,6 +425,7 @@ qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context);
 qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context)
 {
 
+  qmckl_exit_code rc; 
   if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
     return QMCKL_NULL_CONTEXT;
   }
@@ -392,6 +440,14 @@ qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context)
                            NULL);
   }
 
+  rc = qmckl_provide_ao_vgl(context);
+  if (rc != QMCKL_SUCCESS) {
+    return qmckl_failwith( context,
+                           QMCKL_NOT_PROVIDED,
+                           "qmckl_ao_basis",
+                           NULL);
+  }
+
   if(!(ctx->electron.provided)) {
     return qmckl_failwith( context,
                            QMCKL_NOT_PROVIDED,
@@ -425,7 +481,6 @@ qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context)
       ctx->mo_basis.mo_vgl = mo_vgl;
     }
 
-    qmckl_exit_code rc; 
     rc = qmckl_compute_mo_basis_vgl(context,
                                     ctx->ao_basis.ao_num,
                                     ctx->mo_basis.mo_num,
@@ -479,10 +534,6 @@ integer function qmckl_compute_mo_basis_vgl_f(context, &
   double precision,dimension(:,:),allocatable    :: ao_vgl_big
   double precision              :: alpha, beta
   integer                       :: info_qmckl_dgemm_value
-  integer                       :: info_qmckl_dgemm_Gx   
-  integer                       :: info_qmckl_dgemm_Gy   
-  integer                       :: info_qmckl_dgemm_Gz   
-  integer                       :: info_qmckl_dgemm_lap  
   integer*8                     :: M, N, K, LDA, LDB, LDC, i,j
 
   integer*8 :: inucl, iprim, iwalk, ielec, ishell

org/qmckl_trexio.org

@@ -915,6 +915,7 @@ qmckl_trexio_read(const qmckl_context context, const char* file_name)
 
   #+begin_src c :tangle (eval c_test)
 #ifdef HAVE_TREXIO
+#define walk_num 2
 
 qmckl_exit_code rc;
 char fname[256];
@@ -927,6 +928,8 @@ char message[256];
 strncpy(fname, QMCKL_TEST_DIR,255);
 strncat(fname, "/chbrclf", 255);
 printf("Test file: %s\n", fname);
+
+rc = qmckl_set_electron_walk_num(context, walk_num);
 rc = qmckl_trexio_read(context, fname);
 
 if (rc != QMCKL_SUCCESS) {