#+TITLE: Molecular Orbitals #+SETUPFILE: ../tools/theme.setup #+INCLUDE: ../tools/lib.org The molecular orbitals (MOs) are defined in the basis of AOs along with a AO to MO coefficient matrix \[C\]. Using these coefficients (e.g. from Hartree Fock SCF method) the MOs are defined as follows: \[ \phi_i(\mathbf{r}) = C_i * \chi_i (\mathbf{r}) \] In this section we demonstrate how to use the QMCkl specific DGEMM function to calculate the MOs. * Headers :noexport: #+begin_src elisp :noexport :results none (org-babel-lob-ingest "../tools/lib.org") #+end_src #+begin_src c :tangle (eval h_private_func) #ifndef QMCKL_MO_HPF #define QMCKL_MO_HPF #+end_src #+begin_src c :tangle (eval h_private_type) #ifndef QMCKL_MO_HPT #define QMCKL_MO_HPT #include #+end_src #+begin_src c :tangle (eval c_test) :noweb yes #include "qmckl.h" #include "assert.h" #ifdef HAVE_CONFIG_H #include "config.h" #endif #include #include #include "chbrclf.h" #include "qmckl_ao_private_func.h" #include "qmckl_mo_private_func.h" int main() { qmckl_context context; context = qmckl_context_create(); qmckl_exit_code rc; #+end_src #+begin_src c :tangle (eval c) #ifdef HAVE_CONFIG_H #include "config.h" #endif #ifdef HAVE_STDINT_H #include #elif HAVE_INTTYPES_H #include #endif #include #include #include #include #include "qmckl.h" #include "qmckl_context_private_type.h" #include "qmckl_memory_private_type.h" #include "qmckl_memory_private_func.h" #include "qmckl_ao_private_type.h" #include "qmckl_ao_private_func.h" #include "qmckl_mo_private_type.h" #include "qmckl_mo_private_func.h" #+end_src * Context The following arrays are stored in the context: |-----------------+--------------------+----------------------------------------| | ~mo_num~ | | Number of MOs | | ~coefficient~ | ~[mo_num][ao_num]~ | Orbital coefficients | | ~coefficient_t~ | ~[ao_num][mo_num]~ | Transposed of the Orbital coefficients | |-----------------+--------------------+----------------------------------------| Computed data: |-----------------+--------------------------+-------------------------------------------------------------------------------------| | ~mo_value~ | ~[point_num][mo_num]~ | Value of the MOs at point positions | | ~mo_value_date~ | ~uint64_t~ | Late modification date of the value of the MOs at point positions | | ~mo_vgl~ | ~[point_num][5][mo_num]~ | Value, gradients, Laplacian of the MOs at point positions | | ~mo_vgl_date~ | ~uint64_t~ | Late modification date of Value, gradients, Laplacian of the MOs at point positions | |-----------------+--------------------------+-------------------------------------------------------------------------------------| ** Data structure #+begin_src c :comments org :tangle (eval h_private_type) typedef struct qmckl_mo_basis_struct { int64_t mo_num; double * restrict coefficient; double * restrict coefficient_t; double * restrict mo_vgl; double * restrict mo_value; uint64_t mo_vgl_date; uint64_t mo_value_date; int32_t uninitialized; bool provided; } qmckl_mo_basis_struct; #+end_src The ~uninitialized~ integer contains one bit set to one for each initialization function which has not been called. It becomes equal to zero after all initialization functions have been called. The struct is then initialized and ~provided == true~. 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*) 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 qmckl_exit_code qmckl_get_mo_basis_mo_num (const qmckl_context context, int64_t* mo_num); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_mo_basis_mo_num (const qmckl_context context, int64_t* mo_num) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return qmckl_failwith( context, QMCKL_INVALID_CONTEXT, "qmckl_get_mo_basis_mo_num", NULL); } qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int32_t mask = 1; if ( (ctx->mo_basis.uninitialized & mask) != 0) { return qmckl_failwith( context, QMCKL_NOT_PROVIDED, "qmckl_get_mo_basis_mo_num", NULL); } assert (ctx->mo_basis.mo_num > (int64_t) 0); ,*mo_num = ctx->mo_basis.mo_num; return QMCKL_SUCCESS; } #+end_src #+begin_src c :comments org :tangle (eval h_func) :exports none qmckl_exit_code qmckl_get_mo_basis_coefficient (const qmckl_context context, double* const coefficient, const int64_t size_max); #+end_src #+begin_src c :comments org :tangle (eval c) :exports none qmckl_exit_code qmckl_get_mo_basis_coefficient (const qmckl_context context, double* const coefficient, const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return qmckl_failwith( context, QMCKL_INVALID_CONTEXT, "qmckl_get_mo_basis_coefficient", NULL); } qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int32_t mask = 1 << 1; if ( (ctx->ao_basis.uninitialized & mask) != 0) { return qmckl_failwith( context, QMCKL_NOT_PROVIDED, "qmckl_get_mo_basis_coefficient", NULL); } if (coefficient == NULL) { return qmckl_failwith( context, QMCKL_INVALID_ARG_2, "qmckl_get_mo_basis_coefficient", "NULL pointer"); } if (size_max < ctx->ao_basis.ao_num * ctx->mo_basis.mo_num) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_mo_basis_coefficient", "Array too small. Expected mo_num * ao_num"); } assert (ctx->mo_basis.coefficient != NULL); memcpy(coefficient, ctx->mo_basis.coefficient, ctx->ao_basis.ao_num * ctx->mo_basis.mo_num * sizeof(double)); return QMCKL_SUCCESS; } #+end_src When all the data for the AOs have been provided, the following function returns ~true~. #+begin_src c :comments org :tangle (eval h_func) bool qmckl_mo_basis_provided (const qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :exports none bool qmckl_mo_basis_provided(const qmckl_context context) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return false; } qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); return ctx->mo_basis.provided; } #+end_src *** Fortran interfaces #+begin_src f90 :tangle (eval fh_func) :comments org :exports none interface integer(c_int32_t) function qmckl_get_mo_basis_mo_num (context, & mo_num) bind(C) use, intrinsic :: iso_c_binding import implicit none integer (c_int64_t) , intent(in) , value :: context integer (c_int64_t) , intent(out) :: mo_num end function qmckl_get_mo_basis_mo_num end interface interface integer(c_int32_t) function qmckl_get_mo_basis_coefficient(context, & coefficient, size_max) bind(C) use, intrinsic :: iso_c_binding import implicit none integer (c_int64_t) , intent(in) , value :: context double precision, intent(out) :: coefficient(*) integer (c_int64_t) , intent(in), value :: size_max end function qmckl_get_mo_basis_coefficient end interface #+end_src ** Initialization functions To set the basis set, all the following functions need to be called. #+begin_src c :comments org :tangle (eval h_func) qmckl_exit_code qmckl_set_mo_basis_mo_num (qmckl_context context, const int64_t mo_num); qmckl_exit_code qmckl_set_mo_basis_coefficient (qmckl_context context, const double * coefficient); #+end_src #+NAME:pre #+begin_src c :exports none if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_context_struct* const ctx = (qmckl_context_struct*) context; if (mask != 0 && !(ctx->mo_basis.uninitialized & mask)) { return qmckl_failwith( context, QMCKL_ALREADY_SET, "qmckl_set_mo_*", NULL); } #+end_src #+NAME:post #+begin_src c :exports none 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_mo_basis(context); if (rc_ != QMCKL_SUCCESS) return rc_; } return QMCKL_SUCCESS; #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_set_mo_basis_mo_num(qmckl_context context, const int64_t mo_num) { int32_t mask = 1 ; <
>

  if (mo_num <= 0) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_set_mo_basis_mo_num",
                           "mo_num <= 0");
  }

  ctx->mo_basis.mo_num = mo_num;

  <>
}

