diff --git a/QMC.org b/QMC.org
index 63cd513..aab007d 100644
--- a/QMC.org
+++ b/QMC.org
@@ -90,9 +90,14 @@
    :header-args:python: :tangle hydrogen.py
    :header-args:f90: :tangle hydrogen.f90
    :END:
-*** Write a function which computes the potential at $\mathbf{r}$
+
+*** Exercise 1
+
+    #+begin_exercise
+    Write a function which computes the potential at $\mathbf{r}$.
     The function accepts a 3-dimensional vector =r= as input arguments
     and returns the potential.
+    #+end_exercise
 
     $\mathbf{r}=\left( \begin{array}{c} x \\ y\\ z\end{array} \right)$, so
     $$
@@ -117,9 +122,12 @@ double precision function potential(r)
 end function potential
      #+END_SRC
 
-*** Write a function which computes the wave function at $\mathbf{r}$
+*** Exercise 2 
+    #+begin_exercise
+    Write a function which computes the wave function at $\mathbf{r}$.
     The function accepts a scalar =a= and a 3-dimensional vector =r= as
     input arguments, and returns a scalar.
+    #+end_exercise
 
     
     *Python*
@@ -137,9 +145,12 @@ double precision function psi(a, r)
 end function psi
      #+END_SRC
      
-*** Write a function which computes the local kinetic energy at $\mathbf{r}$
+*** Exercise 3
+    #+begin_exercise
+    Write a function which computes the local kinetic energy at $\mathbf{r}$.
     The function accepts =a= and =r= as input arguments and returns the
     local kinetic energy.
+    #+end_exercise
 
     The local kinetic energy is defined as $$-\frac{1}{2}\frac{\Delta \Psi}{\Psi}.$$
      
@@ -187,9 +198,12 @@ double precision function kinetic(a,r)
 end function kinetic
      #+END_SRC
 
-*** Write a function which computes the local energy at $\mathbf{r}$
+*** Exercise 4
+    #+begin_exercise
+    Write a function which computes the local energy at $\mathbf{r}$.
     The function accepts =x,y,z= as input arguments and returns the
     local energy.
+    #+end_exercise
    
     $$
     E_L(\mathbf{r}) = -\frac{1}{2} \frac{\Delta \Psi}{\Psi} (\mathbf{r}) + V(\mathbf{r})
@@ -218,11 +232,15 @@ end function e_loc
    :header-args:f90: :tangle plot_hydrogen.f90
    :END:
 
-   For multiple values of $a$ (0.1, 0.2, 0.5, 1., 1.5, 2.), plot the
-   local energy along the $x$ axis.
+   
+*** Exercise
+    #+begin_exercise
+    For multiple values of $a$ (0.1, 0.2, 0.5, 1., 1.5, 2.), plot the
+    local energy along the $x$ axis.
+    #+end_exercise
 
-   *Python*
-   #+BEGIN_SRC python :results none
+    *Python*
+    #+BEGIN_SRC python :results none
 import numpy as np
 import matplotlib.pyplot as plt
 
@@ -244,16 +262,16 @@ plt.tight_layout()
 plt.legend()
 
 plt.savefig("plot_py.png")
-   #+end_src
+    #+end_src
 
-   #+RESULTS:
+    #+RESULTS:
 
-   [[./plot_py.png]]
+    [[./plot_py.png]]
 
 
     
-   *Fortran*
-   #+begin_src f90 
+    *Fortran*
+    #+begin_src f90 
 program plot
   implicit none
   double precision, external :: e_loc
@@ -282,20 +300,20 @@ program plot
   end do
 
 end program plot
-   #+end_src
+    #+end_src
 
-   To compile and run:
+    To compile and run:
 
-   #+begin_src sh :exports both
+    #+begin_src sh :exports both
 gfortran hydrogen.f90 plot_hydrogen.f90 -o plot_hydrogen
 ./plot_hydrogen > data
-   #+end_src
+    #+end_src
 
-   #+RESULTS:
+    #+RESULTS:
 
-   To plot the data using gnuplot:
+    To plot the data using gnuplot:
 
-   #+begin_src gnuplot :file plot.png :exports both
+    #+begin_src gnuplot :file plot.png :exports both
 set grid
 set xrange [-5:5]
 set yrange [-2:1]
@@ -305,12 +323,12 @@ plot './data' index 0 using 1:2 with lines title 'a=0.1', \
      './data' index 3 using 1:2 with lines title 'a=1.0', \
      './data' index 4 using 1:2 with lines title 'a=1.5', \
      './data' index 5 using 1:2 with lines title 'a=2.0'
-   #+end_src
+    #+end_src
 
-   #+RESULTS:
-   [[file:plot.png]]
+    #+RESULTS:
+    [[file:plot.png]]
 
-** Compute numerically the expectation value of the energy
+** Numerical estimation of the energy
    :PROPERTIES:
    :header-args:python: :tangle energy_hydrogen.py
    :header-args:f90: :tangle energy_hydrogen.f90
@@ -327,18 +345,22 @@ plot './data' index 0 using 1:2 with lines title 'a=0.1', \
    w_i = \left[\Psi(\mathbf{r}_i)\right]^2 \delta \mathbf{r}
    $$
      
-   In this section, we will compute a numerical estimate of the
-   energy in a grid of $50\times50\times50$ points in the range
-   $(-5,-5,-5) \le \mathbf{r} \le (5,5,5)$.
-
    #+begin_note
    The energy is biased because:
    - The volume elements are not infinitely small (discretization error)
    - The energy is evaluated only inside the box (incompleteness of the space)
    #+end_note
 
