From ca15c2866daece2dbbc715ab3dbe9fe6ea5d873a Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 17 Nov 2022 14:49:11 +0100 Subject: [PATCH] Added rescaling MO coefs --- org/qmckl_mo.org | 107 +++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 100 insertions(+), 7 deletions(-) diff --git a/org/qmckl_mo.org b/org/qmckl_mo.org index dd9eca4..1f05334 100644 --- a/org/qmckl_mo.org +++ b/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 ; iao_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];