qmckl_exit_code  qmckl_set_mo_basis_coefficient(qmckl_context context, const double* coefficient) {
  
  int32_t mask = 1 << 1;

  <
>

  if (ctx->mo_basis.coefficient != NULL) {
    qmckl_exit_code rc = qmckl_free(context, ctx->mo_basis.coefficient);
    if (rc != QMCKL_SUCCESS) {
      return qmckl_failwith( context, rc,
                             "qmckl_set_mo_basis_coefficient",
                             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);
  if (new_array == NULL) {
    return qmckl_failwith( context,
                           QMCKL_ALLOCATION_FAILED,
                           "qmckl_set_mo_basis_coefficient",
                           NULL);
  }

  memcpy(new_array, coefficient, mem_info.size);

  ctx->mo_basis.coefficient = new_array;

  <>
}

   #+end_src

 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*) context;
  assert (ctx != 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);
  if (new_array == NULL) {
    return qmckl_failwith( context,
                           QMCKL_ALLOCATION_FAILED,
                           "qmckl_finalize_mo_basis",
                           NULL);
  }

  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];
    }
  }

  ctx->mo_basis.coefficient_t = new_array;
  qmckl_exit_code rc = QMCKL_SUCCESS;
  return rc;
}
   #+end_src

* Computation

** Computation of MOs: values only