-   *Python*
-     #+BEGIN_SRC python :results none
+   
+*** Exercise
+     #+begin_exercise
+    Compute a numerical estimate of the energy in a grid of
+    $50\times50\times50$ points in the range $(-5,-5,-5) \le
+    \mathbf{r} \le (5,5,5)$.
+     #+end_exercise
+
+    *Python*
+      #+BEGIN_SRC python :results none
 import numpy as np
 from hydrogen import e_loc, psi
 
@@ -363,19 +385,19 @@ for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]:
     E = E / norm
     print(f"a = {a} \t E = {E}")                
 
-     #+end_src
+      #+end_src
 
-     #+RESULTS:
-     : a = 0.1 	 E = -0.24518438948809218
-     : a = 0.2 	 E = -0.26966057967803525
-     : a = 0.5 	 E = -0.3856357612517407
-     : a = 0.9 	 E = -0.49435709786716214
-     : a = 1.0 	 E = -0.5
-     : a = 1.5 	 E = -0.39242967082602226
-     : a = 2.0 	 E = -0.08086980667844901
+      #+RESULTS:
+      : a = 0.1 	 E = -0.24518438948809218
+      : a = 0.2 	 E = -0.26966057967803525
+      : a = 0.5 	 E = -0.3856357612517407
+      : a = 0.9 	 E = -0.49435709786716214
+      : a = 1.0 	 E = -0.5
+      : a = 1.5 	 E = -0.39242967082602226
+      : a = 2.0 	 E = -0.08086980667844901
 
-     *Fortran*
-     #+begin_src f90 
+      *Fortran*
+      #+begin_src f90 
 program energy_hydrogen
   implicit none
   double precision, external :: e_loc, psi
@@ -414,24 +436,24 @@ program energy_hydrogen
   end do
 
 end program energy_hydrogen
-     #+end_src
+      #+end_src
 
-     To compile the Fortran and run it:
+      To compile the Fortran and run it:
 
-     #+begin_src sh :results output :exports both
+      #+begin_src sh :results output :exports both
 gfortran hydrogen.f90 energy_hydrogen.f90 -o energy_hydrogen
 ./energy_hydrogen 
-     #+end_src
+      #+end_src
 
-     #+RESULTS:
-     :  a =   0.10000000000000001          E =  -0.24518438948809140     
-     :  a =   0.20000000000000001          E =  -0.26966057967803236     
-     :  a =   0.50000000000000000          E =  -0.38563576125173815     
-     :  a =    1.0000000000000000          E =  -0.50000000000000000     
-     :  a =    1.5000000000000000          E =  -0.39242967082602065     
-     :  a =    2.0000000000000000          E =   -8.0869806678448772E-002
+      #+RESULTS:
+      :  a =   0.10000000000000001          E =  -0.24518438948809140     
+      :  a =   0.20000000000000001          E =  -0.26966057967803236     
+      :  a =   0.50000000000000000          E =  -0.38563576125173815     
+      :  a =    1.0000000000000000          E =  -0.50000000000000000     
+      :  a =    1.5000000000000000          E =  -0.39242967082602065     
+      :  a =    2.0000000000000000          E =   -8.0869806678448772E-002
 
-** Compute the variance of the local energy
+** Variance of the local energy
    :PROPERTIES:
    :header-args:python: :tangle variance_hydrogen.py
    :header-args:f90: :tangle variance_hydrogen.f90
@@ -450,8 +472,13 @@ gfortran hydrogen.f90 energy_hydrogen.f90 -o energy_hydrogen
    $\hat{H}$) the variance is zero, so the variance of the local
    energy can be used as a measure of the quality of a wave function.
 
+*** Exercise
+   #+begin_exercise
    Compute a numerical estimate of the variance of the local energy
-   in a grid of $50\times50\times50$ points in the range  $(-5,-5,-5) \le \mathbf{r} \le (5,5,5)$.
+   in a grid of $50\times50\times50$ points in the range
+   $(-5,-5,-5)
+   \le \mathbf{r} \le (5,5,5)$ for different values of $a$.
+   #+end_exercise
      
    *Python*
    #+begin_src python :results none
@@ -617,8 +644,11 @@ gfortran hydrogen.f90 variance_hydrogen.f90 -o variance_hydrogen
    E \pm \delta E, \text{ where } \delta E = \frac{\sigma}{\sqrt{M}}
    $$
    
+*** Exercise
+   #+begin_exercise
    Write a function returning the average and statistical error of an
    input array.
+   #+end_exercise
 
    *Python*
    #+BEGIN_SRC python :results none
@@ -656,26 +686,38 @@ end subroutine ave_error
    :header-args:f90: :tangle qmc_uniform.f90
    :END:
 
-   In this section we write a function to perform a Monte Carlo
-   calculation of the average energy.
+   We will now do our first Monte Carlo calculation to compute the
+   energy of the hydrogen atom.
+   
    At every Monte Carlo step:
 
-   - Draw 3 uniform random numbers in the interval $(-5,-5,-5) \le
+   - Draw a random point $\mathbf{r}_i$ in the box $(-5,-5,-5) \le
      (x,y,z) \le (5,5,5)$
