From 5666b52d92e590b6759927db1f92f8fdfc618ae8 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 17 Aug 2022 13:42:17 +0200 Subject: [PATCH] Added qmckl_mo_basis_select_mo --- org/qmckl_mo.org | 165 +++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 152 insertions(+), 13 deletions(-) diff --git a/org/qmckl_mo.org b/org/qmckl_mo.org index d88e8c8..b970152 100644 --- a/org/qmckl_mo.org +++ b/org/qmckl_mo.org @@ -254,6 +254,15 @@ qmckl_exit_code qmckl_finalize_mo_basis(qmckl_context context) { qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); + if (ctx->mo_basis.coefficient_t != NULL) { + qmckl_exit_code rc = qmckl_free(context, ctx->mo_basis.coefficient_t); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_finalize_mo_basis", + NULL); + } + } + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = ctx->ao_basis.ao_num * ctx->mo_basis.mo_num * sizeof(double); double* new_array = (double*) qmckl_malloc(context, mem_info); @@ -266,15 +275,6 @@ qmckl_exit_code qmckl_finalize_mo_basis(qmckl_context context) { assert (ctx->mo_basis.coefficient != NULL); - if (ctx->mo_basis.coefficient_t != NULL) { - qmckl_exit_code rc = qmckl_free(context, ctx->mo_basis.coefficient); - if (rc != QMCKL_SUCCESS) { - return qmckl_failwith( context, rc, - "qmckl_finalize_mo_basis", - NULL); - } - } - for (int64_t i=0 ; iao_basis.ao_num ; ++i) { for (int64_t j=0 ; jmo_basis.mo_num ; ++j) { new_array[i*ctx->mo_basis.mo_num + j] = ctx->mo_basis.coefficient[j*ctx->ao_basis.ao_num + i]; @@ -282,8 +282,23 @@ qmckl_exit_code qmckl_finalize_mo_basis(qmckl_context context) { } ctx->mo_basis.coefficient_t = new_array; - qmckl_exit_code rc = QMCKL_SUCCESS; - return rc; + + qmckl_exit_code rc; + if (ctx->mo_basis.mo_vgl != NULL) { + rc = qmckl_free(context, ctx->mo_basis.mo_vgl); + if (rc != QMCKL_SUCCESS) return rc; + ctx->mo_basis.mo_vgl = NULL; + ctx->mo_basis.mo_vgl_date = 0; + } + + if (ctx->mo_basis.mo_value != NULL) { + rc = qmckl_free(context, ctx->mo_basis.mo_value); + if (rc != QMCKL_SUCCESS) return rc; + ctx->mo_basis.mo_value = NULL; + ctx->mo_basis.mo_value_date = 0; + } + + return QMCKL_SUCCESS; } #+end_src @@ -403,7 +418,6 @@ bool qmckl_mo_basis_provided(const qmckl_context context) { #+end_src - *** Fortran interfaces #+begin_src f90 :tangle (eval fh_func) :comments org :exports none @@ -432,6 +446,105 @@ end interface #+end_src +** Update + + Useless MOs can be removed, for instance virtual MOs in a single + determinant calculation. + + To select a subset of MOs that will be kept, create an array of + integers of size =mo_num=. If the integer is zero, the MO is dropped, + otherwise it is kept. + + #+begin_src c :comments org :tangle (eval h_func) +bool qmckl_mo_basis_select_mo (const qmckl_context context, + const int32_t* keep, + const int64_t size_max); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :exports none +bool qmckl_mo_basis_select_mo (const qmckl_context context, + const int32_t* keep, + const int64_t size_max) +{ + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return qmckl_failwith( context, + QMCKL_INVALID_CONTEXT, + "qmckl_get_mo_basis_select_mo", + NULL); + } + + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; + assert (ctx != NULL); + + if ( !(qmckl_mo_basis_provided(context)) ) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_get_mo_basis_select_mo", + NULL); + } + + if (keep == NULL) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_get_mo_basis_select_mo", + "NULL pointer"); + } + + const int64_t mo_num = ctx->mo_basis.mo_num; + const int64_t ao_num = ctx->ao_basis.ao_num; + + if (size_max < mo_num) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_get_mo_basis_select_mo", + "Array too small: expected mo_num."); + } + + int64_t mo_num_new = 0; + for (int64_t i=0 ; imo_basis.coefficient[i*ao_num]), ao_num*sizeof(double)); + ++k; + } + } + + qmckl_exit_code rc = qmckl_free(context, ctx->mo_basis.coefficient); + if (rc != QMCKL_SUCCESS) return rc; + + ctx->mo_basis.coefficient = coefficient; + ctx->mo_basis.mo_num = mo_num_new; + + rc = qmckl_finalize_mo_basis(context); + return rc; +} + + #+end_src + +*** Fortran interface + + #+begin_src f90 :tangle (eval fh_func) :comments org :exports none +interface + integer(c_int32_t) function qmckl_mo_basis_select_mo (context, & + keep, size_max) bind(C) + use, intrinsic :: iso_c_binding + import + implicit none + integer (c_int64_t) , intent(in), value :: context + integer (c_int32_t) , intent(in) :: keep(size_max) + integer (c_int64_t) , intent(in), value :: size_max + end function qmckl_mo_basis_select_mo +end interface + #+end_src + * Computation ** Computation of MOs: values only @@ -1605,7 +1718,7 @@ rc = qmckl_get_ao_basis_ao_vgl(context, &(ao_vgl[0][0][0]), assert (rc == QMCKL_SUCCESS); /* Set up MO data */ -const int64_t mo_num = chbrclf_mo_num; +int64_t mo_num = chbrclf_mo_num; rc = qmckl_set_mo_basis_mo_num(context, mo_num); assert (rc == QMCKL_SUCCESS); @@ -1701,6 +1814,32 @@ printf(" mo_vgl mo_vgl[1][26][223] %25.15e\n", mo_vgl[2][1][3]); printf(" mo_vgl mo_vgl[0][26][224] %25.15e\n", mo_vgl[2][0][3]); printf(" mo_vgl mo_vgl[1][26][224] %25.15e\n", mo_vgl[2][1][3]); printf("\n"); + + +/* Check selection of MOs */ + + int32_t keep[mo_num]; + for (int i=0 ; i