*** Get

    #+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code
qmckl_get_mo_basis_mo_value(qmckl_context context,
                            double* const mo_value,
                            const int64_t size_max);
    #+end_src

    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
qmckl_exit_code
qmckl_get_mo_basis_mo_value(qmckl_context context,
                            double* const mo_value,
                            const int64_t size_max)
{

  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return QMCKL_NULL_CONTEXT;
  }

  qmckl_exit_code rc;

  rc = qmckl_provide_ao_basis_ao_value(context);
  if (rc != QMCKL_SUCCESS) return rc;

  rc = qmckl_provide_mo_value(context);
  if (rc != QMCKL_SUCCESS) return rc;

  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
  assert (ctx != NULL);

  const int64_t sze = ctx->point.num * ctx->mo_basis.mo_num;
  if (size_max < sze) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_3,
                           "qmckl_get_mo_basis_mo_value",
                           "input array too small");
  }
  memcpy(mo_value, ctx->mo_basis.mo_value, sze * sizeof(double));

  return QMCKL_SUCCESS;
}
    #+end_src

    #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
  interface
     integer(c_int32_t) function qmckl_get_mo_basis_mo_value (context, &
          mo_value, size_max) bind(C)
       use, intrinsic :: iso_c_binding
       import
       implicit none

       integer (c_int64_t) , intent(in)  , value :: context
       double precision,     intent(out)         :: mo_value(*)
       integer (c_int64_t) , intent(in)  , value :: size_max
     end function qmckl_get_mo_basis_mo_value
  end interface
    #+end_src

    Uses the given array to compute the values.

    #+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code
qmckl_get_mo_basis_mo_value_inplace (qmckl_context context,
                                     double* const mo_value,
                                     const int64_t size_max);
    #+end_src

    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
qmckl_exit_code
qmckl_get_mo_basis_mo_value_inplace (qmckl_context context,
                                     double* const mo_value,
                                     const int64_t size_max)
{

  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_CONTEXT,
                           "qmckl_get_mo_basis_mo_value",
                           NULL);
  }

  qmckl_exit_code rc;

  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
  assert (ctx != NULL);

  const int64_t sze = ctx->mo_basis.mo_num * ctx->point.num;
  if (size_max < sze) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_3,
                           "qmckl_get_mo_basis_mo_value",
                           "input array too small");
  }

  rc = qmckl_context_touch(context);
  if (rc != QMCKL_SUCCESS) return rc;

  double* old_array = ctx->mo_basis.mo_value;

  ctx->mo_basis.mo_value = mo_value;

  rc = qmckl_provide_mo_value(context);
  if (rc != QMCKL_SUCCESS) return rc;

  ctx->mo_basis.mo_value = old_array;

  return QMCKL_SUCCESS;
}
    #+end_src

    #+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, &
          mo_value, size_max) bind(C)
       use, intrinsic :: iso_c_binding
       import
       implicit none
       integer (c_int64_t) , intent(in)  , value :: context
       double precision,     intent(out)         :: mo_value(*)
       integer (c_int64_t) , intent(in)  , value :: size_max
     end function qmckl_get_mo_basis_mo_value_inplace
  end interface
    #+end_src

*** Provide

    #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_provide_mo_value(qmckl_context context);
    #+end_src

    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