-   - Compute $\Psi^2 \times E_L$ at this point and accumulate the
-     result in E
-   - Compute $\Psi^2$  at this point and accumulate the result in N
+   - Compute $[\Psi(\mathbf{r}_i)]^2$ and accumulate the result in a
+     variable =normalization=
+   - Compute $[\Psi(\mathbf{r}_i)]^2 \times E_L(\mathbf{r}_i)$, and accumulate the
+     result in a variable =energy=
 
-     Once all the steps have been computed, return the average energy
-     computed on the Monte Carlo calculation.
+   One Monte Carlo run will consist of $N$ Monte Carlo steps. Once all the
+   steps have been computed, the run returns the average energy
+   $\bar{E}_k$ over the $N$ steps of the run.
 
-     In the main program, write a loop to perform 30 Monte Carlo runs,
-     and compute the average energy and the associated statistical error.
+   To compute the statistical error, perform $M$ runs. The final
+   estimate of the energy will be the average over the $\bar{E}_k$,
+   and the variance of the $\bar{E}_k$ will be used to compute the
+   statistical error.
+   
+*** Exercise
 
-     Compute the energy of the wave function with $a=0.9$.
+    #+begin_exercise
+    Parameterize the wave function with $a=0.9$.  Perform 30
+    independent Monte Carlo runs, each with 100 000 Monte Carlo
+    steps. Store the final energies of each run and use this array to
+    compute the average energy and the associated error bar.
+    #+end_exercise
 
-   *Python*
-   #+BEGIN_SRC python :results output
+    *Python*
+    #+BEGIN_SRC python :results output
 from hydrogen  import *
 from qmc_stats import *
 
@@ -695,13 +737,13 @@ nmax = 100000
 X = [MonteCarlo(a,nmax) for i in range(30)]
 E, deltaE = ave_error(X)
 print(f"E = {E} +/- {deltaE}")
-   #+END_SRC
+    #+END_SRC
 
-   #+RESULTS:
-   : E = -0.4956255109300764 +/- 0.0007082875482711226
+    #+RESULTS:
+    : E = -0.4956255109300764 +/- 0.0007082875482711226
 
-   *Fortran*
-   #+BEGIN_SRC f90
+    *Fortran*
+    #+BEGIN_SRC f90
 subroutine uniform_montecarlo(a,nmax,energy)
   implicit none
   double precision, intent(in)  :: a
@@ -743,15 +785,15 @@ program qmc
   call ave_error(X,nruns,ave,err)
   print *, 'E = ', ave, '+/-', err
 end program qmc
-   #+END_SRC
+    #+END_SRC
 
-   #+begin_src sh :results output :exports both
+    #+begin_src sh :results output :exports both
 gfortran hydrogen.f90 qmc_stats.f90 qmc_uniform.f90 -o qmc_uniform
 ./qmc_uniform
-   #+end_src
+    #+end_src
 
-   #+RESULTS:
-   :  E =  -0.49588321986667677      +/-   7.1758863546737969E-004
+    #+RESULTS:
+    :  E =  -0.49588321986667677      +/-   7.1758863546737969E-004
 
 ** Gaussian sampling 
    :PROPERTIES:
@@ -773,6 +815,9 @@ gfortran hydrogen.f90 qmc_stats.f90 qmc_uniform.f90 -o qmc_uniform
    z_2 &=& \sqrt{-2 \ln u_1} \sin(2 \pi u_2) 
    \end{eqnarray*}
 
+   Here is a Fortran implementation returning a Gaussian-distributed
+   n-dimensional vector $\mathbf{z}$;
+
    *Fortran*
    #+BEGIN_SRC f90 :tangle qmc_stats.f90
 subroutine random_gauss(z,n)
@@ -805,14 +850,16 @@ end subroutine random_gauss
    #+END_SRC
 
 
-   Now the equation for the energy is changed into
+   Now the sampling probability can be inserted into the equation of the energy:
    
    \[
    E = \frac{\int P(\mathbf{r}) \frac{\left[\Psi(\mathbf{r})\right]^2}{P(\mathbf{r})}\, \frac{\hat{H} \Psi(\mathbf{r})}{\Psi(\mathbf{r})}\,d\mathbf{r}}{\int P(\mathbf{r}) \frac{\left[\Psi(\mathbf{r}) \right]^2}{P(\mathbf{r})} d\mathbf{r}}
    \]
-   with
+
+   with the Gaussian probability
+
    \[
-   P(\mathbf{r}) = \frac{1}{(2 \pi)^{3/2}}\exp\left( -\frac{\mathbf{r}^2}{2} \right)
+   P(\mathbf{r}) = \frac{1}{(2 \pi)^{3/2}}\exp\left( -\frac{\mathbf{r}^2}{2} \right).
    \]
 
    As the coordinates are drawn with probability $P(\mathbf{r})$, the
@@ -823,8 +870,17 @@ end subroutine random_gauss
    w_i = \frac{\left[\Psi(\mathbf{r}_i)\right]^2}{P(\mathbf{r}_i)} \delta \mathbf{r}
    $$
 
-   *Python*
-   #+BEGIN_SRC python :results output
+   
+*** Exercise
+
+    #+begin_exercise
+    Modify the exercise of the previous section to sample with
+    Gaussian-distributed random numbers. Can you see an reduction in
+    the statistical error?
+    #+end_exercise
+
+    *Python*
+    #+BEGIN_SRC python :results output
 from hydrogen  import *
 from qmc_stats import *
 
