diff --git a/org/qmckl_local_energy.org b/org/qmckl_local_energy.org
new file mode 100644
index 0000000..b601321
--- /dev/null
+++ b/org/qmckl_local_energy.org
@@ -0,0 +1,535 @@
+#+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:
+
+*** 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_vgl_alpha,
+                                 ctx->det.det_vgl_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 alpha
+   :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_vgl_alpha[det_num_alpha][walk_num][5][alpha_num][alpha_num]~ | in  | Value, gradients and Laplacian of the Det |
+   | ~double~        | ~det_vgl_beta[det_num_beta][walk_num][5][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_vgl_alpha, det_vgl_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_vgl_alpha(alpha_num, alpha_num, 5, walk_num, det_num_alpha)
+  double precision, intent(in)      :: det_vgl_beta(beta_num, beta_num, 5, 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
+
+  idet = 1
+  e_kin = 0.0d0
+  do iwalk = 1, walk_num
+    ! Alpha part
+    do ielec = 1, alpha_num
+      mo_id = mo_index_beta(ielec, iwalk, idet)
+    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_vgl_alpha,
+	 const double* det_vgl_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_vgl_alpha, &
+		 det_vgl_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_vgl_alpha(alpha_num,alpha_num,5,walk_num,det_num_alpha)
+      real    (c_double ) , intent(in)          :: det_vgl_beta(beta_num,beta_num,5,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_vgl_alpha, &
+		 det_vgl_beta, &
+		 e_kin)
+
+    end function qmckl_compute_kinetic_energy
+    #+end_src
+    
+*** Test
+    
+* 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