qmckl_exit_code qmckl_provide_mo_value(qmckl_context context)
{

  qmckl_exit_code rc;
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return QMCKL_NULL_CONTEXT;
  }

  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
  assert (ctx != NULL);

  if (!ctx->ao_basis.provided) {
    return qmckl_failwith( context,
                           QMCKL_NOT_PROVIDED,
                           "qmckl_ao_basis",
                           NULL);
  }

  rc = qmckl_provide_ao_basis_ao_value(context);
  if (rc != QMCKL_SUCCESS) {
    return qmckl_failwith( context,
                           QMCKL_NOT_PROVIDED,
                           "qmckl_ao_value",
                           NULL);
  }

  if (!ctx->mo_basis.provided) {
    return qmckl_failwith( context,
                           QMCKL_NOT_PROVIDED,
                           "qmckl_mo_basis",
                           NULL);
  }

  /* Compute if necessary */
  if (ctx->point.date > ctx->mo_basis.mo_value_date) {

    /* Allocate array */
    if (ctx->mo_basis.mo_value == NULL) {

      qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
      mem_info.size = ctx->point.num * ctx->mo_basis.mo_num * sizeof(double);
      double* mo_value = (double*) qmckl_malloc(context, mem_info);

      if (mo_value == NULL) {
        return qmckl_failwith( context,
                               QMCKL_ALLOCATION_FAILED,
                               "qmckl_mo_basis_mo_value",
                               NULL);
      }
      ctx->mo_basis.mo_value = mo_value;
    }

    if (ctx->mo_basis.mo_vgl_date == ctx->point.date) {

      // mo_vgl has been computed at this step: Just copy the data.
      
      double * v = &(ctx->mo_basis.mo_value[0]);
      double * vgl = &(ctx->mo_basis.mo_vgl[0]);
      for (int i=0 ; ipoint.num ; ++i) {
        for (int k=0 ; kmo_basis.mo_num ; ++k) {
          v[k] = vgl[k];
        }
        v   += ctx->mo_basis.mo_num;
        vgl += ctx->mo_basis.mo_num * 5;
      }

    } else {

      rc = qmckl_compute_mo_basis_mo_value(context,
                                           ctx->ao_basis.ao_num,
                                           ctx->mo_basis.mo_num,
                                           ctx->point.num,
                                           ctx->mo_basis.coefficient_t,
                                           ctx->ao_basis.ao_value,
                                           ctx->mo_basis.mo_value);

      if (rc != QMCKL_SUCCESS) {
        return rc;
      }

    }

    ctx->mo_basis.mo_value_date = ctx->date;
  }

  return QMCKL_SUCCESS;
}
    #+end_src

*** Compute
   :PROPERTIES:
   :Name:     qmckl_compute_mo_basis_mo_value
   :CRetType: qmckl_exit_code
   :FRetType: qmckl_exit_code
   :END:

    #+NAME: qmckl_mo_basis_mo_value_args
    | Variable        | Type                        | In/Out | Description                                     |
    |-----------------+-----------------------------+--------+-------------------------------------------------|
    | ~context~       | ~qmckl_context~             | in     | Global state                                    |
    | ~ao_num~        | ~int64_t~                   | in     | Number of AOs                                   |
    | ~mo_num~        | ~int64_t~                   | in     | Number of MOs                                   |
    | ~point_num~     | ~int64_t~                   | in     | Number of points                                |
    | ~coefficient_t~ | ~double[mo_num][ao_num]~    | in     | Transpose of the AO to MO transformation matrix |
    | ~ao_value~      | ~double[point_num][ao_num]~ | in     | Value of the AOs                                |
    | ~mo_value~      | ~double[point_num][mo_num]~ | out    | Value of the MOs                                |


    The matrix of AO values is very sparse, so we use a sparse-dense
    matrix multiplication instead of a dgemm, as exposed in
    https://dx.doi.org/10.1007/978-3-642-38718-0_14.



    #+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_mo_basis_mo_value_doc_f(context, &
     ao_num, mo_num, point_num, &
     coefficient_t, ao_value, mo_value) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer*8             , intent(in)  :: ao_num, mo_num
  integer*8             , intent(in)  :: point_num
  double precision      , intent(in)  :: ao_value(ao_num,point_num)
  double precision      , intent(in)  :: coefficient_t(mo_num,ao_num)
  double precision      , intent(out) :: mo_value(mo_num,point_num)
  integer*8 :: i,j,k
  double precision :: c1, c2, c3, c4, c5

  integer*8 :: LDA, LDB, LDC

  info = QMCKL_SUCCESS
  if (.True.)  then    ! fast algorithm
     do j=1,point_num
        mo_value(:,j) = 0.d0
        do k=1,ao_num
           if (ao_value(k,j) /= 0.d0) then
              c1 = ao_value(k,j)
              do i=1,mo_num
                 mo_value(i,j) = mo_value(i,j) + coefficient_t(i,k) * c1
              end do
           end if
        end do
     end do
     
  else ! dgemm

    LDA = size(coefficient_t,1)
    LDB = size(ao_value,1) 
    LDC = size(mo_value,1)

    info = qmckl_dgemm(context,'N', 'N', mo_num, point_num, ao_num, 1.d0,     &
                                    coefficient_t, LDA, ao_value, LDB, &
                                    0.d0, mo_value, LDC)

  end if