@@ -848,14 +904,14 @@ nmax = 100000
 X = [MonteCarlo(a,nmax) for i in range(30)]
 E, deltaE = ave_error(X)
 print(f"E = {E} +/- {deltaE}")
-   #+END_SRC
+    #+END_SRC
 
-   #+RESULTS:
-   : E = -0.49507506093129827 +/- 0.00014164037765553668
+    #+RESULTS:
+    : E = -0.49507506093129827 +/- 0.00014164037765553668
 
    
-   *Fortran*
-   #+BEGIN_SRC f90
+    *Fortran*
+    #+BEGIN_SRC f90
 double precision function gaussian(r)
   implicit none
   double precision, intent(in) :: r(3)
@@ -904,15 +960,15 @@ program qmc
   call ave_error(X,nruns,ave,err)
   print *, 'E = ', ave, '+/-', err
 end program qmc
-   #+END_SRC
+    #+END_SRC
 
-   #+begin_src sh :results output :exports both
+    #+begin_src sh :results output :exports both
 gfortran hydrogen.f90 qmc_stats.f90 qmc_gaussian.f90 -o qmc_gaussian
 ./qmc_gaussian
-   #+end_src
+    #+end_src
 
-   #+RESULTS:
-   :  E =  -0.49606057056767766      +/-   1.3918807547836872E-004
+    #+RESULTS:
+    :  E =  -0.49606057056767766      +/-   1.3918807547836872E-004
 
 ** Sampling with $\Psi^2$
    :PROPERTIES:
@@ -933,88 +989,105 @@ gfortran hydrogen.f90 qmc_stats.f90 qmc_gaussian.f90 -o qmc_gaussian
    E \approx \frac{1}{M}\sum_{i=1}^M E_L(\mathbf{r}_i)
    $$
    
-   To generate the probability density $\Psi^2$, we consider a
-   diffusion process characterized by a time-dependent density
-   $[\Psi(\mathbf{r},t)]^2$, which obeys the Fokker-Planck equation:
-
-   \[
-   \frac{\partial \Psi^2}{\partial t} = \sum_i D
-   \frac{\partial}{\partial \mathbf{r}_i} \left(
-   \frac{\partial}{\partial \mathbf{r}_i} - F_i(\mathbf{r}) \right)
-   [\Psi(\mathbf{r},t)]^2.
-   \]
    
-   $D$ is the diffusion constant and $F_i$ is the i-th component of a
-   drift velocity caused by an external potential. For a stationary
-   density, \( \frac{\partial \Psi^2}{\partial t} = 0 \), so
+*** Importance sampling
 
-   \begin{eqnarray*}
-   0 & = & \sum_i D
-   \frac{\partial}{\partial \mathbf{r}_i} \left(
-   \frac{\partial}{\partial \mathbf{r}_i} - F_i(\mathbf{r}) \right)
-   [\Psi(\mathbf{r})]^2 \\
-   0 & = & \sum_i D
-   \frac{\partial}{\partial \mathbf{r}_i} \left(
-   \frac{\partial [\Psi(\mathbf{r})]^2}{\partial \mathbf{r}_i} -
-   F_i(\mathbf{r})\,[\Psi(\mathbf{r})]^2 \right) \\
-   0 & = &
-   \frac{\partial^2 \Psi^2}{\partial \mathbf{r}_i^2} -
-   \frac{\partial   F_i   }{\partial \mathbf{r}_i}[\Psi(\mathbf{r})]^2  - 
-   \frac{\partial   \Psi^2}{\partial \mathbf{r}_i} F_i(\mathbf{r})
-   \end{eqnarray*}
+    To generate the probability density $\Psi^2$, we consider a
+    diffusion process characterized by a time-dependent density
+    $[\Psi(\mathbf{r},t)]^2$, which obeys the Fokker-Planck equation:
 
-   we search for a drift function which satisfies 
-
-   \[
-   \frac{\partial^2 \Psi^2}{\partial \mathbf{r}_i^2} =
-   [\Psi(\mathbf{r})]^2 \frac{\partial   F_i   }{\partial \mathbf{r}_i} + 
-   \frac{\partial   \Psi^2}{\partial \mathbf{r}_i} F_i(\mathbf{r})
-   \]
-
-   to obtain a second derivative on the left, we need the drift to be
-   of the form
-   \[
-   F_i(\mathbf{r}) = g(\mathbf{r}) \frac{\partial \Psi^2}{\partial \mathbf{r}_i}
-   \]
-
-   \[
-   \frac{\partial^2 \Psi^2}{\partial \mathbf{r}_i^2} =
-   [\Psi(\mathbf{r})]^2 \frac{\partial
-   g(\mathbf{r})}{\partial \mathbf{r}_i}\frac{\partial \Psi^2}{\partial \mathbf{r}_i} + 
-   [\Psi(\mathbf{r})]^2 g(\mathbf{r}) \frac{\partial^2
-   \Psi^2}{\partial \mathbf{r}_i^2} + 
-   \frac{\partial   \Psi^2}{\partial \mathbf{r}_i} 
-   g(\mathbf{r}) \frac{\partial \Psi^2}{\partial \mathbf{r}_i}
-   \]
+    \[
+    \frac{\partial \Psi^2}{\partial t} = \sum_i D
+    \frac{\partial}{\partial \mathbf{r}_i} \left(
+    \frac{\partial}{\partial \mathbf{r}_i} - F_i(\mathbf{r}) \right)
+    [\Psi(\mathbf{r},t)]^2.
+    \]
    
