qmckl/org/qmckl_local_energy.org

#+TITLE: Local Energy
#+SETUPFILE: ../tools/theme.setup
#+INCLUDE: ../tools/lib.org


Here we calculate the final expectation value of the
local energy \[E_L\] as the sum of the kinetic energy and potential energy.

\[
E_L = KE + PE
\]

Where the kinetic energy is given as:

\[
KE = -\frac{1}{2}\frac{\bigtriangleup \Psi}{\Psi} 
\]

The laplacian of the wavefunction in the single-determinant
case is given as follows:

\[
\frac{\bigtriangleup \Psi(r)}{\Psi(r)} = \sum_{j=1}^{N_e} \bigtriangleup \Phi_j(r_i) D_{ji}^{-1}(r)
\]

The potential energy is the sum of all the following terms

\[
PE = \mathcal{V}_{ee} + \mathcal{V}_{en} + \mathcal{V}_{nn}
\]

The potential for is calculated as the sum of single electron
contributions.

\[
\mathcal{V}_{ee} = \sum_{i=1}^{N_e}\sum_{j<i} \frac{1}{r_{ij}}
\]

\[
\mathcal{V}_{en} = - \sum_{i=1}^{N_e}\sum_{A=1}^{N_n}\frac{Z_A}{r_{iA}}
\]

\[
\mathcal{V}_{nn} = \sum_{A=1}^{N_n}\sum_{B<A}\frac{Z_A Z_B}{r_{AB}}
\]

The remaining quantities that are required for the calculation of
a single Monte-Carlo step are as follows:

1. Drift Vector - \[F(x)\]

\[
F(x) = 2\frac{\nabla \Psi(r)}{\Psi(r)}
\]

2. Diffusion move - \[y\]

\[
y = x + D F(x) \delta t + \chi
\]

Where \[\chi\] is a random number with gaussian distribution centered at 0 
and width of \[2D\delta t\]. Here \[D\] is the drift parameter.

3. Acceptance probability - \[min\left[1, q(y,x)\right]\]

\[
q(y,x) = Exp\left[\frac{\delta r}{2} ( F(x) + F(y)) - \frac{D\delta t}{4}\left( F(x)^2 - F(y)^2\right)\right]
\]

With these quantities, a single determinant VMC simulation can be carried out.

* Headers                                                          :noexport:
  #+begin_src elisp :noexport :results none
(org-babel-lob-ingest "../tools/lib.org")
  #+end_src

  #+begin_src c :tangle (eval h_private_type)
#ifndef QMCKL_LOCAL_ENERGY_HPT
#define QMCKL_LOCAL_ENERGY_HPT

#include <stdbool.h>
  #+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 <stdio.h>
#include <math.h>
#include "chbrclf.h"
#include "qmckl_ao_private_func.h"
#include "qmckl_mo_private_func.h"
#include "qmckl_determinant_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 <stdint.h>
#elif HAVE_INTTYPES_H
#include <inttypes.h>
#endif

#include <stdlib.h>
#include <string.h>
#include <stdbool.h>
#include <assert.h>

#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"
#include "qmckl_determinant_private_type.h"
#include "qmckl_determinant_private_func.h"
  #+end_src

* Context

  The following arrays are stored in the context:

  |   |

  Computed data:
  
  |--------------+-----------------+----------------------------|
  | ~e_kin~      | ~[walk_num]~    | total kinetic energy       |
  | ~e_pot~      | ~[walk_num]~    | total potential energy     |
  | ~e_local~    | ~[walk_num]~    | local energy               |
  | ~r_drift~    | ~[3][walk_num]~ | The drift vector           |
  | ~y_move~     | ~[3][walk_num]~ | The diffusion move         |
  | ~accep_prob~ | ~[walk_num]~    | The acceptance probability |
  |--------------+-----------------+----------------------------|
  
** Data structure
   
   #+begin_src c :comments org :tangle (eval h_private_type)
typedef struct qmckl_local_energy_struct {
  double  * e_kin;
  double  * e_pot;
  double  * e_local;
  double  * accep_prob;
  double  * r_drift;
  double  * y_move;
  int64_t   e_kin_date;
  int64_t   e_pot_date;
  int64_t   e_local_date;
  int64_t   accep_prob_date;
  int64_t   r_drift_date;
  int64_t   y_move_date;

  int32_t   uninitialized;
  bool      provided;
} qmckl_local_energy_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.

* Computation
** Kinetic energy
   :PROPERTIES:
   :Name:     qmckl_compute_kinetic_energy
   :CRetType: qmckl_exit_code
   :FRetType: qmckl_exit_code
   :END:

Where the kinetic energy is given as:

\[
KE = -\frac{1}{2}\frac{\bigtriangleup \Psi}{\Psi} 
\]

The laplacian of the wavefunction in the single-determinant
case is given as follows:

\[
\frac{\bigtriangleup \Psi(r)}{\Psi(r)} = \sum_{j=1}^{N_e} \bigtriangleup \Phi_j(r_i) D_{ji}^{-1}(r)
\]

*** Get
    
    #+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_kinetic_energy(qmckl_context context, double* const kinetic_energy);
    #+end_src

    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
qmckl_exit_code qmckl_get_kinetic_energy(qmckl_context context, double * const kinetic_energy) {
  
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return QMCKL_NULL_CONTEXT;
  }

  qmckl_exit_code rc;

  if(!qmckl_electron_provided(context)) return QMCKL_NOT_PROVIDED;

  if(!qmckl_nucleus_provided(context)) return QMCKL_NOT_PROVIDED;

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

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

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

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

  size_t sze = ctx->electron.walk_num * sizeof(double);
  memcpy(kinetic_energy, ctx->local_energy.e_kin, sze);

  return QMCKL_SUCCESS;
}
    #+end_src