end function qmckl_compute_mo_basis_mo_value_doc_f
    #+end_src

    #+CALL: generate_c_header(table=qmckl_mo_basis_mo_value_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_value"))

   #+RESULTS:
   #+begin_src c :tangle (eval h_func) :comments org
   qmckl_exit_code qmckl_compute_mo_basis_mo_value (
         const qmckl_context context,
         const int64_t ao_num,
         const int64_t mo_num,
         const int64_t point_num,
         const double* coefficient_t,
         const double* ao_value,
         double* const mo_value );
   #+end_src

   #+CALL: generate_c_header(table=qmckl_mo_basis_mo_value_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_value_doc"))

   #+RESULTS:
   #+begin_src c :tangle (eval h_func) :comments org
   qmckl_exit_code qmckl_compute_mo_basis_mo_value_doc (
         const qmckl_context context,
         const int64_t ao_num,
         const int64_t mo_num,
         const int64_t point_num,
         const double* coefficient_t,
         const double* ao_value,
         double* const mo_value );
   #+end_src

   #+CALL: generate_c_interface(table=qmckl_mo_basis_mo_value_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_value_doc"))

    #+RESULTS:
    #+begin_src f90 :tangle (eval f) :comments org :exports none
    integer(c_int32_t) function qmckl_compute_mo_basis_mo_value_doc &
        (context, ao_num, mo_num, point_num, coefficient_t, ao_value, mo_value) &
        bind(C) result(info)

      use, intrinsic :: iso_c_binding
      implicit none

      integer (c_int64_t) , intent(in)  , value :: context
      integer (c_int64_t) , intent(in)  , value :: ao_num
      integer (c_int64_t) , intent(in)  , value :: mo_num
      integer (c_int64_t) , intent(in)  , value :: point_num
      real    (c_double ) , intent(in)          :: coefficient_t(ao_num,mo_num)
      real    (c_double ) , intent(in)          :: ao_value(ao_num,point_num)
      real    (c_double ) , intent(out)         :: mo_value(mo_num,point_num)

      integer(c_int32_t), external :: qmckl_compute_mo_basis_mo_value_doc_f
      info = qmckl_compute_mo_basis_mo_value_doc_f &
             (context, ao_num, mo_num, point_num, coefficient_t, ao_value, mo_value)

    end function qmckl_compute_mo_basis_mo_value_doc
    #+end_src

    #+begin_src c :tangle (eval c) :comments org
qmckl_exit_code
qmckl_compute_mo_basis_mo_value (const qmckl_context context,
                                 const int64_t ao_num,
                                 const int64_t mo_num,
                                 const int64_t point_num,
                                 const double* coefficient_t,
                                 const double* ao_value,
                                 double* const mo_value )
{
#ifdef HAVE_HPC
  return qmckl_compute_mo_basis_mo_value_hpc (context, ao_num, mo_num, point_num, coefficient_t, ao_value, mo_value);
#else
  return qmckl_compute_mo_basis_mo_value_doc (context, ao_num, mo_num, point_num, coefficient_t, ao_value, mo_value);
#endif
}
    #+end_src

*** HPC version


    #+begin_src c :tangle (eval h_func) :comments org
#ifdef HAVE_HPC
qmckl_exit_code
qmckl_compute_mo_basis_mo_value_hpc (const qmckl_context context,
                                     const int64_t ao_num,
                                     const int64_t mo_num,
                                     const int64_t point_num,
                                     const double* coefficient_t,
                                     const double* ao_value,
                                     double* const mo_value );
#endif
    #+end_src

    #+begin_src c :tangle (eval c) :comments org