-   $g = 1 / \Psi^2$ satisfies this equation, so 
+    $D$ is the diffusion constant and $F_i$ is the i-th component of a
+    drift velocity caused by an external potential. For a stationary
+    density, \( \frac{\partial \Psi^2}{\partial t} = 0 \), so
 
-   \[
-   F(\mathbf{r}) = \frac{\nabla [\Psi(\mathbf{r})]^2}{[\Psi(\mathbf{r})]^2} = 2 \frac{\nabla
-   \Psi(\mathbf{r})}{\Psi(\mathbf{r})} = 2 \nabla \left( \log \Psi(\mathbf{r}) \right)
-   \]
+    \begin{eqnarray*}
+    0 & = & \sum_i D
+    \frac{\partial}{\partial \mathbf{r}_i} \left(
+    \frac{\partial}{\partial \mathbf{r}_i} - F_i(\mathbf{r}) \right)
+    [\Psi(\mathbf{r})]^2 \\
+    0 & = & \sum_i D
+    \frac{\partial}{\partial \mathbf{r}_i} \left(
+    \frac{\partial [\Psi(\mathbf{r})]^2}{\partial \mathbf{r}_i} -
+    F_i(\mathbf{r})\,[\Psi(\mathbf{r})]^2 \right) \\
+    0 & = &
+    \frac{\partial^2 \Psi^2}{\partial \mathbf{r}_i^2} -
+    \frac{\partial   F_i   }{\partial \mathbf{r}_i}[\Psi(\mathbf{r})]^2  - 
+    \frac{\partial   \Psi^2}{\partial \mathbf{r}_i} F_i(\mathbf{r})
+    \end{eqnarray*}
 
-   following drifted diffusion scheme:
+    we search for a drift function which satisfies 
 
-   \[
-   \mathbf{r}_{n+1} = \mathbf{r}_{n} + \tau \frac{\nabla
-   \Psi(\mathbf{r})}{\Psi(\mathbf{r})} + \eta \sqrt{\tau} 
-   \]
+    \[
+    \frac{\partial^2 \Psi^2}{\partial \mathbf{r}_i^2} =
+    [\Psi(\mathbf{r})]^2 \frac{\partial   F_i   }{\partial \mathbf{r}_i} + 
+    \frac{\partial   \Psi^2}{\partial \mathbf{r}_i} F_i(\mathbf{r})
+    \]
 
-   where $\eta$ is a normally-distributed Gaussian random number.
+    to obtain a second derivative on the left, we need the drift to be
+    of the form
+    \[
+    F_i(\mathbf{r}) = g(\mathbf{r}) \frac{\partial \Psi^2}{\partial \mathbf{r}_i}
+    \]
+
+    \[
+    \frac{\partial^2 \Psi^2}{\partial \mathbf{r}_i^2} =
+    [\Psi(\mathbf{r})]^2 \frac{\partial
+    g(\mathbf{r})}{\partial \mathbf{r}_i}\frac{\partial \Psi^2}{\partial \mathbf{r}_i} + 
+    [\Psi(\mathbf{r})]^2 g(\mathbf{r}) \frac{\partial^2
+    \Psi^2}{\partial \mathbf{r}_i^2} + 
+    \frac{\partial   \Psi^2}{\partial \mathbf{r}_i} 
+    g(\mathbf{r}) \frac{\partial \Psi^2}{\partial \mathbf{r}_i}
+    \]
    
+    $g = 1 / \Psi^2$ satisfies this equation, so 
 
-   First, write a function to compute the drift vector $\frac{\nabla \Psi(\mathbf{r})}{\Psi(\mathbf{r})}$.
+    \[
+    F(\mathbf{r}) = \frac{\nabla [\Psi(\mathbf{r})]^2}{[\Psi(\mathbf{r})]^2} = 2 \frac{\nabla
+    \Psi(\mathbf{r})}{\Psi(\mathbf{r})} = 2 \nabla \left( \log \Psi(\mathbf{r}) \right)
+    \]
+
+    In statistical mechanics, Fokker-Planck trajectories are generated
+    by a Langevin equation:
+
+    \[
+     \frac{\partial \mathbf{r}(t)}{\partial t} = 2D \frac{\nabla
+     \Psi(\mathbf{r}(t))}{\Psi} + \eta
+    \]
+
+    where $\eta$ is a normally-distributed fluctuating random force.
+
+    Discretizing this differential equation gives the following drifted
+    diffusion scheme:
+
+    \[
+    \mathbf{r}_{n+1} = \mathbf{r}_{n} + \tau\, 2D \frac{\nabla
+    \Psi(\mathbf{r})}{\Psi(\mathbf{r})} + \chi 
+    \]
+    where $\chi$ is a Gaussian random variable with zero mean and
+    variance $\tau\,2D$.
    
-   *Python*
-   #+BEGIN_SRC python
+**** Exercise 1
+
+      #+begin_exercise
+      Write a function to compute the drift vector $\frac{\nabla \Psi(\mathbf{r})}{\Psi(\mathbf{r})}$.
+      #+end_exercise
+   
+      *Python*
+      #+BEGIN_SRC python :tangle hydrogen.py
 def drift(a,r):
   ar_inv = -a/np.sqrt(np.dot(r,r))
   return r * ar_inv
-   #+END_SRC
+      #+END_SRC
 