*** Provide

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

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

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

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

  if(!(ctx->nucleus.provided)) {
    return qmckl_failwith( context,
                           QMCKL_NOT_PROVIDED,
                           "qmckl_electron",
                           NULL);
  }
  
  if(!(ctx->electron.provided)) {
    return qmckl_failwith( context,
                           QMCKL_NOT_PROVIDED,
                           "qmckl_electron",
                           NULL);
  }
  
  if (!ctx->ao_basis.provided) {
    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);
  }

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

 /* Compute if necessary */
  if (ctx->electron.coord_new_date > ctx->local_energy.e_kin_date) {

 /* Allocate array */
    if (ctx->local_energy.e_kin == NULL) {

      qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
      mem_info.size = ctx->electron.walk_num * sizeof(double);
      double* e_kin = (double*) qmckl_malloc(context, mem_info);

      if (e_kin == NULL) {
        return qmckl_failwith( context,
                               QMCKL_ALLOCATION_FAILED,
                               "qmckl_e_kin",
                               NULL);
      }
      ctx->local_energy.e_kin = e_kin;
    }

    qmckl_exit_code rc; 
    if (ctx->det.type == 'G') {
      rc = qmckl_compute_kinetic_energy(context,
                                 ctx->det.walk_num,
                                 ctx->det.det_num_alpha,
                                 ctx->det.det_num_beta,
                                 ctx->electron.up_num,
                                 ctx->electron.down_num,
                                 ctx->electron.num,
                                 ctx->det.mo_index_alpha,
                                 ctx->det.mo_index_beta,
                                 ctx->mo_basis.mo_num,
                                 ctx->mo_basis.mo_vgl,
                                 ctx->det.det_adj_matrix_alpha,
                                 ctx->det.det_adj_matrix_beta,
                                 ctx->local_energy.e_kin);
    } else {
      return qmckl_failwith( context,
                             QMCKL_FAILURE,
                             "compute_kinetic_energy",
                             "Not yet implemented");
    } 
    if (rc != QMCKL_SUCCESS) {
      return rc;
    }

    ctx->local_energy.e_kin_date = ctx->date;
  }

  return QMCKL_SUCCESS;
}
    #+end_src

*** Compute kinetic enregy
   :PROPERTIES:
   :Name:     qmckl_compute_kinetic_energy
   :CRetType: qmckl_exit_code
   :FRetType: qmckl_exit_code
   :END:

   #+NAME: qmckl_compute_kinetic_energy_args
   | ~qmckl_context~ | ~context~                                                             | in  | Global state                              |
   | ~int64_t~       | ~walk_num~                                                            | in  | Number of walkers                         |
   | ~int64_t~       | ~det_num_alpha~                                                       | in  | Number of determinants                    |
   | ~int64_t~       | ~det_num_beta~                                                        | in  | Number of determinants                    |
   | ~int64_t~       | ~alpha_num~                                                           | in  | Number of electrons                       |
   | ~int64_t~       | ~beta_num~                                                            | in  | Number of electrons                       |
   | ~int64_t~       | ~elec_num~                                                            | in  | Number of electrons                       |
   | ~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][walk_num][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~        | ~e_kin[walk_num]~                                                     | out | Kinetic energy                            |
    
    #+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_kinetic_energy_f(context, walk_num, &
     det_num_alpha, det_num_beta, alpha_num, beta_num, elec_num, mo_index_alpha, mo_index_beta, &
     mo_num, mo_vgl, det_adj_matrix_alpha, det_adj_matrix_beta, e_kin) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context)  , intent(in)  :: context
  integer*8, intent(in)             :: walk_num
  integer*8, intent(in)             :: det_num_alpha
  integer*8, intent(in)             :: det_num_beta
  integer*8, intent(in)             :: alpha_num
  integer*8, intent(in)             :: beta_num
  integer*8, intent(in)             :: elec_num
  integer*8, intent(in)             :: mo_num
  integer*8, intent(in)             :: mo_index_alpha(alpha_num, walk_num, det_num_alpha)
  integer*8, intent(in)             :: mo_index_beta(beta_num, walk_num, det_num_beta)
  double precision, intent(in)      :: mo_vgl(mo_num, elec_num, walk_num, 5)
  double precision, intent(in)      :: det_adj_matrix_alpha(alpha_num, alpha_num, walk_num, det_num_alpha)
  double precision, intent(in)      :: det_adj_matrix_beta(beta_num, beta_num, walk_num, det_num_beta)
  double precision, intent(inout)   :: e_kin(walk_num)
  integer*8 :: idet, iwalk, ielec, mo_id, imo

  info = QMCKL_SUCCESS

  if (context == QMCKL_NULL_CONTEXT) then
     info = QMCKL_INVALID_CONTEXT
     return
  endif

  if (walk_num <= 0) then
     info = QMCKL_INVALID_ARG_2
     return
  endif

  if (alpha_num <= 0) then
     info = QMCKL_INVALID_ARG_3
     return
  endif

  if (beta_num < 0) then
     info = QMCKL_INVALID_ARG_4
     return
  endif

  if (elec_num <= 0) then
     info = QMCKL_INVALID_ARG_5
     return
  endif

  e_kin = 0.0d0
  do idet = 1, det_num_alpha
  do iwalk = 1, walk_num
    ! Alpha part
    do imo = 1, alpha_num
    do ielec = 1, alpha_num
      mo_id = mo_index_alpha(ielec, iwalk, idet)
      e_kin(iwalk) = e_kin(iwalk) - 0.5d0 * det_adj_matrix_alpha(imo, ielec, iwalk, idet) * &
                                    mo_vgl(mo_id, ielec, iwalk, 5)
    end do
    end do
    ! Beta part
    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) * &
                                    mo_vgl(mo_id, alpha_num + ielec, iwalk, 5)
    end do
    end do
  end do
  end do