#ifdef HAVE_HPC
qmckl_exit_code
qmckl_compute_mo_basis_mo_value_hpc (const qmckl_context context,
                                     const int64_t ao_num,
                                     const int64_t mo_num,
                                     const int64_t point_num,
                                     const double* restrict coefficient_t,
                                     const double* restrict ao_value,
                                     double* restrict const mo_value )
{
  assert (context != QMCKL_NULL_CONTEXT);

#ifdef HAVE_OPENMP
  #pragma omp parallel for
#endif
  for (int64_t ipoint=0 ; ipoint < point_num ; ++ipoint) {
    double* restrict const vgl1  = &(mo_value[ipoint*mo_num]);
    const double* restrict avgl1 = &(ao_value[ipoint*ao_num]);

    for (int64_t i=0 ; ipoint.num * 5 * ctx->mo_basis.mo_num;
  if (size_max < sze) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_3,
                           "qmckl_get_mo_basis_mo_vgl",
                           "input array too small");
  }
  memcpy(mo_vgl, ctx->mo_basis.mo_vgl, sze * sizeof(double));

  return QMCKL_SUCCESS;
}
    #+end_src

    #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
  interface
     integer(c_int32_t) function qmckl_get_mo_basis_mo_vgl (context, &
          mo_vgl, size_max) bind(C)
       use, intrinsic :: iso_c_binding
       import
       implicit none

       integer (c_int64_t) , intent(in)  , value :: context
       double precision,     intent(out)         :: mo_vgl(*)
       integer (c_int64_t) , intent(in)  , value :: size_max
     end function qmckl_get_mo_basis_mo_vgl
  end interface
    #+end_src

    Uses the given array to compute the VGL.

    #+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code
qmckl_get_mo_basis_mo_vgl_inplace (qmckl_context context,
                                   double* const mo_vgl,
                                   const int64_t size_max);
    #+end_src

    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
qmckl_exit_code
qmckl_get_mo_basis_mo_vgl_inplace (qmckl_context context,
                                   double* const mo_vgl,
                                   const int64_t size_max)
{

  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_CONTEXT,
                           "qmckl_get_mo_basis_mo_vgl",
                           NULL);
  }

  qmckl_exit_code rc;

  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
  assert (ctx != NULL);

  const int64_t sze = ctx->mo_basis.mo_num * 5 * ctx->point.num;
  if (size_max < sze) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_3,
                           "qmckl_get_mo_basis_mo_vgl",
                           "input array too small");
  }

  rc = qmckl_context_touch(context);
  if (rc != QMCKL_SUCCESS) return rc;

  double* old_array = ctx->mo_basis.mo_vgl;

  ctx->mo_basis.mo_vgl = mo_vgl;

  rc = qmckl_provide_mo_vgl(context);
  if (rc != QMCKL_SUCCESS) return rc;

  ctx->mo_basis.mo_vgl = old_array;

  return QMCKL_SUCCESS;
}
    #+end_src

    #+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, &
          mo_vgl, size_max) bind(C)
       use, intrinsic :: iso_c_binding
       import
       implicit none
       integer (c_int64_t) , intent(in)  , value :: context
       double precision,     intent(out)         :: mo_vgl(*)
       integer (c_int64_t) , intent(in)  , value :: size_max
     end function qmckl_get_mo_basis_mo_vgl_inplace
  end interface
    #+end_src

*** Provide

    #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context);
    #+end_src

    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
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;
  }

  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
  assert (ctx != NULL);

  if (!ctx->ao_basis.provided) {
    return qmckl_failwith( context,
                           QMCKL_NOT_PROVIDED,
                           "qmckl_ao_basis",
                           NULL);
  }

  rc = qmckl_provide_ao_vgl(context);
  if (rc != QMCKL_SUCCESS) {
    return qmckl_failwith( context,
                           QMCKL_NOT_PROVIDED,
                           "qmckl_ao_basis",
                           NULL);
  }

  if (!ctx->mo_basis.provided) {
    return qmckl_failwith( context,
                           QMCKL_NOT_PROVIDED,
                           "qmckl_mo_basis",
                           NULL);
  }

  /* Compute if necessary */
  if (ctx->point.date > ctx->mo_basis.mo_vgl_date) {

    /* Allocate array */
    if (ctx->mo_basis.mo_vgl == NULL) {

      qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
      mem_info.size = 5 * ctx->point.num * ctx->mo_basis.mo_num * sizeof(double);
      double* mo_vgl = (double*) qmckl_malloc(context, mem_info);

      if (mo_vgl == NULL) {
        return qmckl_failwith( context,
                               QMCKL_ALLOCATION_FAILED,
                               "qmckl_mo_basis_mo_vgl",
                               NULL);
      }
      ctx->mo_basis.mo_vgl = mo_vgl;
    }

    rc = qmckl_compute_mo_basis_mo_vgl(context,
                                       ctx->ao_basis.ao_num,
                                       ctx->mo_basis.mo_num,
                                       ctx->point.num,
                                       ctx->mo_basis.coefficient_t,
                                       ctx->ao_basis.ao_vgl,
                                       ctx->mo_basis.mo_vgl);
    if (rc != QMCKL_SUCCESS) {
      return rc;
    }

    ctx->mo_basis.mo_vgl_date = ctx->date;
  }

  return QMCKL_SUCCESS;
}
    #+end_src