-   *Fortran*
-   #+BEGIN_SRC f90
+      *Fortran*
+      #+BEGIN_SRC f90 :tangle hydrogen.f90
 subroutine drift(a,r,b)
   implicit none
   double precision, intent(in)  :: a, r(3)
@@ -1023,54 +1096,50 @@ subroutine drift(a,r,b)
   ar_inv = -a / dsqrt(r(1)*r(1) + r(2)*r(2) + r(3)*r(3))
   b(:) = r(:) * ar_inv
 end subroutine drift
-   #+END_SRC
+      #+END_SRC
 
+**** TODO Exercise 2
 
-   Now we can write the Monte Carlo sampling:
+     #+begin_exercise
+     Sample $\Psi^2$ approximately using the drifted diffusion scheme,
+     with a diffusion constant $D=1/2$. You can use a time step of
+     0.001 a.u.
+     #+end_exercise
    
-   *Python*
-   #+BEGIN_SRC python 
+      *Python*
+      #+BEGIN_SRC python :results output :tangle vmc.py 
+from hydrogen  import *
+from qmc_stats import *
 def MonteCarlo(a,tau,nmax):
     E = 0.
     N = 0.
-    sq_tau = sqrt(tau)
+    sq_tau = np.sqrt(tau)
     r_old = np.random.normal(loc=0., scale=1.0, size=(3))
     d_old = drift(a,r_old)
     d2_old = np.dot(d_old,d_old)
     psi_old = psi(a,r_old)
     for istep in range(nmax):
-        eta = np.random.normal(loc=0., scale=1.0, size=(3))
-        r_new = r_old + tau * d_old + sq_tau * eta
-        d_new = drift(a,r_new)
-        d2_new = np.dot(d_new,d_new)
-        psi_new = psi(a,r_new)
-        # Metropolis
-        prod = np.dot((d_new + d_old), (r_new - r_old))
-        argexpo = 0.5 * (d2_new - d2_old)*tau + prod
-        q = psi_new / psi_old
-        q = np.exp(-argexpo) * q*q
-        if np.random.uniform() < q:
-            r_old = r_new
-            d_old = d_new
-            d2_old = d2_new
-            psi_old = psi_new
-            N += 1.
-            E += e_loc(a,r_old)
+        chi = np.random.normal(loc=0., scale=1.0, size=(3))
+        r_new = r_old + tau * d_old + sq_tau * chi
+        N += 1.
+        E += e_loc(a,r_new)
+        r_old = r_new
     return E/N
 
 
+a = 0.9
 nmax = 100000
-tau = 0.1
+tau = 0.001
 X = [MonteCarlo(a,tau,nmax) for i in range(30)]
 E, deltaE = ave_error(X)
 print(f"E = {E} +/- {deltaE}")
-   #+END_SRC
+      #+END_SRC
 
-   #+RESULTS:
-   : E = -0.4951783346213532 +/- 0.00022067316984271938
+      #+RESULTS:
+      : E = -0.4112049153828464 +/- 0.00027934927432953063
    
-   *Fortran*
-   #+BEGIN_SRC f90
+      *Fortran*
+      #+BEGIN_SRC f90
 subroutine variational_montecarlo(a,nmax,energy)
   implicit none
   double precision, intent(in)  :: a
@@ -1111,15 +1180,153 @@ program qmc
   call ave_error(X,nruns,ave,err)
   print *, 'E = ', ave, '+/-', err
 end program qmc
-   #+END_SRC
+      #+END_SRC
 
-   #+begin_src sh :results output :exports both
+      #+begin_src sh :results output :exports both
 gfortran hydrogen.f90 qmc_stats.f90 vmc.f90 -o vmc
 ./vmc
-   #+end_src
+      #+end_src
+     
+*** Metropolis algorithm
 