end function qmckl_compute_kinetic_energy_f
    #+end_src

   #+CALL: generate_c_header(table=qmckl_compute_kinetic_energy_args,rettyp=get_value("CRetType"),fname="qmckl_compute_kinetic_energy"))

   #+RESULTS:
   #+begin_src c :tangle (eval h_func) :comments org
   qmckl_exit_code qmckl_compute_kinetic_energy (
	 const qmckl_context context,
	 const int64_t walk_num,
	 const int64_t det_num_alpha,
	 const int64_t det_num_beta,
	 const int64_t alpha_num,
	 const int64_t beta_num,
	 const int64_t elec_num,
	 const int64_t* mo_index_alpha,
	 const int64_t* mo_index_beta,
	 const int64_t mo_num,
	 const double* mo_vgl,
	 const double* det_adj_matrix_alpha,
	 const double* det_adj_matrix_beta,
	 double* const e_kin ); 
   #+end_src

    #+CALL: generate_c_interface(table=qmckl_compute_kinetic_energy_args,rettyp=get_value("CRetType"),fname="qmckl_compute_kinetic_energy"))

    #+RESULTS:
    #+begin_src f90 :tangle (eval f) :comments org :exports none
    integer(c_int32_t) function qmckl_compute_kinetic_energy &
	(context, &
		 walk_num, &
		 det_num_alpha, &
		 det_num_beta, &
		 alpha_num, &
		 beta_num, &
		 elec_num, &
		 mo_index_alpha, &
		 mo_index_beta, &
		 mo_num, &
		 mo_vgl, &
		 det_adj_matrix_alpha, &
		 det_adj_matrix_beta, &
		 e_kin) &
	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 :: walk_num
      integer (c_int64_t) , intent(in)  , value :: det_num_alpha
      integer (c_int64_t) , intent(in)  , value :: det_num_beta
      integer (c_int64_t) , intent(in)  , value :: alpha_num
      integer (c_int64_t) , intent(in)  , value :: beta_num
      integer (c_int64_t) , intent(in)  , value :: elec_num
      integer (c_int64_t) , intent(in)          :: mo_index_alpha(alpha_num,walk_num,det_num_alpha)
      integer (c_int64_t) , intent(in)          :: mo_index_beta(beta_num,walk_num,det_num_beta)
      integer (c_int64_t) , intent(in)  , value :: mo_num
      real    (c_double ) , intent(in)          :: mo_vgl(mo_num,elec_num,walk_num,5)
      real    (c_double ) , intent(in)          :: det_adj_matrix_alpha(alpha_num,alpha_num,walk_num,det_num_alpha)
      real    (c_double ) , intent(in)          :: det_adj_matrix_beta(beta_num,beta_num,walk_num,det_num_beta)
      real    (c_double ) , intent(out)         :: e_kin(walk_num)

      integer(c_int32_t), external :: qmckl_compute_kinetic_energy_f
      info = qmckl_compute_kinetic_energy_f &
	     (context, &
		 walk_num, &
		 det_num_alpha, &
		 det_num_beta, &
		 alpha_num, &
		 beta_num, &
		 elec_num, &
		 mo_index_alpha, &
		 mo_index_beta, &
		 mo_num, &
		 mo_vgl, &
		 det_adj_matrix_alpha, &
		 det_adj_matrix_beta, &
		 e_kin)

    end function qmckl_compute_kinetic_energy
    #+end_src
    
*** Test
     #+begin_src c :tangle (eval c_test) :exports none

#define walk_num chbrclf_walk_num
#define elec_num chbrclf_elec_num
#define shell_num chbrclf_shell_num
#define ao_num chbrclf_ao_num

int64_t elec_up_num   = chbrclf_elec_up_num;
int64_t elec_dn_num   = chbrclf_elec_dn_num;
double* elec_coord    = &(chbrclf_elec_coord[0][0][0]);
const int64_t   nucl_num      = chbrclf_nucl_num;
const double*   nucl_charge   = chbrclf_charge;
const double*   nucl_coord    = &(chbrclf_nucl_coord[0][0]);

rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num);  
assert (rc == QMCKL_SUCCESS);

rc = qmckl_set_electron_walk_num (context, walk_num);
assert (rc == QMCKL_SUCCESS);

assert(qmckl_electron_provided(context));

rc = qmckl_set_electron_coord (context, 'N', elec_coord);                                     
assert(rc == QMCKL_SUCCESS);

rc = qmckl_set_nucleus_num (context, nucl_num);
assert(rc == QMCKL_SUCCESS);

rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]));
assert(rc == QMCKL_SUCCESS);

rc = qmckl_set_nucleus_charge(context, nucl_charge);
assert(rc == QMCKL_SUCCESS);

assert(qmckl_nucleus_provided(context));

const int64_t *  nucleus_index     =  &(chbrclf_basis_nucleus_index[0]);
const int64_t *  nucleus_shell_num =  &(chbrclf_basis_nucleus_shell_num[0]);
const int32_t *  shell_ang_mom     =  &(chbrclf_basis_shell_ang_mom[0]);
const int64_t *  shell_prim_num    =  &(chbrclf_basis_shell_prim_num[0]);
const int64_t *  shell_prim_index  =  &(chbrclf_basis_shell_prim_index[0]);
const double  *  shell_factor      =  &(chbrclf_basis_shell_factor[0]);
const double  *  exponent          =  &(chbrclf_basis_exponent[0]);
const double  *  coefficient       =  &(chbrclf_basis_coefficient[0]);
const double  *  prim_factor       =  &(chbrclf_basis_prim_factor[0]);
const double  *  ao_factor         =  &(chbrclf_basis_ao_factor[0]);

const char typ = 'G';

assert(!qmckl_ao_basis_provided(context));

rc = qmckl_set_ao_basis_type (context, typ);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));

rc = qmckl_set_ao_basis_shell_num (context, chbrclf_shell_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));

rc = qmckl_set_ao_basis_prim_num (context, chbrclf_prim_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));

rc = qmckl_set_ao_basis_nucleus_index (context, nucleus_index);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));

rc = qmckl_set_ao_basis_nucleus_shell_num (context, nucleus_shell_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));

rc = qmckl_set_ao_basis_shell_ang_mom (context, shell_ang_mom);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));

rc = qmckl_set_ao_basis_shell_factor  (context, shell_factor);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));

rc = qmckl_set_ao_basis_shell_prim_num (context, shell_prim_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));

rc = qmckl_set_ao_basis_shell_prim_index (context, shell_prim_index);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));

rc = qmckl_set_ao_basis_exponent      (context, exponent);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));

rc = qmckl_set_ao_basis_coefficient   (context, coefficient);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));

rc = qmckl_set_ao_basis_prim_factor (context, prim_factor);
assert(rc == QMCKL_SUCCESS);