*** Compute
   :PROPERTIES:
   :Name:     qmckl_compute_mo_basis_mo_vgl
   :CRetType: qmckl_exit_code
   :FRetType: qmckl_exit_code
   :END:

    #+NAME: qmckl_mo_basis_mo_vgl_args
    | Variable            | Type                           | In/Out | Description                                     |
    |---------------------+--------------------------------+--------+-------------------------------------------------|
    | ~context~           | ~qmckl_context~                | in     | Global state                                    |
    | ~ao_num~            | ~int64_t~                      | in     | Number of AOs                                   |
    | ~mo_num~            | ~int64_t~                      | in     | Number of MOs                                   |
    | ~point_num~         | ~int64_t~                      | in     | Number of points                                |
    | ~coefficient_t~     | ~double[mo_num][ao_num]~       | in     | Transpose of the AO to MO transformation matrix |
    | ~ao_vgl~            | ~double[point_num][5][ao_num]~ | in     | Value, gradients and Laplacian of the AOs       |
    | ~mo_vgl~            | ~double[point_num][5][mo_num]~ | out    | Value, gradients and Laplacian of the MOs       |


    The matrix of AO values is very sparse, so we use a sparse-dense
    matrix multiplication instead of a dgemm, as exposed in
    https://dx.doi.org/10.1007/978-3-642-38718-0_14.



    #+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_mo_basis_mo_vgl_doc_f(context, &
     ao_num, mo_num, point_num, &
     coefficient_t, ao_vgl, mo_vgl) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer*8             , intent(in)  :: ao_num, mo_num
  integer*8             , intent(in)  :: point_num
  double precision      , intent(in)  :: ao_vgl(ao_num,5,point_num)
  double precision      , intent(in)  :: coefficient_t(mo_num,ao_num)
  double precision      , intent(out) :: mo_vgl(mo_num,5,point_num)
  integer*8 :: i,j,k
  double precision :: c1, c2, c3, c4, c5

  do j=1,point_num
     mo_vgl(:,:,j) = 0.d0
     do k=1,ao_num
        if (ao_vgl(k,1,j) /= 0.d0) then
           c1 = ao_vgl(k,1,j)
           c2 = ao_vgl(k,2,j)
           c3 = ao_vgl(k,3,j)
           c4 = ao_vgl(k,4,j)
           c5 = ao_vgl(k,5,j)
           do i=1,mo_num
              mo_vgl(i,1,j) = mo_vgl(i,1,j) + coefficient_t(i,k) * c1
              mo_vgl(i,2,j) = mo_vgl(i,2,j) + coefficient_t(i,k) * c2
              mo_vgl(i,3,j) = mo_vgl(i,3,j) + coefficient_t(i,k) * c3
              mo_vgl(i,4,j) = mo_vgl(i,4,j) + coefficient_t(i,k) * c4
              mo_vgl(i,5,j) = mo_vgl(i,5,j) + coefficient_t(i,k) * c5
           end do
        end if
     end do
  end do
  info = QMCKL_SUCCESS

! info = qmckl_dgemm(context,'N', 'N', mo_num, point_num, ao_num, 1.d0, &
!      coefficient_t, int(size(coefficient_t,1),8),      &
!      ao_vgl, int(size(ao_vgl,1),8), 0.d0,                  &
!      mo_vgl, int(size(mo_vgl,1),8))