+    Discretizing the differential equation to generate the desired
+    probability density will suffer from a discretization error
+    leading to biases in the averages. The [[https://en.wikipedia.org/wiki/Metropolis%E2%80%93Hastings_algorithm][Metropolis-Hastings
+    sampling algorithm]] removes exactly the discretization errors, so
+    large time steps can be employed.
+
+    After the new position $\mathbf{r}_{n+1}$ has been computed, an
+    additional accept/reject step is performed. The acceptance
+    probability $A$ is chosen so that it is consistent with the
+    probability of leaving $\mathbf{r}_n$ and the probability of
+    entering $\mathbf{r}_{n+1}$:
+
+    \[ A(\mathbf{r}_{n} \rightarrow \mathbf{r}_{n+1}) = \min \left( 1,
+    \frac{g(\mathbf{r}_{n+1} \rightarrow \mathbf{r}_{n}) P(\mathbf{r}_{n+1})}
+    {g(\mathbf{r}_{n} \rightarrow \mathbf{r}_{n+1}) P(\mathbf{r}_{n})}
+    \right)
+    \]
+
+    In our Hydrogen atom example, $P$ is $\Psi^2$ and $g$ is a
+    solution of the discretized Fokker-Planck equation:
+    
+    \begin{eqnarray*}
+    P(r_{n}) &=& \Psi^2(\mathbf{r}_n) \\
+    g(\mathbf{r}_{n} \rightarrow \mathbf{r}_{n+1}) & = &
+    \frac{1}{(4\pi\,D\,\tau)^{3/2}} \exp \left[ - \frac{\left(
+    \mathbf{r}_{n+1} - \mathbf{r}_{n} - 2D \frac{\nabla
+    \Psi(\mathbf{r}_n)}{\Psi(\mathbf{r}_n)} \right)^2}{4D\,\tau} \right]
+    \end{eqnarray*}
+    
+    The accept/reject step is the following:
+    - Compute $A(\mathbf{r}_{n} \rightarrow \mathbf{r}_{n+1})$.
+    - Draw a uniform random number $u$
+    - if $u \le A(\mathbf{r}_{n} \rightarrow \mathbf{r}_{n+1})$, accept
+      the move
+    - if $u>A(\mathbf{r}_{n} \rightarrow \mathbf{r}_{n+1})$, reject
+      the move: set $\mathbf{r}_{n+1} = \mathbf{r}_{n}$, but *don't
+      remove the sample from the average!*
+
+
+**** TODO Exercise 
+
+     #+begin_exercise
+     Modify the previous program to introduce the accept/reject step.
+     You should recover the unbiased result.
+     #+end_exercise
+   
+      *Python*
+      #+BEGIN_SRC python 
+def MonteCarlo(a,tau,nmax):
+    E = 0.
+    N = 0.
+    sq_tau = np.sqrt(tau)
+    r_old = np.random.normal(loc=0., scale=1.0, size=(3))
+    d_old = drift(a,r_old)
+    d2_old = np.dot(d_old,d_old)
+    psi_old = psi(a,r_old)
+    for istep in range(nmax):
+        eta = np.random.normal(loc=0., scale=1.0, size=(3))
+        r_new = r_old + tau * d_old + sq_tau * eta
+        d_new = drift(a,r_new)
+        d2_new = np.dot(d_new,d_new)
+        psi_new = psi(a,r_new)
+        # Metropolis
+        prod = np.dot((d_new + d_old), (r_new - r_old))
+        argexpo = 0.5 * (d2_new - d2_old)*tau + prod
+        q = psi_new / psi_old
+        q = np.exp(-argexpo) * q*q
+        if np.random.uniform() < q:
+            r_old = r_new
+            d_old = d_new
+            d2_old = d2_new
+            psi_old = psi_new
+            N += 1.
+            E += e_loc(a,r_old)
+    return E/N
+
+
+nmax = 100000
+tau = 0.1
+X = [MonteCarlo(a,tau,nmax) for i in range(30)]
+E, deltaE = ave_error(X)
+print(f"E = {E} +/- {deltaE}")
+      #+END_SRC
+
+      #+RESULTS:
+      : E = -0.4951783346213532 +/- 0.00022067316984271938
+   
+      *Fortran*
+      #+BEGIN_SRC f90
+subroutine variational_montecarlo(a,nmax,energy)
+  implicit none
+  double precision, intent(in)  :: a
+  integer         , intent(in)  :: nmax 
+  double precision, intent(out) :: energy
+
+  integer*8 :: istep
+
+  double precision :: norm, r(3), w
+
+  double precision, external :: e_loc, psi, gaussian
+
+  energy = 0.d0
+  norm   = 0.d0
+  do istep = 1,nmax
+     call random_gauss(r,3)
+     w = psi(a,r) 
+     w = w*w / gaussian(r)
+     norm = norm + w
+     energy = energy + w * e_loc(a,r)
+  end do
+  energy = energy / norm
+end subroutine variational_montecarlo
+
+program qmc
+  implicit none
+  double precision, parameter :: a = 0.9
+  integer         , parameter :: nmax = 100000
+  integer         , parameter :: nruns = 30
+
+  integer :: irun
+  double precision :: X(nruns)
+  double precision :: ave, err
+
+  do irun=1,nruns
+     call gaussian_montecarlo(a,nmax,X(irun))
+  enddo
+  call ave_error(X,nruns,ave,err)
+  print *, 'E = ', ave, '+/-', err
+end program qmc
+      #+END_SRC
+
+      #+begin_src sh :results output :exports both
+gfortran hydrogen.f90 qmc_stats.f90 vmc.f90 -o vmc
+./vmc
+      #+end_src
+     
   
-* Diffusion Monte Carlo
+* TODO Diffusion Monte Carlo
 
   We will now consider the H_2 molecule in a minimal basis composed of the
   $1s$ orbitals of the hydrogen atoms:
@@ -1133,3 +1340,4 @@ gfortran hydrogen.f90 qmc_stats.f90 vmc.f90 -o vmc
   the nuclei.
 
 
+  
diff --git a/energy_hydrogen.f90 b/energy_hydrogen.f90
index 8b3fe7f..db75ab0 100644
--- a/energy_hydrogen.f90
+++ b/energy_hydrogen.f90
@@ -26,7 +26,6 @@ program energy_hydrogen
               r(3) = x(l)
               w = psi(a(j),r)
               w = w * w * delta
-
               energy = energy + w * e_loc(a(j), r)
               norm   = norm   + w 
            end do
diff --git a/energy_hydrogen.py b/energy_hydrogen.py
index 719ba7b..969ed4d 100644
--- a/energy_hydrogen.py
+++ b/energy_hydrogen.py
@@ -1,23 +1,23 @@
 import numpy as np
 from hydrogen import e_loc, psi
 
- interval = np.linspace(-5,5,num=50)
- delta = (interval[1]-interval[0])**3
+interval = np.linspace(-5,5,num=50)
+delta = (interval[1]-interval[0])**3
 
- r = np.array([0.,0.,0.])
+r = np.array([0.,0.,0.])
 