rc = qmckl_set_ao_basis_ao_num(context, chbrclf_ao_num);
assert(rc == QMCKL_SUCCESS);

rc = qmckl_set_ao_basis_ao_factor (context, ao_factor);
assert(rc == QMCKL_SUCCESS);

assert(qmckl_ao_basis_provided(context));


double ao_vgl[5][walk_num][elec_num][chbrclf_ao_num];

rc = qmckl_get_ao_vgl(context, &(ao_vgl[0][0][0][0]));
assert (rc == QMCKL_SUCCESS);

/* Set up MO data */
const int64_t mo_num = chbrclf_mo_num;
rc = qmckl_set_mo_basis_mo_num(context, mo_num);
assert (rc == QMCKL_SUCCESS);

const double  * mo_coefficient          =  &(chbrclf_mo_coef[0]);

rc = qmckl_set_mo_basis_coefficient(context, mo_coefficient);
assert (rc == QMCKL_SUCCESS);

assert(qmckl_mo_basis_provided(context));

double mo_vgl[5][elec_num][chbrclf_mo_num];
rc = qmckl_get_mo_basis_vgl(context, &(mo_vgl[0][0][0]));
assert (rc == QMCKL_SUCCESS);

/* Set up determinant data */

const int64_t det_num_alpha = 1;
const int64_t det_num_beta  = 1;
int64_t mo_index_alpha[det_num_alpha][walk_num][elec_up_num];
int64_t mo_index_beta[det_num_alpha][walk_num][elec_dn_num];

int i, j, k;
for(k = 0; k < det_num_alpha; ++k)
  for(i = 0; i < walk_num; ++i)
    for(j = 0; j < elec_up_num; ++j)
      mo_index_alpha[k][i][j] = j + 1;
for(k = 0; k < det_num_beta; ++k)
  for(i = 0; i < walk_num; ++i)
    for(j = 0; j < elec_up_num; ++j)
      mo_index_beta[k][i][j] = j + 1;

rc = qmckl_set_determinant_type (context, typ);
assert(rc == QMCKL_SUCCESS);

rc = qmckl_set_determinant_walk_num (context, walk_num);
assert (rc == QMCKL_SUCCESS);

rc = qmckl_set_determinant_det_num_alpha (context, det_num_alpha);
assert (rc == QMCKL_SUCCESS);

rc = qmckl_set_determinant_det_num_beta (context, det_num_beta);
assert (rc == QMCKL_SUCCESS);

rc = qmckl_set_determinant_mo_index_alpha (context, &(mo_index_alpha[0][0][0]));
assert (rc == QMCKL_SUCCESS);

rc = qmckl_set_determinant_mo_index_beta (context, &(mo_index_beta[0][0][0]));
assert (rc == QMCKL_SUCCESS);

// Get alpha determinant

double det_vgl_alpha[det_num_alpha][walk_num][5][elec_up_num][elec_up_num];
double det_vgl_beta[det_num_beta][walk_num][5][elec_dn_num][elec_dn_num];

rc = qmckl_get_det_vgl_alpha(context, &(det_vgl_alpha[0][0][0][0][0]));
assert (rc == QMCKL_SUCCESS);

rc = qmckl_get_det_vgl_beta(context, &(det_vgl_beta[0][0][0][0][0]));
assert (rc == QMCKL_SUCCESS);

// Get adjoint of the slater-determinant

double det_inv_matrix_alpha[det_num_alpha][walk_num][elec_up_num][elec_up_num];
double det_inv_matrix_beta[det_num_beta][walk_num][elec_dn_num][elec_dn_num];

rc = qmckl_get_det_inv_matrix_alpha(context, &(det_inv_matrix_alpha[0][0][0][0]));
assert (rc == QMCKL_SUCCESS);

rc = qmckl_get_det_inv_matrix_beta(context, &(det_inv_matrix_beta[0][0][0][0]));
assert (rc == QMCKL_SUCCESS);

// Calculate the Kinetic energy

double kinetic_energy[walk_num];

rc = qmckl_get_kinetic_energy(context, &(kinetic_energy[0]));
assert (rc == QMCKL_SUCCESS);

     #+end_src
    
** Potential energy
   :PROPERTIES:
   :Name:     qmckl_compute_potential_energy
   :CRetType: qmckl_exit_code
   :FRetType: qmckl_exit_code
   :END:

The potential energy is the sum of all the following terms

\[
PE = \mathcal{V}_{ee} + \mathcal{V}_{en} + \mathcal{V}_{nn}
\]

The potential for is calculated as the sum of single electron
contributions.

\[
\mathcal{V}_{ee} = \sum_{i=1}^{N_e}\sum_{j<i} \frac{1}{r_{ij}}
\]

\[
\mathcal{V}_{en} = - \sum_{i=1}^{N_e}\sum_{A=1}^{N_n}\frac{Z_A}{r_{iA}}
\]

\[
\mathcal{V}_{nn} = \sum_{A=1}^{N_n}\sum_{B<A}\frac{Z_A Z_B}{r_{AB}}
\]

*** Get
    
    #+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_potential_energy(qmckl_context context, double* const potential_energy);
    #+end_src

    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
qmckl_exit_code qmckl_get_potential_energy(qmckl_context context, double * const potential_energy) {
  
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return QMCKL_NULL_CONTEXT;
  }

  qmckl_exit_code rc;

  if(!qmckl_electron_provided(context)) return QMCKL_NOT_PROVIDED;

  if(!qmckl_nucleus_provided(context)) return QMCKL_NOT_PROVIDED;

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

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

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

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

  size_t sze = ctx->electron.walk_num * sizeof(double);
  memcpy(potential_energy, ctx->local_energy.e_pot, sze);

  return QMCKL_SUCCESS;
}
    #+end_src