end function qmckl_compute_mo_basis_mo_vgl_doc_f
    #+end_src

    #+CALL: generate_c_header(table=qmckl_mo_basis_mo_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_vgl"))

   #+RESULTS:
   #+begin_src c :tangle (eval h_func) :comments org
   qmckl_exit_code qmckl_compute_mo_basis_mo_vgl (
         const qmckl_context context,
         const int64_t ao_num,
         const int64_t mo_num,
         const int64_t point_num,
         const double* coefficient_t,
         const double* ao_vgl,
         double* const mo_vgl );
   #+end_src

   #+CALL: generate_c_header(table=qmckl_mo_basis_mo_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_vgl_doc"))

   #+RESULTS:
   #+begin_src c :tangle (eval h_func) :comments org
   qmckl_exit_code qmckl_compute_mo_basis_mo_vgl_doc (
         const qmckl_context context,
         const int64_t ao_num,
         const int64_t mo_num,
         const int64_t point_num,
         const double* coefficient_t,
         const double* ao_vgl,
         double* const mo_vgl );
   #+end_src

   #+CALL: generate_c_interface(table=qmckl_mo_basis_mo_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_vgl_doc"))

    #+RESULTS:
    #+begin_src f90 :tangle (eval f) :comments org :exports none
    integer(c_int32_t) function qmckl_compute_mo_basis_mo_vgl_doc &
        (context, ao_num, mo_num, point_num, coefficient_t, ao_vgl, mo_vgl) &
        bind(C) result(info)

      use, intrinsic :: iso_c_binding
      implicit none

      integer (c_int64_t) , intent(in)  , value :: context
      integer (c_int64_t) , intent(in)  , value :: ao_num
      integer (c_int64_t) , intent(in)  , value :: mo_num
      integer (c_int64_t) , intent(in)  , value :: point_num
      real    (c_double ) , intent(in)          :: coefficient_t(ao_num,mo_num)
      real    (c_double ) , intent(in)          :: ao_vgl(ao_num,5,point_num)
      real    (c_double ) , intent(out)         :: mo_vgl(mo_num,5,point_num)

      integer(c_int32_t), external :: qmckl_compute_mo_basis_mo_vgl_doc_f
      info = qmckl_compute_mo_basis_mo_vgl_doc_f &
             (context, ao_num, mo_num, point_num, coefficient_t, ao_vgl, mo_vgl)

    end function qmckl_compute_mo_basis_mo_vgl_doc
    #+end_src

    #+begin_src c :tangle (eval c) :comments org
qmckl_exit_code
qmckl_compute_mo_basis_mo_vgl (const qmckl_context context,
                            const int64_t ao_num,
                            const int64_t mo_num,
                            const int64_t point_num,
                            const double* coefficient_t,
                            const double* ao_vgl,
                            double* const mo_vgl )
{
#ifdef HAVE_HPC
  return qmckl_compute_mo_basis_mo_vgl_hpc (context, ao_num, mo_num, point_num, coefficient_t, ao_vgl, mo_vgl);
#else
  return qmckl_compute_mo_basis_mo_vgl_doc (context, ao_num, mo_num, point_num, coefficient_t, ao_vgl, mo_vgl);
#endif
}
    #+end_src

*** HPC version


    #+begin_src c :tangle (eval h_func) :comments org
#ifdef HAVE_HPC
qmckl_exit_code
qmckl_compute_mo_basis_mo_vgl_hpc (const qmckl_context context,
                                   const int64_t ao_num,
                                   const int64_t mo_num,
                                   const int64_t point_num,
                                   const double* coefficient_t,
                                   const double* ao_vgl,
                                   double* const mo_vgl );
#endif
    #+end_src

    #+begin_src c :tangle (eval c) :comments org
#ifdef HAVE_HPC
qmckl_exit_code
qmckl_compute_mo_basis_mo_vgl_hpc (const qmckl_context context,
                                   const int64_t ao_num,
                                   const int64_t mo_num,
                                   const int64_t point_num,
                                   const double* restrict coefficient_t,
                                   const double* restrict ao_vgl,
                                   double* restrict const mo_vgl )
{
  assert (context != QMCKL_NULL_CONTEXT);

#ifdef HAVE_OPENMP
  #pragma omp parallel for
#endif
  for (int64_t ipoint=0 ; ipoint < point_num ; ++ipoint) {
    double* restrict const vgl1 = &(mo_vgl[ipoint*5*mo_num]);
    double* restrict const vgl2 =  vgl1 + mo_num;
    double* restrict const vgl3 =  vgl1 + (mo_num << 1);
    double* restrict const vgl4 =  vgl1 + (mo_num << 1) + mo_num;
    double* restrict const vgl5 =  vgl1 + (mo_num << 2);

    const double* restrict avgl1 = &(ao_vgl[ipoint*5*ao_num]);
    const double* restrict avgl2 = avgl1 + ao_num;
    const double* restrict avgl3 = avgl1 + (ao_num << 1);
    const double* restrict avgl4 = avgl1 + (ao_num << 1) + ao_num;
    const double* restrict avgl5 = avgl1 + (ao_num << 2);

    for (int64_t i=0 ; i