- for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]:
-     E = 0.
-     norm = 0.
-       for x in interval:
-           r[0] = x
-             for y in interval:
-                 r[1] = y
-                   for z in interval:
-                       r[2] = z
-                       w = psi(a,r)
-                       w = w * w * delta
-                       E    += w * e_loc(a,r)
-                       norm += w 
-                       E = E / norm
-                       print(f"a = {a} \t E = {E}")
+for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]:
+    E = 0.
+    norm = 0.
+      for x in interval:
+          r[0] = x
+            for y in interval:
+                r[1] = y
+                  for z in interval:
+                      r[2] = z
+                      w = psi(a,r)
+                      w = w * w * delta
+                      E    += w * e_loc(a,r)
+                      norm += w 
+    E = E / norm
+    print(f"a = {a} \t E = {E}")
diff --git a/hydrogen.f90 b/hydrogen.f90
index eccca10..b074930 100644
--- a/hydrogen.f90
+++ b/hydrogen.f90
@@ -23,3 +23,12 @@ double precision function e_loc(a,r)
   double precision, external   :: kinetic, potential
   e_loc = kinetic(a,r) + potential(r)
 end function e_loc
+
+subroutine drift(a,r,b)
+  implicit none
+  double precision, intent(in)  :: a, r(3)
+  double precision, intent(out) :: b(3)
+  double precision :: ar_inv
+  ar_inv = -a / dsqrt(r(1)*r(1) + r(2)*r(2) + r(3)*r(3))
+  b(:) = r(:) * ar_inv
+end subroutine drift
diff --git a/hydrogen.py b/hydrogen.py
index ff8fa96..7a78652 100644
--- a/hydrogen.py
+++ b/hydrogen.py
@@ -11,3 +11,7 @@ def kinetic(a,r):
 
 def e_loc(a,r):
     return kinetic(a,r) + potential(r)
+
+def drift(a,r):
+  ar_inv = -a/np.sqrt(np.dot(r,r))
+  return r * ar_inv
diff --git a/qmc_uniform.py b/qmc_uniform.py
index 616e72d..36a41ec 100644
--- a/qmc_uniform.py
+++ b/qmc_uniform.py
@@ -1,19 +1,19 @@
 from hydrogen  import *
 from qmc_stats import *
 
-   def MonteCarlo(a, nmax):
-       E = 0.
-       N = 0.
-       for istep in range(nmax):
-           r = np.random.uniform(-5., 5., (3))
-           w = psi(a,r)
-           w = w*w
-           N += w
-           E += w * e_loc(a,r)
-       return E/N
+def MonteCarlo(a, nmax):
+     E = 0.
+     N = 0.
+     for istep in range(nmax):
+          r = np.random.uniform(-5., 5., (3))
+          w = psi(a,r)
+          w = w*w
+          N += w
+          E += w * e_loc(a,r)
+   return E/N
 
-   a = 0.9
-   nmax = 100000
-   X = [MonteCarlo(a,nmax) for i in range(30)]
-   E, deltaE = ave_error(X)
-   print(f"E = {E} +/- {deltaE}")
+a = 0.9
+nmax = 100000
+X = [MonteCarlo(a,nmax) for i in range(30)]
+E, deltaE = ave_error(X)
+print(f"E = {E} +/- {deltaE}")
diff --git a/variance_hydrogen.f90 b/variance_hydrogen.f90
index 63a4fb7..95eb20a 100644
--- a/variance_hydrogen.f90
+++ b/variance_hydrogen.f90
@@ -26,7 +26,6 @@ program variance_hydrogen
               r(3) = x(l)
               w = psi(a(j),r)
               w = w * w * delta
-
               energy = energy + w * e_loc(a(j), r)
               norm   = norm   + w 
            end do
@@ -44,7 +43,6 @@ program variance_hydrogen
               r(3) = x(l)
               w = psi(a(j),r)
               w = w * w * delta
-
               s2 = s2 + w * ( e_loc(a(j), r) - energy )**2
               norm   = norm   + w 
            end do
diff --git a/variance_hydrogen.py b/variance_hydrogen.py
index 071d6aa..4de7300 100644
--- a/variance_hydrogen.py
+++ b/variance_hydrogen.py
@@ -20,8 +20,8 @@ for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]:
                 El = e_loc(a, r)
                 E  += w * El
                 norm += w 
-                E = E / norm
-                s2 = 0.
+    E = E / norm
+    s2 = 0.
     for x in interval:
         r[0] = x
         for y in interval:
@@ -32,5 +32,5 @@ for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]:
                 w = w * w * delta
                 El = e_loc(a, r)
                 s2 += w * (El - E)**2
-                s2 = s2 / norm
-                print(f"a = {a} \t E = {E:10.8f}  \t  \sigma^2 = {s2:10.8f}")
+    s2 = s2 / norm
+    print(f"a = {a} \t E = {E:10.8f}  \t  \sigma^2 = {s2:10.8f}")
diff --git a/worg.css b/worg.css
index 50e4846..9f23339 100644
--- a/worg.css
+++ b/worg.css
@@ -306,6 +306,10 @@
 	/* font-lock-warning-face */
 	background-color: #e3e3f7;
     }
+    .exercise {
+	/* font-lock-warning-face */
+	background-color: #e3f7e3;
+    }
     .note {
 	/* font-lock-warning-face */
 	background-color: #f7f7d9;