*** Provide

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

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

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

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

  qmckl_exit_code rc; 
  if(!(ctx->nucleus.provided)) {
    return qmckl_failwith( context,
                           QMCKL_NOT_PROVIDED,
                           "qmckl_electron",
                           NULL);
  }
  
  if(!(ctx->electron.provided)) {
    return qmckl_failwith( context,
                           QMCKL_NOT_PROVIDED,
                           "qmckl_electron",
                           NULL);
  }
  
  if (!ctx->ao_basis.provided) {
    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);
  }

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

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

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

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

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

 /* Compute if necessary */
  if (ctx->electron.coord_new_date > ctx->local_energy.e_pot_date) {

 /* Allocate array */
    if (ctx->local_energy.e_pot == NULL) {

      qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
      mem_info.size = ctx->electron.walk_num * sizeof(double);
      double* e_pot = (double*) qmckl_malloc(context, mem_info);

      if (e_pot == NULL) {
        return qmckl_failwith( context,
                               QMCKL_ALLOCATION_FAILED,
                               "qmckl_e_pot",
                               NULL);
      }
      ctx->local_energy.e_pot = e_pot;
    }

    if (ctx->det.type == 'G') {
      rc = qmckl_compute_potential_energy(context,
                                 ctx->det.walk_num,
                                 ctx->electron.num,
                                 ctx->nucleus.num,
                                 ctx->electron.ee_pot,
                                 ctx->electron.en_pot,
                                 ctx->nucleus.repulsion,
                                 ctx->local_energy.e_pot);
    } else {
      return qmckl_failwith( context,
                             QMCKL_FAILURE,
                             "compute_potential_energy",
                             "Not yet implemented");
    } 
    if (rc != QMCKL_SUCCESS) {
      return rc;
    }

    ctx->local_energy.e_pot_date = ctx->date;
  }

  return QMCKL_SUCCESS;
}
    #+end_src

*** Compute potential enregy
   :PROPERTIES:
   :Name:     qmckl_compute_potential_energy
   :CRetType: qmckl_exit_code
   :FRetType: qmckl_exit_code
   :END:

    #+NAME: qmckl_compute_potential_energy_args
   | ~qmckl_context~ | ~context~          | in  | Global state        |
   | ~int64_t~       | ~walk_num~         | in  | Number of walkers   |
   | ~int64_t~       | ~elec_num~         | in  | Number of electrons |
   | ~int64_t~       | ~nucl_num~         | in  | Number of MOs       |
   | ~double~        | ~ee_pot[walk_num]~ | in  | ee potential        |
   | ~double~        | ~en_pot[walk_num]~ | in  | en potential        |
   | ~double~        | ~repulsion~        | in  | en potential        |
   | ~double~        | ~e_pot[walk_num]~  | out | Potential energy    |
    
    #+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_potential_energy_f(context, walk_num, &
     elec_num, nucl_num, ee_pot, en_pot, repulsion, e_pot) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context)  , intent(in)  :: context
  integer*8, intent(in)             :: walk_num
  integer*8, intent(in)             :: elec_num
  integer*8, intent(in)             :: nucl_num
  double precision, intent(in)      :: ee_pot(walk_num)
  double precision, intent(in)      :: en_pot(walk_num)
  double precision, intent(in)      :: repulsion
  double precision, intent(inout)   :: e_pot(walk_num)
  integer*8 :: idet, iwalk, ielec, mo_id, imo

  info = QMCKL_SUCCESS

  if (context == QMCKL_NULL_CONTEXT) then
     info = QMCKL_INVALID_CONTEXT
     return
  endif

  if (walk_num <= 0) then
     info = QMCKL_INVALID_ARG_2
     return
  endif

  if (elec_num <= 0) then
     info = QMCKL_INVALID_ARG_3
     return
  endif

  e_pot = 0.0d0 + repulsion
  do iwalk = 1, walk_num
    e_pot(iwalk) = e_pot(iwalk) + ee_pot(iwalk) + en_pot(iwalk)
  end do

end function qmckl_compute_potential_energy_f
    #+end_src

    #+CALL: generate_c_header(table=qmckl_compute_potential_energy_args,rettyp=get_value("CRetType"),fname="qmckl_compute_potential_energy"))

   #+RESULTS:
   #+begin_src c :tangle (eval h_func) :comments org
   qmckl_exit_code qmckl_compute_potential_energy (
	 const qmckl_context context,
	 const int64_t walk_num,
	 const int64_t elec_num,
	 const int64_t nucl_num,
	 const double* ee_pot,
	 const double* en_pot,
	 const double repulsion,
	 double* const e_pot ); 
   #+end_src

    #+CALL: generate_c_interface(table=qmckl_compute_potential_energy_args,rettyp=get_value("CRetType"),fname="qmckl_compute_potential_energy"))

    #+RESULTS:
    #+begin_src f90 :tangle (eval f) :comments org :exports none
    integer(c_int32_t) function qmckl_compute_potential_energy &
	(context, walk_num, elec_num, nucl_num, ee_pot, en_pot, repulsion, e_pot) &
	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 :: walk_num
      integer (c_int64_t) , intent(in)  , value :: elec_num
      integer (c_int64_t) , intent(in)  , value :: nucl_num
      real    (c_double ) , intent(in)          :: ee_pot(walk_num)
      real    (c_double ) , intent(in)          :: en_pot(walk_num)
      real    (c_double ) , intent(in)  , value :: repulsion
      real    (c_double ) , intent(out)         :: e_pot(walk_num)

      integer(c_int32_t), external :: qmckl_compute_potential_energy_f
      info = qmckl_compute_potential_energy_f &
	     (context, walk_num, elec_num, nucl_num, ee_pot, en_pot, repulsion, e_pot)

    end function qmckl_compute_potential_energy
    #+end_src
    
*** Test
     #+begin_src c :tangle (eval c_test) :exports none
// Calculate the Potential energy

double potential_energy[walk_num];

rc = qmckl_get_potential_energy(context, &(potential_energy[0]));
assert (rc == QMCKL_SUCCESS);

     #+end_src

** Local energy
   :PROPERTIES:
   :Name:     qmckl_compute_local_energy
   :CRetType: qmckl_exit_code
   :FRetType: qmckl_exit_code
   :END:

The local energy is the sum of kinetic and potential energies.

\[
E_L = KE + PE
\]


*** Get
    
    #+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_local_energy(qmckl_context context, double* const local_energy);
    #+end_src

    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
