Added rescaling MO coefs

org/qmckl_mo.org

@@ -422,7 +422,7 @@ bool qmckl_mo_basis_provided(const qmckl_context context) {
 
    #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
 interface
-  integer(c_int32_t) function qmckl_get_mo_basis_mo_num (context, &
+  integer(qmckl_exit_code) function qmckl_get_mo_basis_mo_num (context, &
        mo_num) bind(C)
     use, intrinsic :: iso_c_binding
     import
@@ -433,7 +433,7 @@ interface
 end interface
 
 interface
-  integer(c_int32_t) function qmckl_get_mo_basis_coefficient(context, &
+  integer(qmckl_exit_code) function qmckl_get_mo_basis_coefficient(context, &
        coefficient, size_max) bind(C)
     use, intrinsic :: iso_c_binding
     import
@@ -533,7 +533,7 @@ bool qmckl_mo_basis_select_mo (const qmckl_context context,
 
    #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
 interface
-  integer(c_int32_t) function qmckl_mo_basis_select_mo (context, &
+  integer(qmckl_exit_code) function qmckl_mo_basis_select_mo (context, &
        keep, size_max) bind(C)
     use, intrinsic :: iso_c_binding
     import
@@ -592,7 +592,7 @@ qmckl_get_mo_basis_mo_value(qmckl_context context,
 
     #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
   interface
-     integer(c_int32_t) function qmckl_get_mo_basis_mo_value (context, &
+     integer(qmckl_exit_code) function qmckl_get_mo_basis_mo_value (context, &
           mo_value, size_max) bind(C)
        use, intrinsic :: iso_c_binding
        import
@@ -659,7 +659,7 @@ qmckl_get_mo_basis_mo_value_inplace (qmckl_context context,
 
     #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
   interface
-     integer(c_int32_t) function qmckl_get_mo_basis_mo_value_inplace (context, &
+     integer(qmckl_exit_code) function qmckl_get_mo_basis_mo_value_inplace (context, &
           mo_value, size_max) bind(C)
        use, intrinsic :: iso_c_binding
        import
@@ -1074,7 +1074,7 @@ qmckl_get_mo_basis_mo_vgl(qmckl_context context,
 
     #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
   interface
-     integer(c_int32_t) function qmckl_get_mo_basis_mo_vgl (context, &
+     integer(qmckl_exit_code) function qmckl_get_mo_basis_mo_vgl (context, &
           mo_vgl, size_max) bind(C)
        use, intrinsic :: iso_c_binding
        import
@@ -1141,7 +1141,7 @@ qmckl_get_mo_basis_mo_vgl_inplace (qmckl_context context,
 
     #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
   interface
-     integer(c_int32_t) function qmckl_get_mo_basis_mo_vgl_inplace (context, &
+     integer(qmckl_exit_code) function qmckl_get_mo_basis_mo_vgl_inplace (context, &
           mo_vgl, size_max) bind(C)
        use, intrinsic :: iso_c_binding
        import
@@ -1541,6 +1541,80 @@ qmckl_compute_mo_basis_mo_vgl_hpc (const qmckl_context context,
 #endif
     #+end_src
 
+** Rescaling of MO coefficients
+
+   When evaluating Slater determinants, the value of the determinants
+   may get out of the range of double precision. A simple fix is to
+   rescale the MO coefficients to put back the determinants in the
+   correct range.
+   
+    #+begin_src c :comments org :tangle (eval h_func) :noweb yes
+qmckl_exit_code
+qmckl_mo_basis_rescale(qmckl_context context,
+                          const double scaling_factor);
+    #+end_src
+
+    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
+qmckl_exit_code
+qmckl_mo_basis_rescale(qmckl_context context,
+                          const double scaling_factor)
+{
+
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+      return qmckl_failwith( context,
+                             QMCKL_NULL_CONTEXT,
+                             "qmckl_mo_basis_rescale",
+                             NULL);
+  }
+
+  if (scaling_factor == 0.) {
+      return qmckl_failwith( context,
+                             QMCKL_INVALID_ARG_2,
+                             "qmckl_mo_basis_rescale",
+                             "scaling factor can't be zero");
+  }
+
+  qmckl_exit_code rc;
+
+  rc = qmckl_provide_mo_basis_mo_vgl(context);
+  if (rc != QMCKL_SUCCESS) return rc;
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
+  assert (ctx != NULL);
+
+  if (!ctx->mo_basis.provided) {
+    return qmckl_failwith( context,
+                           QMCKL_NOT_PROVIDED,
+                           "qmckl_mo_basis_rescale",
+                           NULL);
+  }
+
+  for (int64_t i=0 ; i<ctx->ao_basis.ao_num * ctx->mo_basis.mo_num ; ++i) {
+    ctx->mo_basis.coefficient[i] *= scaling_factor;
+    ctx->mo_basis.coefficient_t[i] *= scaling_factor;
+  }
+  rc = qmckl_context_touch(context);
+
+
+  return rc;
+}
+    #+end_src
+
+*** Fortran interface
+
+   #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
+interface
+  integer(qmckl_exit_code) function qmckl_mo_basis_rescale (context, &
+       scaling_factor) bind(C)
+    use, intrinsic :: iso_c_binding
+    import
+    implicit none
+    integer (c_int64_t) , intent(in), value :: context
+    real    (c_double)  , intent(in), value :: scaling_factor
+  end function qmckl_mo_basis_rescale
+end interface
+    #+end_src
+
 ** Test
 
    #+begin_src python :results output :exports none
@@ -1758,6 +1832,25 @@ for (int i=0 ; i< point_num; ++i) {
   }
  }
 
+rc = qmckl_mo_basis_rescale(context, 0.);
+assert (rc != QMCKL_SUCCESS);
+
+rc = qmckl_mo_basis_rescale(context, 2.);
+assert (rc == QMCKL_SUCCESS);
+
+
+rc = qmckl_get_mo_basis_mo_value(context, &(mo_value[0][0]), point_num*chbrclf_mo_num);
+assert (rc == QMCKL_SUCCESS);
+
+for (int i=0 ; i< point_num; ++i) {
+  for (int k=0 ; k< chbrclf_mo_num ; ++k) {
+   assert(fabs(2.*mo_vgl[i][0][k] - mo_value[i][k]) < 1.e-12) ;
+  }
+ }
+
+rc = qmckl_mo_basis_rescale(context, 0.5);
+assert (rc == QMCKL_SUCCESS);
+
 
 // Test overlap of MO
 //double point_x[10];