qmckl_exit_code qmckl_get_local_energy(qmckl_context context, double * const local_energy) {
  
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return QMCKL_NULL_CONTEXT;
  }

  qmckl_exit_code rc;

  if(!qmckl_electron_provided(context)) return QMCKL_NOT_PROVIDED;

  if(!qmckl_nucleus_provided(context)) return QMCKL_NOT_PROVIDED;

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

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

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

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

  size_t sze = ctx->electron.walk_num * sizeof(double);
  memcpy(local_energy, ctx->local_energy.e_local, sze);

  return QMCKL_SUCCESS;
}
    #+end_src

*** Provide

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

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

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

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

  if(!(ctx->nucleus.provided)) {
    return qmckl_failwith( context,
                           QMCKL_NOT_PROVIDED,
                           "qmckl_electron",
                           NULL);
  }
  
  if(!(ctx->electron.provided)) {
    return qmckl_failwith( context,
                           QMCKL_NOT_PROVIDED,
                           "qmckl_electron",
                           NULL);
  }
  
  if (!ctx->ao_basis.provided) {
    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);
  }

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

  qmckl_exit_code rc;
  rc = qmckl_provide_kinetic_energy(context);
  if(rc != QMCKL_SUCCESS){
    return qmckl_failwith( context,
                           QMCKL_NOT_PROVIDED,
                           "qmckl_kinetic_energy",
                           NULL);
  }

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

 /* Compute if necessary */
  if (ctx->electron.coord_new_date > ctx->local_energy.e_local_date) {

 /* Allocate array */
    if (ctx->local_energy.e_local == NULL) {

      qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
      mem_info.size = ctx->electron.walk_num * sizeof(double);
      double* local_energy = (double*) qmckl_malloc(context, mem_info);

      if (local_energy == NULL) {
        return qmckl_failwith( context,
                               QMCKL_ALLOCATION_FAILED,
                               "qmckl_local_energy",
                               NULL);
      }
      ctx->local_energy.e_local = local_energy;
    }

    if (ctx->det.type == 'G') {
      rc = qmckl_compute_local_energy(context,
                                 ctx->det.walk_num,
                                 ctx->local_energy.e_kin,
                                 ctx->local_energy.e_pot,
                                 ctx->local_energy.e_local);
    } else {
      return qmckl_failwith( context,
                             QMCKL_FAILURE,
                             "compute_local_energy",
                             "Not yet implemented");
    } 
    if (rc != QMCKL_SUCCESS) {
      return rc;
    }

    ctx->local_energy.e_local_date = ctx->date;
  }

  return QMCKL_SUCCESS;
}
    #+end_src

*** Compute local enregy
   :PROPERTIES:
   :Name:     qmckl_compute_local_energy
   :CRetType: qmckl_exit_code
   :FRetType: qmckl_exit_code
   :END:

    #+NAME: qmckl_compute_local_energy_args
   | ~qmckl_context~ | ~context~           | in  | Global state      |
   | ~int64_t~       | ~walk_num~          | in  | Number of walkers |
   | ~double~        | ~e_kin[walk_num]~   | in  | e kinetic         |
   | ~double~        | ~e_pot[walk_num]~   | in  | e potential       |
   | ~double~        | ~e_local[walk_num]~ | out | local energy      |
    
    #+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_local_energy_f(context, walk_num, &
     e_kin, e_pot, e_local) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context)  , intent(in)  :: context
  integer*8, intent(in)             :: walk_num
  double precision, intent(in)      :: e_kin(walk_num)
  double precision, intent(in)      :: e_pot(walk_num)
  double precision, intent(inout)   :: e_local(walk_num)
  integer*8 :: idet, iwalk, ielec, mo_id, imo

  info = QMCKL_SUCCESS

  if (context == QMCKL_NULL_CONTEXT) then
     info = QMCKL_INVALID_CONTEXT
     return
  endif

  if (walk_num <= 0) then
     info = QMCKL_INVALID_ARG_2
     return
  endif

  do iwalk = 1, walk_num
    e_local(iwalk) = e_local(iwalk) + e_kin(iwalk) + e_pot(iwalk)
  end do

end function qmckl_compute_local_energy_f
    #+end_src

    #+CALL: generate_c_header(table=qmckl_compute_local_energy_args,rettyp=get_value("CRetType"),fname="qmckl_compute_local_energy"))

   #+RESULTS:
   #+begin_src c :tangle (eval h_func) :comments org
   qmckl_exit_code qmckl_compute_local_energy (
	 const qmckl_context context,
	 const int64_t walk_num,
	 const double* e_kin,
	 const double* e_pot,
	 double* const e_local ); 
   #+end_src

   #+CALL: generate_c_interface(table=qmckl_compute_local_energy_args,rettyp=get_value("CRetType"),fname="qmckl_compute_local_energy"))

    #+RESULTS:
    #+begin_src f90 :tangle (eval f) :comments org :exports none
    integer(c_int32_t) function qmckl_compute_local_energy &
	(context, walk_num, e_kin, e_pot, e_local) &
	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 :: walk_num
      real    (c_double ) , intent(in)          :: e_kin(walk_num)
      real    (c_double ) , intent(in)          :: e_pot(walk_num)
      real    (c_double ) , intent(out)         :: e_local(walk_num)

      integer(c_int32_t), external :: qmckl_compute_local_energy_f
      info = qmckl_compute_local_energy_f &
	     (context, walk_num, e_kin, e_pot, e_local)

    end function qmckl_compute_local_energy
    #+end_src
    
*** Test
     #+begin_src c :tangle (eval c_test) :exports none
// Calculate the Local energy

double local_energy[walk_num];

rc = qmckl_get_local_energy(context, &(local_energy[0]));
assert (rc == QMCKL_SUCCESS);

     #+end_src

** Drift vector
   :PROPERTIES:
   :Name:     qmckl_compute_drift_vector
   :CRetType: qmckl_exit_code
   :FRetType: qmckl_exit_code
   :END:

The local energy is the sum of kinetic and potential energies.

\[
E_L = KE + PE
\]


*** Get
    
    #+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_drift_vector(qmckl_context context, double* const drift_vector);
    #+end_src

    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
qmckl_exit_code qmckl_get_drift_vector(qmckl_context context, double * const drift_vector) {
  
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return QMCKL_NULL_CONTEXT;
  }

  qmckl_exit_code rc;

  if(!qmckl_electron_provided(context)) return QMCKL_NOT_PROVIDED;

  if(!qmckl_nucleus_provided(context)) return QMCKL_NOT_PROVIDED;

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

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

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

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

  size_t sze = ctx->electron.walk_num * 3 * sizeof(double);
  memcpy(drift_vector, ctx->local_energy.r_drift, sze);

  return QMCKL_SUCCESS;
}
    #+end_src

*** Provide

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

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

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

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

  if(!(ctx->nucleus.provided)) {
    return qmckl_failwith( context,
                           QMCKL_NOT_PROVIDED,
                           "qmckl_electron",
                           NULL);
  }
  
  if(!(ctx->electron.provided)) {
    return qmckl_failwith( context,
                           QMCKL_NOT_PROVIDED,
                           "qmckl_electron",
                           NULL);
  }
  
  if (!ctx->ao_basis.provided) {
    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);
  }

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

 /* Compute if necessary */
  if (ctx->electron.coord_new_date > ctx->local_energy.r_drift_date) {

 /* Allocate array */
    if (ctx->local_energy.r_drift == NULL) {

      qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
      mem_info.size = ctx->electron.walk_num * 3 * sizeof(double);
      double* r_drift = (double*) qmckl_malloc(context, mem_info);

      if (r_drift == NULL) {
        return qmckl_failwith( context,
                               QMCKL_ALLOCATION_FAILED,
                               "qmckl_r_drift",
                               NULL);
      }
      ctx->local_energy.r_drift = r_drift;
    }

    qmckl_exit_code rc; 
    if (ctx->det.type == 'G') {
      rc = qmckl_compute_drift_vector(context,
                                 ctx->det.walk_num,
                                 ctx->det.det_num_alpha,
                                 ctx->det.det_num_beta,
                                 ctx->electron.up_num,
                                 ctx->electron.down_num,
                                 ctx->electron.num,
                                 ctx->det.mo_index_alpha,
                                 ctx->det.mo_index_beta,
                                 ctx->mo_basis.mo_num,
                                 ctx->mo_basis.mo_vgl,
                                 ctx->det.det_adj_matrix_alpha,
                                 ctx->det.det_adj_matrix_beta,
                                 ctx->local_energy.r_drift);
    } else {
      return qmckl_failwith( context,
                             QMCKL_FAILURE,
                             "compute_drift_vector",
                             "Not yet implemented");
    } 
    if (rc != QMCKL_SUCCESS) {
      return rc;
    }

    ctx->local_energy.r_drift_date = ctx->date;
  }

  return QMCKL_SUCCESS;
}
    #+end_src

*** Compute drift vector
   :PROPERTIES:
   :Name:     qmckl_compute_drift_vector
   :CRetType: qmckl_exit_code
   :FRetType: qmckl_exit_code
   :END:

    #+NAME: qmckl_compute_drift_vector_args
   | ~qmckl_context~ | ~context~                                                             | in  | Global state                              |
   | ~int64_t~       | ~walk_num~                                                            | in  | Number of walkers                         |
   | ~int64_t~       | ~det_num_alpha~                                                       | in  | Number of determinants                    |
   | ~int64_t~       | ~det_num_beta~                                                        | in  | Number of determinants                    |
   | ~int64_t~       | ~alpha_num~                                                           | in  | Number of electrons                       |
   | ~int64_t~       | ~beta_num~                                                            | in  | Number of electrons                       |
   | ~int64_t~       | ~elec_num~                                                            | in  | Number of electrons                       |
   | ~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][walk_num][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                            |
    
    #+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_drift_vector_f(context, walk_num, &
     det_num_alpha, det_num_beta, alpha_num, beta_num, elec_num, mo_index_alpha, mo_index_beta, &
     mo_num, mo_vgl, det_adj_matrix_alpha, det_adj_matrix_beta, r_drift) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context)  , intent(in)  :: context
  integer*8, intent(in)             :: walk_num
  integer*8, intent(in)             :: det_num_alpha
  integer*8, intent(in)             :: det_num_beta
  integer*8, intent(in)             :: alpha_num
  integer*8, intent(in)             :: beta_num
  integer*8, intent(in)             :: elec_num
  integer*8, intent(in)             :: mo_num
  integer*8, intent(in)             :: mo_index_alpha(alpha_num, walk_num, det_num_alpha)
  integer*8, intent(in)             :: mo_index_beta(beta_num, walk_num, det_num_beta)
  double precision, intent(in)      :: mo_vgl(mo_num, elec_num, walk_num, 5)
  double precision, intent(in)      :: det_adj_matrix_alpha(alpha_num, alpha_num, walk_num, det_num_alpha)
  double precision, intent(in)      :: det_adj_matrix_beta(beta_num, beta_num, walk_num, det_num_beta)
  double precision, intent(inout)   :: r_drift(3,walk_num)
  integer*8 :: idet, iwalk, ielec, mo_id, imo

  info = QMCKL_SUCCESS

  if (context == QMCKL_NULL_CONTEXT) then
     info = QMCKL_INVALID_CONTEXT
     return
  endif

  if (walk_num <= 0) then
     info = QMCKL_INVALID_ARG_2
     return
  endif

  if (alpha_num <= 0) then
     info = QMCKL_INVALID_ARG_3
     return
  endif

  if (beta_num < 0) then
     info = QMCKL_INVALID_ARG_4
     return
  endif

  if (elec_num <= 0) then
     info = QMCKL_INVALID_ARG_5
     return
  endif

  r_drift = 0.0d0
  do idet = 1, det_num_alpha
  do iwalk = 1, walk_num
    ! Alpha part
    do imo = 1, alpha_num
    do ielec = 1, alpha_num
      mo_id = mo_index_alpha(ielec, iwalk, idet)
      r_drift(1,iwalk) = r_drift(1,iwalk) + 2.0d0 * det_adj_matrix_alpha(imo, ielec, iwalk, idet) * &
                                    mo_vgl(mo_id, ielec, iwalk, 2)
      r_drift(2,iwalk) = r_drift(2,iwalk) + 2.0d0 * det_adj_matrix_alpha(imo, ielec, iwalk, idet) * &
                                    mo_vgl(mo_id, ielec, iwalk, 3)
      r_drift(3,iwalk) = r_drift(3,iwalk) + 2.0d0 * det_adj_matrix_alpha(imo, ielec, iwalk, idet) * &
                                    mo_vgl(mo_id, ielec, iwalk, 4)
    end do
    end do
    ! Beta part
    do imo = 1, beta_num
    do ielec = 1, beta_num
      mo_id = mo_index_beta(ielec, iwalk, idet)
      r_drift(1,iwalk) = r_drift(1,iwalk) + 2.0d0 * det_adj_matrix_beta(imo, ielec, iwalk, idet) * &
                                    mo_vgl(mo_id, alpha_num + ielec, iwalk, 2)
      r_drift(2,iwalk) = r_drift(2,iwalk) + 2.0d0 * det_adj_matrix_beta(imo, ielec, iwalk, idet) * &
                                    mo_vgl(mo_id, alpha_num + ielec, iwalk, 3)
      r_drift(3,iwalk) = r_drift(3,iwalk) + 2.0d0 * det_adj_matrix_beta(imo, ielec, iwalk, idet) * &
                                    mo_vgl(mo_id, alpha_num + ielec, iwalk, 4)
    end do
    end do
  end do
  end do

end function qmckl_compute_drift_vector_f
    #+end_src

    #+CALL: generate_c_header(table=qmckl_compute_drift_vector_args,rettyp=get_value("CRetType"),fname="qmckl_compute_drift_vector"))

   #+RESULTS:
   #+begin_src c :tangle (eval h_func) :comments org
   qmckl_exit_code qmckl_compute_drift_vector (
	 const qmckl_context context,
	 const int64_t walk_num,
	 const int64_t det_num_alpha,
	 const int64_t det_num_beta,
	 const int64_t alpha_num,
	 const int64_t beta_num,
	 const int64_t elec_num,
	 const int64_t* mo_index_alpha,
	 const int64_t* mo_index_beta,
	 const int64_t mo_num,
	 const double* mo_vgl,
	 const double* det_adj_matrix_alpha,
	 const double* det_adj_matrix_beta,
	 double* const r_drift ); 
   #+end_src

   #+CALL: generate_c_interface(table=qmckl_compute_drift_vector_args,rettyp=get_value("CRetType"),fname="qmckl_compute_drift_vector"))

    #+RESULTS:
    #+begin_src f90 :tangle (eval f) :comments org :exports none
    integer(c_int32_t) function qmckl_compute_drift_vector &
	(context, &
		 walk_num, &
		 det_num_alpha, &
		 det_num_beta, &
		 alpha_num, &
		 beta_num, &
		 elec_num, &
		 mo_index_alpha, &
		 mo_index_beta, &
		 mo_num, &
		 mo_vgl, &
		 det_adj_matrix_alpha, &
		 det_adj_matrix_beta, &
		 r_drift) &
	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 :: walk_num
      integer (c_int64_t) , intent(in)  , value :: det_num_alpha
      integer (c_int64_t) , intent(in)  , value :: det_num_beta
      integer (c_int64_t) , intent(in)  , value :: alpha_num
      integer (c_int64_t) , intent(in)  , value :: beta_num
      integer (c_int64_t) , intent(in)  , value :: elec_num
      integer (c_int64_t) , intent(in)          :: mo_index_alpha(alpha_num,walk_num,det_num_alpha)
      integer (c_int64_t) , intent(in)          :: mo_index_beta(beta_num,walk_num,det_num_beta)
      integer (c_int64_t) , intent(in)  , value :: mo_num
      real    (c_double ) , intent(in)          :: mo_vgl(mo_num,elec_num,walk_num,5)
      real    (c_double ) , intent(in)          :: det_adj_matrix_alpha(alpha_num,alpha_num,walk_num,det_num_alpha)
      real    (c_double ) , intent(in)          :: det_adj_matrix_beta(beta_num,beta_num,walk_num,det_num_beta)
      real    (c_double ) , intent(out)         :: r_drift(3,walk_num)

      integer(c_int32_t), external :: qmckl_compute_drift_vector_f
      info = qmckl_compute_drift_vector_f &
	     (context, &
		 walk_num, &
		 det_num_alpha, &
		 det_num_beta, &
		 alpha_num, &
		 beta_num, &
		 elec_num, &
		 mo_index_alpha, &
		 mo_index_beta, &
		 mo_num, &
		 mo_vgl, &
		 det_adj_matrix_alpha, &
		 det_adj_matrix_beta, &
		 r_drift)

    end function qmckl_compute_drift_vector
    #+end_src
    
*** Test
     #+begin_src c :tangle (eval c_test) :exports none
// Calculate the Drift vector

double drift_vector[walk_num][3];

rc = qmckl_get_drift_vector(context, &(drift_vector[0][0]));
assert (rc == QMCKL_SUCCESS);

     #+end_src
* End of files                                                     :noexport:

  #+begin_src c :tangle (eval h_private_type)
#endif
  #+end_src

*** Test
  #+begin_src c :tangle (eval c_test)
    rc = qmckl_context_destroy(context);
    assert (rc == QMCKL_SUCCESS);

    return 0;
}
  #+end_src

*** Compute file names
    #+begin_src emacs-lisp
; The following is required to compute the file names

(setq pwd (file-name-directory buffer-file-name))
(setq name (file-name-nondirectory (substring buffer-file-name 0 -4)))
(setq f  (concat pwd name "_f.f90"))
(setq fh (concat pwd name "_fh.f90"))
(setq c  (concat pwd name ".c"))
(setq h  (concat name ".h"))
(setq h_private  (concat name "_private.h"))
(setq c_test  (concat pwd "test_" name ".c"))
(setq f_test  (concat pwd "test_" name "_f.f90"))

; Minted
(require 'ox-latex)
(setq org-latex-listings 'minted)
(add-to-list 'org-latex-packages-alist '("" "listings"))
(add-to-list 'org-latex-packages-alist '("" "color"))

    #+end_src


# -*- mode: org -*-
# vim: syntax=c