From 09179024ea19b9207d566ecf1c3017cf1d2b7f85 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Fri, 29 Jan 2021 13:23:00 +0100
Subject: [PATCH] PDMC codes OK

---
 QMC.org               | 990 +++++++++++++++++++++++++++++-------------
 energy_hydrogen.f90   |  14 +-
 energy_hydrogen.py    |   6 +-
 hydrogen.f90          |  27 +-
 hydrogen.py           |   3 +-
 qmc_metropolis.f90    |  37 +-
 qmc_metropolis.py     |  35 +-
 qmc_stats.f90         |  26 +-
 qmc_stats.py          |   8 +-
 qmc_uniform.f90       |  24 +-
 qmc_uniform.py        |  14 +-
 variance_hydrogen.f90 |  31 +-
 variance_hydrogen.py  |  26 +-
 vmc.f90               |  14 +-
 vmc.py                |  10 +-
 vmc_metropolis.f90    |  72 ++-
 vmc_metropolis.py     |  55 ++-
 17 files changed, 951 insertions(+), 441 deletions(-)

diff --git a/QMC.org b/QMC.org
index 1fc2883..8264023 100644
--- a/QMC.org
+++ b/QMC.org
@@ -9,8 +9,7 @@
 
 #+OPTIONS: H:4 num:t toc:t \n:nil @:t ::t |:t ^:t -:t f:t *:t <:t
 #+OPTIONS: TeX:t LaTeX:t skip:nil d:nil todo:t pri:nil tags:not-in-toc
-# EXCLUDE_TAGS: Python solution
-# EXCLUDE_TAGS: Fortran solution
+# EXCLUDE_TAGS: solution
 
   #+BEGIN_SRC elisp :output none :exports none
 (setq org-latex-listings 'minted
@@ -152,7 +151,9 @@ def potential(r):
 double precision function potential(r)
   implicit none
   double precision, intent(in) :: r(3)
+
   ! TODO
+
 end function potential
      #+END_SRC
 
@@ -172,13 +173,17 @@ def potential(r):
 double precision function potential(r)
   implicit none
   double precision, intent(in) :: r(3)
-  double precision :: distance
+
+  double precision             :: distance
+
   distance = dsqrt( r(1)*r(1) + r(2)*r(2) + r(3)*r(3) )
+
   if (distance > 0.d0) then
-     potential = -1.d0 / dsqrt( r(1)*r(1) + r(2)*r(2) + r(3)*r(3) )
+     potential = -1.d0 / distance
   else
      stop 'potential at r=0.d0 diverges'
   end if
+
 end function potential
      #+END_SRC
 
@@ -200,7 +205,9 @@ def psi(a, r):
 double precision function psi(a, r)
   implicit none
   double precision, intent(in) :: a, r(3)
+
   ! TODO
+
 end function psi
      #+END_SRC
      
@@ -216,6 +223,7 @@ def psi(a, r):
 double precision function psi(a, r)
   implicit none
   double precision, intent(in) :: a, r(3)
+
   psi = dexp(-a * dsqrt( r(1)*r(1) + r(2)*r(2) + r(3)*r(3) ))
 end function psi
      #+END_SRC
@@ -268,7 +276,9 @@ def kinetic(a,r):
 double precision function kinetic(a,r)
   implicit none
   double precision, intent(in) :: a, r(3)
+
   ! TODO
+
 end function kinetic
      #+END_SRC
 
@@ -278,7 +288,8 @@ end function kinetic
 def kinetic(a,r):
     distance = np.sqrt(np.dot(r,r))
     assert (distance > 0.)
-    return -0.5 * (a**2 - (2.*a)/distance)
+
+    return a * (1./distance - 0.5 * a)
      #+END_SRC
 
     *Fortran*
@@ -286,13 +297,19 @@ def kinetic(a,r):
 double precision function kinetic(a,r)
   implicit none
   double precision, intent(in) :: a, r(3)
-  double precision :: distance
+
+  double precision             :: distance
+
   distance = dsqrt( r(1)*r(1) + r(2)*r(2) + r(3)*r(3) ) 
+
   if (distance > 0.d0) then
-     kinetic = -0.5d0 * (a*a - (2.d0*a) / distance)
+
+     kinetic =  a * (1.d0 / distance - 0.5d0 * a)
+
   else
      stop 'kinetic energy diverges at r=0'
   end if
+
 end function kinetic
      #+END_SRC
 
@@ -320,7 +337,9 @@ def e_loc(a,r):
 double precision function e_loc(a,r)
   implicit none
   double precision, intent(in) :: a, r(3)
+
   ! TODO
+
 end function e_loc
      #+END_SRC
    
@@ -336,8 +355,11 @@ def e_loc(a,r):
 double precision function e_loc(a,r)
   implicit none
   double precision, intent(in) :: a, r(3)
+
   double precision, external   :: kinetic, potential
+
   e_loc = kinetic(a,r) + potential(r)
+
 end function e_loc
      #+END_SRC
    
@@ -576,7 +598,9 @@ program energy_hydrogen
   end do
 
   do j=1,size(a)
+
      ! TODO
+
      print *, 'a = ', a(j), '    E = ', energy
   end do
 
@@ -586,7 +610,6 @@ end program energy_hydrogen
      To compile the Fortran and run it:
 
      #+begin_src sh :results output :exports code
-
 gfortran hydrogen.f90 energy_hydrogen.f90 -o energy_hydrogen
 ./energy_hydrogen 
      #+end_src
@@ -603,18 +626,22 @@ delta = (interval[1]-interval[0])**3
 r = np.array([0.,0.,0.])
 
 for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]:
-    E = 0.
+    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}")                
 
@@ -635,7 +662,7 @@ program energy_hydrogen
   implicit none
   double precision, external :: e_loc, psi
   double precision :: x(50), w, delta, energy, dx, r(3), a(6), norm
-  integer :: i, k, l, j
+  integer          :: i, k, l, j
 
   a = (/ 0.1d0, 0.2d0, 0.5d0, 1.d0, 1.5d0, 2.d0 /)
 
@@ -650,20 +677,28 @@ program energy_hydrogen
 
   do j=1,size(a)
      energy = 0.d0
-     norm = 0.d0
+     norm   = 0.d0
+     
      do i=1,size(x)
         r(1) = x(i)
+
         do k=1,size(x)
            r(2) = x(k)
+
            do l=1,size(x)
               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
+
         end do
+
      end do
+
      energy = energy / norm
      print *, 'a = ', a(j), '    E = ', energy
   end do
@@ -738,102 +773,179 @@ gfortran hydrogen.f90 energy_hydrogen.f90 -o energy_hydrogen
      #+begin_src python :results none :tangle none
 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.])
 
 for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]:
+
     # TODO
+
     print(f"a = {a} \t E = {E:10.8f} \t \sigma^2 = {s2:10.8f}")
     #+end_src
 
     *Fortran*
      #+begin_src f90 :tangle none
-program variance_hydrogen implicit none double precision, external ::
-  e_loc, psi double precision :: x(50), w, delta, energy, dx, r(3),
-  a(6), norm, s2 double precision :: e, energy2 integer :: i, k, l, j
+program variance_hydrogen
+  implicit none
+
+  double precision :: x(50), w, delta, energy, energy2
+  double precision :: dx, r(3), a(6), norm, e_tmp, s2
+  integer          :: i, k, l, j
+
+  double precision, external :: e_loc, psi
 
   a = (/ 0.1d0, 0.2d0, 0.5d0, 1.d0, 1.5d0, 2.d0 /)
-  
-  dx = 10.d0/(size(x)-1) do i=1,size(x) x(i) = -5.d0 + (i-1)*dx end do
-  
-  delta = dx**3
-  
-  r(:) = 0.d0
-  
-  do j=1,size(a) ! TODO print *, 'a = ', a(j), ' E = ', energy, ' s2 =
-  ', s2 end do
 
-end program variance_hydrogen #+end_src
+  dx = 10.d0/(size(x)-1)
+  do i=1,size(x)
+     x(i) = -5.d0 + (i-1)*dx
+  end do
+
+  do j=1,size(a)
+
+     ! TODO
+
+     print *, 'a = ', a(j), '    E = ', energy
+  end do
+
+end program variance_hydrogen
+     #+end_src
 
      To compile and run:
 
      #+begin_src sh :results output :exports both
 gfortran hydrogen.f90 variance_hydrogen.f90 -o variance_hydrogen
-./variance_hydrogen #+end_src
-
+./variance_hydrogen
+     #+end_src
+    
 **** Solution                                                      :solution:
     *Python*
-     #+begin_src python :results none :exports both
-import numpy as np from hydrogen import e_loc, psi
+     #+BEGIN_SRC python :results none :exports both
+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.])
 
-for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]: E = 0.  E2 = 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 El = e_loc(a,
-    r) E += w * El E2 += w * El*El norm += w E = E / norm E2 = E2 /
-    norm s2 = E2 - E*E print(f"a = {a} \t E = {E:10.8f} \t \sigma^2 =
-    {s2:10.8f}") #+end_src
+for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]:
+    E    = 0.
+    E2   = 0.
+    norm = 0.
 
-     #+RESULTS:
-     : a = 0.1 	 E = -0.24518439  	  \sigma^2 = 0.02696522
-     : a = 0.2 	 E = -0.26966058  	  \sigma^2 = 0.03719707
-     : a = 0.5 	 E = -0.38563576  	  \sigma^2 = 0.05318597
-     : a = 0.9 	 E = -0.49435710  	  \sigma^2 = 0.00577812
-     : a = 1.0 	 E = -0.50000000  	  \sigma^2 = 0.00000000
-     : a = 1.5 	 E = -0.39242967  	  \sigma^2 = 0.31449671
-     : a = 2.0 	 E = -0.08086981  	  \sigma^2 = 1.80688143
+    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_tmp = e_loc(a,r)
+                E    += w * e_tmp
+                E2   += w * e_tmp * e_tmp
+                norm += w 
+
+    E  = E  / norm
+    E2 = E2 / norm
+
+    s2 = E2 - E**2
+    print(f"a = {a} \t E = {E:10.8f} \t \sigma^2 = {s2:10.8f}")
+
+     #+end_src
      
+     #+RESULTS:
+     : a = 0.1 	 E = -0.24518439 	 \sigma^2 = 0.02696522
+     : a = 0.2 	 E = -0.26966058 	 \sigma^2 = 0.03719707
+     : a = 0.5 	 E = -0.38563576 	 \sigma^2 = 0.05318597
+     : a = 0.9 	 E = -0.49435710 	 \sigma^2 = 0.00577812
+     : a = 1.0 	 E = -0.50000000 	 \sigma^2 = 0.00000000
+     : a = 1.5 	 E = -0.39242967 	 \sigma^2 = 0.31449671
+     : a = 2.0 	 E = -0.08086981 	 \sigma^2 = 1.80688143
+
     *Fortran*
-     #+begin_src f90 
-program variance_hydrogen implicit none double precision, external ::
-  e_loc, psi double precision :: x(50), w, delta, energy, dx, r(3),
-  a(6), norm, s2 double precision :: e, energy2 integer :: i, k, l, j
+
+    #+begin_src f90 
+program variance_hydrogen
+  implicit none
+
+  double precision :: x(50), w, delta, energy, energy2
+  double precision :: dx, r(3), a(6), norm, e_tmp, s2
+  integer          :: i, k, l, j
+
+  double precision, external :: e_loc, psi
 
   a = (/ 0.1d0, 0.2d0, 0.5d0, 1.d0, 1.5d0, 2.d0 /)
 
-  dx = 10.d0/(size(x)-1) do i=1,size(x) x(i) = -5.d0 + (i-1)*dx end do
+  dx = 10.d0/(size(x)-1)
+  do i=1,size(x)
+     x(i) = -5.d0 + (i-1)*dx
+  end do
 
   delta = dx**3
 
   r(:) = 0.d0
 
-  do j=1,size(a) energy = 0.d0 energy2 = 0.d0 norm = 0.d0 do
-     i=1,size(x) r(1) = x(i) do k=1,size(x) r(2) = x(k) do l=1,size(x)
-     r(3) = x(l) w = psi(a(j),r) w = w * w * delta e = e_loc(a(j), r)
-     energy = energy + w * e energy2 = energy2 + w * e * e norm =
-     norm + w end do end do end do energy = energy / norm energy2 =
-     energy2 / norm s2 = energy2 - energy*energy print *, 'a = ',
-     a(j), ' E = ', energy, ' s2 = ', s2 end do
+  do j=1,size(a)
+     energy  = 0.d0
+     energy2 = 0.d0
+     norm    = 0.d0
 
-end program variance_hydrogen #+end_src
+     do i=1,size(x)
+        r(1) = x(i)
 
-     #+begin_src sh :results output :exports results
+        do k=1,size(x)
+           r(2) = x(k)
+
+           do l=1,size(x)
+              r(3) = x(l)
+
+              w = psi(a(j),r)
+              w = w * w * delta
+
+              e_tmp = e_loc(a(j), r)
+
+              energy  = energy  + w * e_tmp
+              energy2 = energy2 + w * e_tmp * e_tmp
+              norm   = norm     + w 
+           end do
+
+        end do
+
+     end do
+
+     energy  = energy  / norm
+     energy2 = energy2 / norm
+
+     s2 = energy2 - energy*energy
+
+     print *, 'a = ', a(j), ' E = ', energy, ' s2 = ', s2
+  end do
+
+end program variance_hydrogen
+    #+end_src
+
+    #+begin_src sh :results output :exports results
 gfortran hydrogen.f90 variance_hydrogen.f90 -o variance_hydrogen
-./variance_hydrogen #+end_src
+./variance_hydrogen
+    #+end_src
 
      #+RESULTS:
-     :  a =   0.10000000000000001       E =  -0.24518438948809140       s2 =    2.6965218719722767E-002
-     :  a =   0.20000000000000001       E =  -0.26966057967803236       s2 =    3.7197072370201284E-002
-     :  a =   0.50000000000000000       E =  -0.38563576125173815       s2 =    5.3185967578480653E-002
-     :  a =    1.0000000000000000       E =  -0.50000000000000000       s2 =    0.0000000000000000     
-     :  a =    1.5000000000000000       E =  -0.39242967082602065       s2 =   0.31449670909172917     
-     :  a =    2.0000000000000000       E =   -8.0869806678448772E-002  s2 =    1.8068814270846534     
+     :  a =   0.10000000000000001      E =  -0.24518438948809140       s2 =    2.6965218719722767E-002
+     :  a =   0.20000000000000001      E =  -0.26966057967803236       s2 =    3.7197072370201284E-002
+     :  a =   0.50000000000000000      E =  -0.38563576125173815       s2 =    5.3185967578480653E-002
+     :  a =    1.0000000000000000      E =  -0.50000000000000000       s2 =    0.0000000000000000     
+     :  a =    1.5000000000000000      E =  -0.39242967082602065       s2 =   0.31449670909172917     
+     :  a =    2.0000000000000000      E =   -8.0869806678448772E-002  s2 =    1.8068814270846534     
 
 * Variational Monte Carlo
 
@@ -890,15 +1002,17 @@ def ave_error(arr):
      #+END_SRC
 
     *Fortran*
-        #+BEGIN_SRC f90
+    #+BEGIN_SRC f90 :tangle none
 subroutine ave_error(x,n,ave,err)
   implicit none
   integer, intent(in)           :: n 
   double precision, intent(in)  :: x(n) 
   double precision, intent(out) :: ave, err
+
   ! TODO
+
 end subroutine ave_error
-        #+END_SRC
+    #+END_SRC
    
 **** Solution                                                      :solution:
     *Python*
@@ -907,32 +1021,43 @@ from math import sqrt
 def ave_error(arr):
     M = len(arr)
     assert(M>0)
+
     if M == 1:
-        return (arr[0], 0.)
+        average = arr[0]
+        error   = 0.
+
     else:
         average = sum(arr)/M
         variance = 1./(M-1) * sum( [ (x - average)**2 for x in arr ] )
         error = sqrt(variance/M)
-        return (average, error)
+
+    return (average, error)
      #+END_SRC
 
     *Fortran*
         #+BEGIN_SRC f90 :exports both
 subroutine ave_error(x,n,ave,err)
   implicit none
+
   integer, intent(in)           :: n 
   double precision, intent(in)  :: x(n) 
   double precision, intent(out) :: ave, err
-  double precision :: variance
+
+  double precision              :: variance
+
   if (n < 1) then
      stop 'n<1 in ave_error'
+
   else if (n == 1) then
      ave = x(1)
      err = 0.d0
+
   else
-     ave = sum(x(:)) / dble(n)
-     variance = sum( (x(:) - ave)**2 ) / dble(n-1)
-     err = dsqrt(variance/dble(n))
+     ave      = sum(x(:)) / dble(n)
+
+     variance = sum((x(:) - ave)**2) / dble(n-1)
+     err      = dsqrt(variance/dble(n))
+
   endif
 end subroutine ave_error
         #+END_SRC
@@ -986,9 +1111,11 @@ from qmc_stats import *
 def MonteCarlo(a, nmax):
      # TODO
 
-a = 0.9
+a    = 0.9
 nmax = 100000
+
 #TODO
+
 print(f"E = {E} +/- {deltaE}")
      #+END_SRC
 
@@ -1012,8 +1139,7 @@ subroutine uniform_montecarlo(a,nmax,energy)
   integer*8       , intent(in)  :: nmax 
   double precision, intent(out) :: energy
 
-  integer*8 :: istep
-
+  integer*8        :: istep
   double precision :: norm, r(3), w
 
   double precision, external :: e_loc, psi
@@ -1027,12 +1153,14 @@ program qmc
   integer*8       , parameter :: nmax = 100000
   integer         , parameter :: nruns = 30
 
-  integer :: irun
+  integer          :: irun
   double precision :: X(nruns)
   double precision :: ave, err
 
   !TODO
+
   print *, 'E = ', ave, '+/-', err
+
 end program qmc
      #+END_SRC
 
@@ -1050,18 +1178,24 @@ from qmc_stats import *
 def MonteCarlo(a, nmax):
      energy = 0.
      normalization = 0.
+
      for istep in range(nmax):
           r = np.random.uniform(-5., 5., (3))
+
           w = psi(a,r)
           w = w*w
-          normalization += w
-          energy += w * e_loc(a,r)
-     return energy/normalization
 
-a = 0.9
+          energy        += w * e_loc(a,r)
+          normalization += w
+
+     return energy / normalization
+
+a    = 0.9
 nmax = 100000
+
 X = [MonteCarlo(a,nmax) for i in range(30)]
 E, deltaE = ave_error(X)
+
 print(f"E = {E} +/- {deltaE}")
      #+END_SRC
 
@@ -1083,39 +1217,47 @@ subroutine uniform_montecarlo(a,nmax,energy)
   integer*8       , intent(in)  :: nmax 
   double precision, intent(out) :: energy
 
-  integer*8 :: istep
-
+  integer*8        :: istep
   double precision :: norm, r(3), w
 
   double precision, external :: e_loc, psi
 
   energy = 0.d0
   norm   = 0.d0
+
   do istep = 1,nmax
+
      call random_number(r)
      r(:) = -5.d0 + 10.d0*r(:)
+
      w = psi(a,r)
      w = w*w
-     norm = norm + w
+
      energy = energy + w * e_loc(a,r)
+     norm   = norm   + w
+
   end do
+
   energy = energy / norm
+
 end subroutine uniform_montecarlo
 
 program qmc
   implicit none
-  double precision, parameter :: a = 0.9
-  integer*8       , parameter :: nmax = 100000
+  double precision, parameter :: a     = 0.9
+  integer*8       , parameter :: nmax  = 100000
   integer         , parameter :: nruns = 30
 
-  integer :: irun
+  integer          :: irun
   double precision :: X(nruns)
   double precision :: ave, err
 
   do irun=1,nruns
-     call uniform_montecarlo(a,nmax,X(irun))
+     call uniform_montecarlo(a, nmax, X(irun))
   enddo
-  call ave_error(X,nruns,ave,err)
+
+  call ave_error(X, nruns, ave, err)
+
   print *, 'E = ', ave, '+/-', err
 end program qmc
      #+END_SRC
@@ -1126,7 +1268,7 @@ gfortran hydrogen.f90 qmc_stats.f90 qmc_uniform.f90 -o qmc_uniform
      #+end_src
 
      #+RESULTS:
-     :  E =  -0.49588321986667677      +/-   7.1758863546737969E-004
+     :  E =  -0.49518773675598715      +/-   5.2391494923686175E-004
 
 ** Metropolis sampling with $\Psi^2$
    :PROPERTIES:
@@ -1212,13 +1354,17 @@ from hydrogen  import *
 from qmc_stats import *
 
 def MonteCarlo(a,nmax,dt):
+
     # TODO
+
     return energy/nmax, N_accep/nmax
 
+
 # Run simulation
-a = 0.9
+a    = 0.9
 nmax = 100000
-dt = 1.3
+dt   = #TODO
+
 X0 = [ MonteCarlo(a,nmax,dt) for i in range(30)]
 
 # Energy
@@ -1242,24 +1388,25 @@ subroutine metropolis_montecarlo(a,nmax,dt,energy,accep)
   double precision, intent(out) :: energy
   double precision, intent(out) :: accep
 
-  integer*8 :: istep
-
+  integer*8        :: istep
+  integer*8        :: n_accep
   double precision :: r_old(3), r_new(3), psi_old, psi_new
   double precision :: v, ratio
-  integer*8        :: n_accep
+
   double precision, external :: e_loc, psi, gaussian
 
   ! TODO
+
 end subroutine metropolis_montecarlo
 
 program qmc
   implicit none
-  double precision, parameter :: a = 0.9d0
-  double precision, parameter :: dt = 1.3d0
-  integer*8       , parameter :: nmax = 100000
+  double precision, parameter :: a     = 0.9d0
+  double precision, parameter :: dt    = ! TODO
+  integer*8       , parameter :: nmax  = 100000
   integer         , parameter :: nruns = 30
 
-  integer :: irun
+  integer          :: irun
   double precision :: X(nruns), Y(nruns)
   double precision :: ave, err
 
@@ -1272,6 +1419,7 @@ program qmc
 
   call ave_error(Y,nruns,ave,err)
   print *, 'A = ', ave, '+/-', err
+
 end program qmc
      #+END_SRC
 
@@ -1287,26 +1435,33 @@ from hydrogen  import *
 from qmc_stats import *
 
 def MonteCarlo(a,nmax,dt):
-    energy = 0.
+    energy  = 0.
     N_accep = 0
+
     r_old = np.random.uniform(-dt, dt, (3))
     psi_old = psi(a,r_old)
+
     for istep in range(nmax):
+        energy += e_loc(a,r_old)
+
         r_new = r_old + np.random.uniform(-dt,dt,(3))
         psi_new = psi(a,r_new)
+
         ratio = (psi_new / psi_old)**2
-        v = np.random.uniform()
-        if v <= ratio:
+
+        if np.random.uniform() <= ratio:
             N_accep += 1
-            r_old = r_new
+
+            r_old   = r_new
             psi_old = psi_new
-        energy += e_loc(a,r_old)
+
     return energy/nmax, N_accep/nmax
 
 # Run simulation
-a = 0.9
+a    = 0.9
 nmax = 100000
-dt = 1.3
+dt   = 1.3
+
 X0 = [ MonteCarlo(a,nmax,dt) for i in range(30)]
 
 # Energy
@@ -1334,33 +1489,45 @@ subroutine metropolis_montecarlo(a,nmax,dt,energy,accep)
   double precision, intent(out) :: energy
   double precision, intent(out) :: accep
 
-  integer*8 :: istep
-
   double precision :: r_old(3), r_new(3), psi_old, psi_new
   double precision :: v, ratio
   integer*8        :: n_accep
+  integer*8        :: istep
+
   double precision, external :: e_loc, psi, gaussian
 
-  energy = 0.d0
+  energy  = 0.d0
   n_accep = 0_8
+
   call random_number(r_old)
   r_old(:) = dt * (2.d0*r_old(:) - 1.d0)
   psi_old = psi(a,r_old)
+
   do istep = 1,nmax
+     energy = energy + e_loc(a,r_old)
+
      call random_number(r_new)
-     r_new(:) = r_old(:) + dt * (2.d0*r_new(:) - 1.d0)
+     r_new(:) = r_old(:) + dt*(2.d0*r_new(:) - 1.d0)
+
      psi_new = psi(a,r_new)
+
      ratio = (psi_new / psi_old)**2
      call random_number(v)
+
      if (v <= ratio) then
+
+        n_accep = n_accep + 1_8
+
         r_old(:) = r_new(:)
         psi_old = psi_new
-        n_accep = n_accep + 1_8
+
      endif
-     energy = energy + e_loc(a,r_old)
+
   end do
+
   energy = energy / dble(nmax)
   accep  = dble(n_accep) / dble(nmax)
+
 end subroutine metropolis_montecarlo
 
 program qmc
@@ -1370,7 +1537,7 @@ program qmc
   integer*8       , parameter :: nmax = 100000
   integer         , parameter :: nruns = 30
 
-  integer :: irun
+  integer          :: irun
   double precision :: X(nruns), Y(nruns)
   double precision :: ave, err
 
@@ -1383,6 +1550,7 @@ program qmc
 
   call ave_error(Y,nruns,ave,err)
   print *, 'A = ', ave, '+/-', err
+
 end program qmc
      #+END_SRC
 
@@ -1391,8 +1559,8 @@ gfortran hydrogen.f90 qmc_stats.f90 qmc_metropolis.f90 -o qmc_metropolis
 ./qmc_metropolis
      #+end_src
      #+RESULTS:
-     :  E =  -0.49515370205041676      +/-   1.7660819245720729E-004
-     :  A =   0.51713866666666664      +/-   3.7072551835783688E-004
+     :  E =  -0.49503130891988767      +/-   1.7160104275040037E-004
+     :  A =   0.51695266666666673      +/-   4.0445505648997396E-004
 
 ** Gaussian random number generator
    
@@ -1419,6 +1587,7 @@ subroutine random_gauss(z,n)
   integer :: i
 
   call random_number(u)
+
   if (iand(n,1) == 0) then
      ! n is even
      do i=1,n,2
@@ -1426,6 +1595,7 @@ subroutine random_gauss(z,n)
         z(i+1) = z(i) * dsin( two_pi*u(i+1) )
         z(i)   = z(i) * dcos( two_pi*u(i+1) )
      end do
+
   else
      ! n is odd
      do i=1,n-1,2
@@ -1433,9 +1603,12 @@ subroutine random_gauss(z,n)
         z(i+1) = z(i) * dsin( two_pi*u(i+1) )
         z(i)   = z(i) * dcos( two_pi*u(i+1) )
      end do
+
      z(n)   = dsqrt(-2.d0*dlog(u(n))) 
      z(n)   = z(n) * dcos( two_pi*u(n+1) )
+
   end if
+
 end subroutine random_gauss
    #+END_SRC
 
@@ -1532,7 +1705,9 @@ subroutine drift(a,r,b)
   implicit none
   double precision, intent(in)  :: a, r(3)
   double precision, intent(out) :: b(3)
+
   ! TODO
+
 end subroutine drift
      #+END_SRC
 
@@ -1550,9 +1725,12 @@ 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
+  b(:)   = r(:) * ar_inv
+
 end subroutine drift
      #+END_SRC
 
@@ -1572,9 +1750,10 @@ def MonteCarlo(a,nmax,dt):
    # TODO
 
 # Run simulation
-a = 0.9
+a    = 0.9
 nmax = 100000
-dt = 1.3
+dt   = # TODO
+
 X0 = [ MonteCarlo(a,nmax,dt) for i in range(30)]
 
 # Energy
@@ -1590,17 +1769,19 @@ print(f"A = {A} +/- {deltaA}")
 
     *Fortran*
      #+BEGIN_SRC f90 :tangle none
-subroutine variational_montecarlo(a,dt,nmax,energy,accep_rate)
+subroutine variational_montecarlo(a,dt,nmax,energy,accep)
   implicit none
   double precision, intent(in)  :: a, dt
   integer*8       , intent(in)  :: nmax 
-  double precision, intent(out) :: energy, accep_rate
+  double precision, intent(out) :: energy, accep
 
-  integer*8 :: istep
+  integer*8        :: istep
+  integer*8        :: n_accep
   double precision :: sq_dt, chi(3)
   double precision :: psi_old, psi_new
   double precision :: r_old(3), r_new(3)
   double precision :: d_old(3), d_new(3)
+
   double precision, external :: e_loc, psi
 
   ! TODO
@@ -1609,22 +1790,25 @@ end subroutine variational_montecarlo
 
 program qmc
   implicit none
-  double precision, parameter :: a = 0.9
-  double precision, parameter :: dt = 1.0
-  integer*8       , parameter :: nmax = 100000
+  double precision, parameter :: a     = 0.9
+  double precision, parameter :: dt    = ! TODO
+  integer*8       , parameter :: nmax  = 100000
   integer         , parameter :: nruns = 30
 
-  integer :: irun
+  integer          :: irun
   double precision :: X(nruns), accep(nruns)
   double precision :: ave, err
 
   do irun=1,nruns
      call variational_montecarlo(a,dt,nmax,X(irun),accep(irun))
   enddo
+
   call ave_error(X,nruns,ave,err)
   print *, 'E = ', ave, '+/-', err
+
   call ave_error(accep,nruns,ave,err)
   print *, 'A = ', ave, '+/-', err
+
 end program qmc
      #+END_SRC
 
@@ -1640,38 +1824,49 @@ from hydrogen  import *
 from qmc_stats import *
 
 def MonteCarlo(a,nmax,dt):
-    energy = 0.
-    accep_rate = 0.
     sq_dt = np.sqrt(dt)
-    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)
+
+    energy  = 0.
+    N_accep = 0
+
+    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):
         chi = np.random.normal(loc=0., scale=1.0, size=(3))
-        r_new = r_old + dt * d_old + sq_dt * chi
-        d_new = drift(a,r_new)
-        d2_new = np.dot(d_new,d_new)
+
+        energy += e_loc(a,r_old)
+
+        r_new   = r_old + dt * d_old + sq_dt * chi
+        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))
+        prod    = np.dot((d_new + d_old), (r_new - r_old))
         argexpo = 0.5 * (d2_new - d2_old)*dt + prod
+
         q = psi_new / psi_old
         q = np.exp(-argexpo) * q*q
-        if np.random.uniform() < q:
-            accep_rate += 1.
-            r_old = r_new
-            d_old = d_new
-            d2_old = d2_new
+
+        if np.random.uniform() <= q:
+            N_accep += 1
+
+            r_old   = r_new
+            d_old   = d_new
+            d2_old  = d2_new
             psi_old = psi_new
-        energy += e_loc(a,r_old)
+
     return energy/nmax, accep_rate/nmax
 
 
 # Run simulation
-a = 0.9
+a    = 0.9
 nmax = 100000
-dt = 1.3
+dt   = 1.3
+
 X0 = [ MonteCarlo(a,nmax,dt) for i in range(30)]
 
 # Energy
@@ -1691,74 +1886,100 @@ print(f"A = {A} +/- {deltaA}")
    
     *Fortran*
      #+BEGIN_SRC f90
-subroutine variational_montecarlo(a,dt,nmax,energy,accep_rate)
+subroutine variational_montecarlo(a,dt,nmax,energy,accep)
   implicit none
   double precision, intent(in)  :: a, dt
   integer*8       , intent(in)  :: nmax 
-  double precision, intent(out) :: energy, accep_rate
+  double precision, intent(out) :: energy, accep
 
-  integer*8 :: istep
+  integer*8        :: istep
+  integer*8        :: n_accep
   double precision :: sq_dt, chi(3), d2_old, prod, u
   double precision :: psi_old, psi_new, d2_new, argexpo, q
   double precision :: r_old(3), r_new(3)
   double precision :: d_old(3), d_new(3)
+
   double precision, external :: e_loc, psi
 
   sq_dt = dsqrt(dt)
 
   ! Initialization
-  energy = 0.d0
-  accep_rate = 0.d0
+  energy  = 0.d0
+  n_accep = 0_8
+
   call random_gauss(r_old,3)
+
   call drift(a,r_old,d_old)
-  d2_old = d_old(1)*d_old(1) + d_old(2)*d_old(2) + d_old(3)*d_old(3)
+  d2_old  = d_old(1)*d_old(1) + &
+            d_old(2)*d_old(2) + &
+            d_old(3)*d_old(3)
+
   psi_old = psi(a,r_old)
 
   do istep = 1,nmax
+     energy = energy + e_loc(a,r_old)
+
      call random_gauss(chi,3)
-     r_new(:) = r_old(:) + dt * d_old(:) + chi(:)*sq_dt
+     r_new(:) = r_old(:) + dt*d_old(:) + chi(:)*sq_dt
+
      call drift(a,r_new,d_new)
-     d2_new = d_new(1)*d_new(1) + d_new(2)*d_new(2) + d_new(3)*d_new(3)
+     d2_new = d_new(1)*d_new(1) + &
+              d_new(2)*d_new(2) + &
+              d_new(3)*d_new(3)
+
      psi_new = psi(a,r_new)
+
      ! Metropolis
      prod = (d_new(1) + d_old(1))*(r_new(1) - r_old(1)) + &
             (d_new(2) + d_old(2))*(r_new(2) - r_old(2)) + &
             (d_new(3) + d_old(3))*(r_new(3) - r_old(3))
+
      argexpo = 0.5d0 * (d2_new - d2_old)*dt + prod
+
      q = psi_new / psi_old
      q = dexp(-argexpo) * q*q
+
      call random_number(u)
-     if (u<q) then
-        accep_rate = accep_rate + 1.d0
+
+     if (u <= q) then
+
+        n_accep = n_accep + 1_8
+
         r_old(:) = r_new(:)
         d_old(:) = d_new(:)
-        d2_old = d2_new
-        psi_old = psi_new
+        d2_old   = d2_new
+        psi_old  = psi_new
+
      end if
-     energy = energy + e_loc(a,r_old)
+
   end do
+
   energy = energy / dble(nmax)
-  accep_rate = dble(accep_rate) / dble(nmax)
+  accep  = dble(n_accep) / dble(nmax)
+
 end subroutine variational_montecarlo
 
 program qmc
   implicit none
-  double precision, parameter :: a = 0.9
-  double precision, parameter :: dt = 1.0
-  integer*8       , parameter :: nmax = 100000
+  double precision, parameter :: a     = 0.9
+  double precision, parameter :: dt    = 1.0
+  integer*8       , parameter :: nmax  = 100000
   integer         , parameter :: nruns = 30
 
-  integer :: irun
+  integer          :: irun
   double precision :: X(nruns), accep(nruns)
   double precision :: ave, err
 
   do irun=1,nruns
      call variational_montecarlo(a,dt,nmax,X(irun),accep(irun))
   enddo
+
   call ave_error(X,nruns,ave,err)
   print *, 'E = ', ave, '+/-', err
+
   call ave_error(accep,nruns,ave,err)
   print *, 'A = ', ave, '+/-', err
+
 end program qmc
      #+END_SRC
 
@@ -1768,15 +1989,10 @@ gfortran hydrogen.f90 qmc_stats.f90 vmc_metropolis.f90 -o vmc_metropolis
      #+end_src
 
      #+RESULTS:
-     :  E =  -0.49495906384751226      +/-   1.5257646086088266E-004
-     :  A =   0.78861366666666666      +/-   3.7855335138754813E-004
+     :  E =  -0.49497258331144794      +/-   1.0973395750688713E-004
+     :  A =   0.78839866666666658      +/-   3.2503783452043152E-004
      
-* Diffusion Monte Carlo
-   :PROPERTIES:
-   :header-args:python: :tangle dmc.py
-   :header-args:f90: :tangle dmc.f90
-   :END:
-   
+* Diffusion Monte Carlo                                            :solution:
 
    
 ** Schrödinger equation in imaginary time
@@ -1861,16 +2077,14 @@ gfortran hydrogen.f90 qmc_stats.f90 vmc_metropolis.f90 -o vmc_metropolis
     practice.
     Fortunately, if we multiply the Schrödinger equation by a chosen
     /trial wave function/ $\Psi_T(\mathbf{r})$ (Hartree-Fock, Kohn-Sham
-    determinant, CI wave function, /etc/), one obtains (see appendix
-    for details)
+    determinant, CI wave function, /etc/), one obtains
 
   \[
     -\frac{\partial \psi(\mathbf{r},\tau)}{\partial \tau} \Psi_T(\mathbf{r}) =
     \left[ -\frac{1}{2} \Delta \psi(\mathbf{r},\tau) + V(\mathbf{r}) \psi(\mathbf{r},\tau) \right] \Psi_T(\mathbf{r}) 
   \]
 
-  Defining $\Pi(\mathbf{r},\tau) = \psi(\mathbf{r},\tau)
-  \Psi_T(\mathbf{r})$, 
+  Defining $\Pi(\mathbf{r},\tau) = \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r})$, (see appendix for details)
 
   \[
   -\frac{\partial \Pi(\mathbf{r},\tau)}{\partial \tau}
@@ -1892,12 +2106,11 @@ gfortran hydrogen.f90 qmc_stats.f90 vmc_metropolis.f90 -o vmc_metropolis
   This equation generates the /N/-electron density $\Pi$, which is the
   product of the ground state with the trial wave function. It
   introduces the constraint that $\Pi(\mathbf{r},\tau)=0$ where
-  $\Psi_T(\mathbf{r})=0$. If the wave function has the same sign
-  everywhere, as in the hydrogen atom or the H_2 molecule, this
+  $\Psi_T(\mathbf{r})=0$. In the few cases where the wave function has no nodes,
+  such as in the hydrogen atom or the H_2 molecule, this
   constraint is harmless and we can obtain the exact energy. But for
-  systems where the wave function has nodes (systems with 3 or more
-  electrons), this scheme introduces an error known as the /fixed node
-  error/.
+  systems where the wave function has nodes, this scheme introduces an
+  error known as the /fixed node error/.
 
 *** Appendix : Details of the Derivation
     
@@ -1958,7 +2171,6 @@ gfortran hydrogen.f90 qmc_stats.f90 vmc_metropolis.f90 -o vmc_metropolis
    \infty} E(\tau).
    \]
    
-   Indeed, this is equivalent to
 
    \[
    E(\tau)   =   \frac{\int \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) E_L(\mathbf{r}) d\mathbf{r}}
@@ -1982,23 +2194,26 @@ gfortran hydrogen.f90 qmc_stats.f90 vmc_metropolis.f90 -o vmc_metropolis
    \]
 
    So computing the energy within DMC consists in generating the
-   density $\Pi$, and simply averaging the local energies computed
-   with the trial wave function.
+   density $\Pi$ with random walks, and simply averaging the local
+   energies computed with the trial wave function.
    
-** Pure Diffusion Monte Carlo
+** Pure Diffusion Monte Carlo (PDMC)
 
    Instead of having a variable number of particles to simulate the
    branching process, one can choose to sample $[\Psi_T(\mathbf{r})]^2$ instead of
    $\psi(\mathbf{r},\tau) \Psi_T(\mathbf{r})$, and consider the term
    $\exp \left( -\delta t\,(  E_L(\mathbf{r}) - E_{\text{ref}}  \right)$ as a
-   cumulative product of branching weights:
+   cumulative product of weights:
 
    \[
    W(\mathbf{r}_n, \tau)
    = \exp \left( \int_0^\tau - (E_L(\mathbf{r}_t) - E_{\text{ref}}) dt \right)
-   \approx \prod_{i=1}^{n} \exp \left( -\delta t\, (E_L(\mathbf{r}_i) - E_{\text{ref}}) \right).
+   \approx \prod_{i=1}^{n} \exp \left( -\delta t\,
+   (E_L(\mathbf{r}_i) - E_{\text{ref}}) \right) = 
+   \prod_{i=1}^{n} w(\mathbf{r}_i)
    \]
 
+   where $\mathbf{r}_i$ are the coordinates along the trajectory.
    The wave function becomes
 
    \[
@@ -2017,9 +2232,14 @@ gfortran hydrogen.f90 qmc_stats.f90 vmc_metropolis.f90 -o vmc_metropolis
    This algorithm is less stable than the branching algorithm: it 
    requires to have a value of $E_\text{ref}$ which is close to the
    fixed-node energy, and a good trial wave function. Its big
-   advantage is that it is very easy to program starting from a VMC code.
+   advantage is that it is very easy to program starting from a VMC
+   code, so this is what we will do in the next section.
    
-** TODO Hydrogen atom
+** Hydrogen atom
+  :PROPERTIES:
+  :header-args:python: :tangle pdmc.py
+  :header-args:f90: :tangle pdmc.f90
+  :END:
    
 *** Exercise 
 
@@ -2029,145 +2249,298 @@ gfortran hydrogen.f90 qmc_stats.f90 vmc_metropolis.f90 -o vmc_metropolis
      energy of H for any value of $a$.
      #+end_exercise
    
-**** Python
+    *Python*
          #+BEGIN_SRC python :results output
 from hydrogen  import *
 from qmc_stats import *
 
-def MonteCarlo(a,nmax,tau,Eref):
-    energy = 0.
-    normalization = 0.
-    accep_rate = 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)
-    w = 1.0
-    for istep in range(nmax):
-        chi = np.random.normal(loc=0., scale=1.0, size=(3))
-        el = e_loc(a,r_old)
-        w *= np.exp(-tau*(el - Eref))
-        normalization += w
-        energy += w * el
+def MonteCarlo(a, nmax, dt, Eref):
+    # TODO
 
-        r_new = r_old + tau * d_old + sq_tau * chi
-        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
-        # PDMC weight
-        if np.random.uniform() < q:
-            accep_rate += 1
-            r_old = r_new
-            d_old = d_new
-            d2_old = d2_new
-            psi_old = psi_new
-    return energy/normalization, accep_rate/nmax
-
-
-a = 0.9
-nmax = 10000
-tau = .1
+# Run simulation
+a     = 0.9
+nmax  = 100000
+dt    = 0.01
 E_ref = -0.5
-X = [MonteCarlo(a,nmax,tau,E_ref) for i in range(30)]
-E, deltaE = ave_error([x[0] for x in X])
-A, deltaA = ave_error([x[1] for x in X])
-print(f"E = {E} +/- {deltaE}\nA = {A} +/- {deltaA}")
+
+X0 = [ MonteCarlo(a, nmax, dt, E_ref) for i in range(30)]
+
+# Energy
+X = [ x for (x, _) in X0 ]
+E, deltaE = ave_error(X)
+print(f"E = {E} +/- {deltaE}")
+
+# Acceptance rate
+X = [ x for (_, x) in X0 ]
+A, deltaA = ave_error(X)
+print(f"A = {A} +/- {deltaA}")
          #+END_SRC
 
-         #+RESULTS:
-         : E = -0.49654807434947584 +/- 0.0006868522447409156
-         : A = 0.9876193891840709 +/- 0.00041857361650995804
-   
-**** Fortran
-     #+BEGIN_SRC f90
-subroutine variational_montecarlo(a,tau,nmax,energy,accep_rate)
+    *Fortran*
+     #+BEGIN_SRC f90 :tangle none
+subroutine pdmc(a, dt, nmax, energy, accep, tau, E_ref)
   implicit none
-  double precision, intent(in)  :: a, tau
+  double precision, intent(in)  :: a, dt, tau
   integer*8       , intent(in)  :: nmax 
-  double precision, intent(out) :: energy, accep_rate
+  double precision, intent(out) :: energy, accep
+  double precision, intent(in)  :: E_ref
 
-  integer*8 :: istep
-  double precision :: norm, sq_tau, chi(3), d2_old, prod, u
-  double precision :: psi_old, psi_new, d2_new, argexpo, q
+  integer*8        :: istep
+  integer*8        :: n_accep
+  double precision :: sq_dt, chi(3)
+  double precision :: psi_old, psi_new
   double precision :: r_old(3), r_new(3)
   double precision :: d_old(3), d_new(3)
+
   double precision, external :: e_loc, psi
 
-  sq_tau = dsqrt(tau)
-  
-  ! Initialization
-  energy = 0.d0
-  norm   = 0.d0
-  accep_rate = 0.d0
-  call random_gauss(r_old,3)
-  call drift(a,r_old,d_old)
-  d2_old = d_old(1)*d_old(1) + d_old(2)*d_old(2) + d_old(3)*d_old(3)
-  psi_old = psi(a,r_old)
+  ! TODO
 
-  do istep = 1,nmax
-     call random_gauss(chi,3)
-     r_new(:) = r_old(:) + tau * d_old(:) + chi(:)*sq_tau
-     call drift(a,r_new,d_new)
-     d2_new = d_new(1)*d_new(1) + d_new(2)*d_new(2) + d_new(3)*d_new(3)
-     psi_new = psi(a,r_new)
-     ! Metropolis
-     prod = (d_new(1) + d_old(1))*(r_new(1) - r_old(1)) + &
-            (d_new(2) + d_old(2))*(r_new(2) - r_old(2)) + &
-            (d_new(3) + d_old(3))*(r_new(3) - r_old(3))
-     argexpo = 0.5d0 * (d2_new - d2_old)*tau + prod
-     q = psi_new / psi_old
-     q = dexp(-argexpo) * q*q
-     call random_number(u)
-     if (u<q) then
-        accep_rate = accep_rate + 1.d0
-        r_old(:) = r_new(:)
-        d_old(:) = d_new(:)
-        d2_old = d2_new
-        psi_old = psi_new
-     end if
-     norm = norm + 1.d0
-     energy = energy + e_loc(a,r_old)
-  end do
-  energy = energy / norm
-  accep_rate = accep_rate / norm
-end subroutine variational_montecarlo
+end subroutine pdmc
 
 program qmc
   implicit none
-  double precision, parameter :: a = 0.9
-  double precision, parameter :: tau = 1.0
-  integer*8       , parameter :: nmax = 100000
+  double precision, parameter :: a     = 0.9
+  double precision, parameter :: dt    = 0.1d0
+  double precision, parameter :: E_ref = -0.5d0
+  double precision, parameter :: tau   = 10.d0
+  integer*8       , parameter :: nmax  = 100000
   integer         , parameter :: nruns = 30
 
-  integer :: irun
+  integer          :: irun
   double precision :: X(nruns), accep(nruns)
   double precision :: ave, err
 
   do irun=1,nruns
-     call variational_montecarlo(a,tau,nmax,X(irun),accep(irun))
+     call pdmc(a, dt, nmax, X(irun), accep(irun), tau, E_ref)
   enddo
+
   call ave_error(X,nruns,ave,err)
   print *, 'E = ', ave, '+/-', err
+
   call ave_error(accep,nruns,ave,err)
   print *, 'A = ', ave, '+/-', err
+
 end program qmc
-         #+END_SRC
+     #+END_SRC
 
-         #+begin_src sh :results output :exports both
-gfortran hydrogen.f90 qmc_stats.f90 vmc_metropolis.f90 -o vmc_metropolis
-./vmc_metropolis
-         #+end_src
+     #+begin_src sh :results output :exports code
+gfortran hydrogen.f90 qmc_stats.f90 pdmc.f90 -o pdmc
+./pdmc
+     #+end_src
 
-         #+RESULTS:
-         :  E =  -0.49499990423528023      +/-   1.5958250761863871E-004
-         :  A =   0.78861366666666655      +/-   3.5096729498002445E-004
+**** Solution                                                      :solution:
+
+     *Python*
+     #+BEGIN_SRC python :results output
+from hydrogen  import *
+from qmc_stats import *
+
+def MonteCarlo(a, nmax, dt, tau, Eref):
+    sq_dt = np.sqrt(dt)
+
+    energy        = 0.
+    N_accep       = 0
+    normalization = 0.
+
+    w = 1.
+    tau_current = 0.
+
+    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):
+        el = e_loc(a,r_old)
+        w *= np.exp(-dt*(el - Eref))
+
+        normalization += w
+        energy        += w * el
+
+        tau_current += dt
+
+        # Reset when tau is reached
+        if tau_current >= tau:
+            w = 1.
+            tau_current = 0.
+
+        chi = np.random.normal(loc=0., scale=1.0, size=(3))
+
+        r_new = r_old + dt * d_old + sq_dt * chi
+        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)*dt + prod
+
+        q = psi_new / psi_old
+        q = np.exp(-argexpo) * q*q
+
+        if np.random.uniform() <= q:
+            N_accep += 1
+            r_old   = r_new
+            d_old   = d_new
+            d2_old  = d2_new
+            psi_old = psi_new
+
+    return energy/normalization, N_accep/nmax
+
+
+# Run simulation
+a     = 0.9
+nmax  = 100000
+dt    = 0.1
+tau   = 10.
+E_ref = -0.5
+
+X0 = [ MonteCarlo(a, nmax, dt, tau, E_ref) for i in range(30)]
+
+# Energy
+X = [ x for (x, _) in X0 ]
+E, deltaE = ave_error(X)
+print(f"E = {E} +/- {deltaE}")
+
+# Acceptance rate
+X = [ x for (_, x) in X0 ]
+A, deltaA = ave_error(X)
+print(f"A = {A} +/- {deltaA}")
+     #+END_SRC
+
+     #+RESULTS:
+     : E = -0.5006126340155459 +/- 0.00042555955652480285
+     : A = 0.9878509999999998 +/- 6.920791596475006e-05
+
+     *Fortran* 
+     #+BEGIN_SRC f90
+subroutine pdmc(a, dt, nmax, energy, accep, tau, E_ref)
+  implicit none
+  double precision, intent(in)  :: a, dt, tau
+  integer*8       , intent(in)  :: nmax 
+  double precision, intent(out) :: energy, accep
+  double precision, intent(in)  :: E_ref
+
+  integer*8        :: istep
+  integer*8        :: n_accep
+  double precision :: sq_dt, chi(3), d2_old, prod, u
+  double precision :: psi_old, psi_new, d2_new, argexpo, q
+  double precision :: r_old(3), r_new(3)
+  double precision :: d_old(3), d_new(3)
+  double precision :: e, w, norm, tau_current
+
+  double precision, external :: e_loc, psi
+
+  sq_dt = dsqrt(dt)
+
+  ! Initialization
+  energy  = 0.d0
+  n_accep = 0_8
+  norm    = 0.d0
+
+  w           = 1.d0
+  tau_current = 0.d0
+
+  call random_gauss(r_old,3)
+
+  call drift(a,r_old,d_old)
+  d2_old  = d_old(1)*d_old(1) + &
+            d_old(2)*d_old(2) + &
+            d_old(3)*d_old(3)
+
+  psi_old = psi(a,r_old)
+
+  do istep = 1,nmax
+     e = e_loc(a,r_old)
+     w = w * dexp(-dt*(e - E_ref))
+
+     energy = energy + w*e
+     norm   = norm + w
      
+     tau_current = tau_current + dt
+
+     ! Reset when tau is reached
+     if (tau_current >= tau) then
+        w = 1.d0
+        tau_current = 0.d0
+     endif
+
+     call random_gauss(chi,3)
+     r_new(:) = r_old(:) + dt*d_old(:) + chi(:)*sq_dt
+
+     call drift(a,r_new,d_new)
+     d2_new = d_new(1)*d_new(1) + &
+              d_new(2)*d_new(2) + &
+              d_new(3)*d_new(3)
+
+     psi_new = psi(a,r_new)
+
+     ! Metropolis
+     prod = (d_new(1) + d_old(1))*(r_new(1) - r_old(1)) + &
+            (d_new(2) + d_old(2))*(r_new(2) - r_old(2)) + &
+            (d_new(3) + d_old(3))*(r_new(3) - r_old(3))
+
+     argexpo = 0.5d0 * (d2_new - d2_old)*dt + prod
+
+     q = psi_new / psi_old
+     q = dexp(-argexpo) * q*q
+
+     call random_number(u)
+
+     if (u <= q) then
+
+        n_accep = n_accep + 1_8
+
+        r_old(:) = r_new(:)
+        d_old(:) = d_new(:)
+        d2_old   = d2_new
+        psi_old  = psi_new
+
+     end if
+
+  end do
+
+  energy = energy / norm
+  accep  = dble(n_accep) / dble(nmax)
+
+end subroutine pdmc
+
+program qmc
+  implicit none
+  double precision, parameter :: a     = 0.9
+  double precision, parameter :: dt    = 0.1d0
+  double precision, parameter :: E_ref = -0.5d0
+  double precision, parameter :: tau   = 10.d0
+  integer*8       , parameter :: nmax  = 100000
+  integer         , parameter :: nruns = 30
+
+  integer          :: irun
+  double precision :: X(nruns), accep(nruns)
+  double precision :: ave, err
+
+  do irun=1,nruns
+     call pdmc(a, dt, nmax, X(irun), accep(irun), tau, E_ref)
+  enddo
+
+  call ave_error(X,nruns,ave,err)
+  print *, 'E = ', ave, '+/-', err
+
+  call ave_error(accep,nruns,ave,err)
+  print *, 'A = ', ave, '+/-', err
+
+end program qmc
+     #+END_SRC
+   
+     #+begin_src sh :results output :exports results
+gfortran hydrogen.f90 qmc_stats.f90 pdmc.f90 -o pdmc
+./pdmc
+     #+end_src
+
+     #+RESULTS:
+     :  E =  -0.50067519934141380      +/-   7.9390940184720371E-004
+     :  A =   0.98788066666666663      +/-   7.2889356133441110E-005
+
 
 ** TODO H_2
 
@@ -2185,6 +2558,11 @@ gfortran hydrogen.f90 qmc_stats.f90 vmc_metropolis.f90 -o vmc_metropolis
    
 * Old sections to be removed                                       :noexport:
 
+   :PROPERTIES:
+   :header-args:python: :tangle none
+   :header-args:f90: :tangle none
+   :END:
+   
 ** Gaussian sampling                                               :noexport:
    :PROPERTIES:
    :header-args:python: :tangle qmc_gaussian.py
@@ -2399,11 +2777,11 @@ gfortran hydrogen.f90 qmc_stats.f90 qmc_gaussian.f90 -o qmc_gaussian
     diffusion scheme:
 
     \[
-    \mathbf{r}_{n+1} = \mathbf{r}_{n} + \tau\, 2D \frac{\nabla
+    \mathbf{r}_{n+1} = \mathbf{r}_{n} + \delta t\, 2D \frac{\nabla
     \Psi(\mathbf{r})}{\Psi(\mathbf{r})} + \chi 
     \]
     where $\chi$ is a Gaussian random variable with zero mean and
-    variance $\tau\,2D$.
+    variance $\delta t\,2D$.
    
 **** Exercise 2
 
@@ -2418,8 +2796,8 @@ gfortran hydrogen.f90 qmc_stats.f90 qmc_gaussian.f90 -o qmc_gaussian
 from hydrogen  import *
 from qmc_stats import *
 
-def MonteCarlo(a,tau,nmax):
-    sq_tau = np.sqrt(tau)
+def MonteCarlo(a,dt,nmax):
+    sq_dt = np.sqrt(dt)
 
     # Initialization
     E = 0.
@@ -2429,7 +2807,7 @@ def MonteCarlo(a,tau,nmax):
     for istep in range(nmax):
         d_old = drift(a,r_old)
         chi = np.random.normal(loc=0., scale=1.0, size=(3))
-        r_new = r_old + tau * d_old + chi*sq_tau
+        r_new = r_old + dt * d_old + chi*sq_dt
         N += 1.
         E += e_loc(a,r_new)
         r_old = r_new
@@ -2438,8 +2816,8 @@ def MonteCarlo(a,tau,nmax):
 
 a = 0.9
 nmax = 100000
-tau = 0.2
-X = [MonteCarlo(a,tau,nmax) for i in range(30)]
+dt = 0.2
+X = [MonteCarlo(a,dt,nmax) for i in range(30)]
 E, deltaE = ave_error(X)
 print(f"E = {E} +/- {deltaE}")
       #+END_SRC
@@ -2449,17 +2827,17 @@ print(f"E = {E} +/- {deltaE}")
 
       *Fortran*
       #+BEGIN_SRC f90
-subroutine variational_montecarlo(a,tau,nmax,energy)
+subroutine variational_montecarlo(a,dt,nmax,energy)
   implicit none
-  double precision, intent(in)  :: a, tau
+  double precision, intent(in)  :: a, dt
   integer*8       , intent(in)  :: nmax 
   double precision, intent(out) :: energy
 
   integer*8 :: istep
-  double precision :: norm, r_old(3), r_new(3), d_old(3), sq_tau, chi(3)
+  double precision :: norm, r_old(3), r_new(3), d_old(3), sq_dt, chi(3)
   double precision, external :: e_loc
 
-  sq_tau = dsqrt(tau)
+  sq_dt = dsqrt(dt)
   
   ! Initialization
   energy = 0.d0
@@ -2469,7 +2847,7 @@ subroutine variational_montecarlo(a,tau,nmax,energy)
   do istep = 1,nmax
      call drift(a,r_old,d_old)
      call random_gauss(chi,3)
-     r_new(:) = r_old(:) + tau * d_old(:) + chi(:)*sq_tau
+     r_new(:) = r_old(:) + dt * d_old(:) + chi(:)*sq_dt
      norm = norm + 1.d0
      energy = energy + e_loc(a,r_new)
      r_old(:) = r_new(:)
@@ -2480,7 +2858,7 @@ end subroutine variational_montecarlo
 program qmc
   implicit none
   double precision, parameter :: a = 0.9
-  double precision, parameter :: tau = 0.2
+  double precision, parameter :: dt = 0.2
   integer*8       , parameter :: nmax = 100000
   integer         , parameter :: nruns = 30
 
@@ -2489,7 +2867,7 @@ program qmc
   double precision :: ave, err
 
   do irun=1,nruns
-     call variational_montecarlo(a,tau,nmax,X(irun))
+     call variational_montecarlo(a,dt,nmax,X(irun))
   enddo
   call ave_error(X,nruns,ave,err)
   print *, 'E = ', ave, '+/-', err
diff --git a/energy_hydrogen.f90 b/energy_hydrogen.f90
index 2f65f5d..f875c4b 100644
--- a/energy_hydrogen.f90
+++ b/energy_hydrogen.f90
@@ -12,7 +12,9 @@ program energy_hydrogen
   end do
 
   do j=1,size(a)
+
      ! TODO
+
      print *, 'a = ', a(j), '    E = ', energy
   end do
 
@@ -22,7 +24,7 @@ program energy_hydrogen
   implicit none
   double precision, external :: e_loc, psi
   double precision :: x(50), w, delta, energy, dx, r(3), a(6), norm
-  integer :: i, k, l, j
+  integer          :: i, k, l, j
 
   a = (/ 0.1d0, 0.2d0, 0.5d0, 1.d0, 1.5d0, 2.d0 /)
 
@@ -37,20 +39,28 @@ program energy_hydrogen
 
   do j=1,size(a)
      energy = 0.d0
-     norm = 0.d0
+     norm   = 0.d0
+     
      do i=1,size(x)
         r(1) = x(i)
+
         do k=1,size(x)
            r(2) = x(k)
+
            do l=1,size(x)
               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
+
         end do
+
      end do
+
      energy = energy / norm
      print *, 'a = ', a(j), '    E = ', energy
   end do
diff --git a/energy_hydrogen.py b/energy_hydrogen.py
index c6342b7..14be354 100644
--- a/energy_hydrogen.py
+++ b/energy_hydrogen.py
@@ -7,17 +7,21 @@ delta = (interval[1]-interval[0])**3
 r = np.array([0.,0.,0.])
 
 for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]:
-    E = 0.
+    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 527e690..338f42b 100644
--- a/hydrogen.f90
+++ b/hydrogen.f90
@@ -1,45 +1,62 @@
 double precision function potential(r)
   implicit none
   double precision, intent(in) :: r(3)
-  double precision :: distance
+
+  double precision             :: distance
+
   distance = dsqrt( r(1)*r(1) + r(2)*r(2) + r(3)*r(3) )
+
   if (distance > 0.d0) then
-     potential = -1.d0 / dsqrt( r(1)*r(1) + r(2)*r(2) + r(3)*r(3) )
+     potential = -1.d0 / distance
   else
      stop 'potential at r=0.d0 diverges'
   end if
+
 end function potential
 
 double precision function psi(a, r)
   implicit none
   double precision, intent(in) :: a, r(3)
+
   psi = dexp(-a * dsqrt( r(1)*r(1) + r(2)*r(2) + r(3)*r(3) ))
 end function psi
 
 double precision function kinetic(a,r)
   implicit none
   double precision, intent(in) :: a, r(3)
-  double precision :: distance
+
+  double precision             :: distance
+
   distance = dsqrt( r(1)*r(1) + r(2)*r(2) + r(3)*r(3) ) 
+
   if (distance > 0.d0) then
-     kinetic = -0.5d0 * (a*a - (2.d0*a) / distance)
+
+     kinetic =  a * (1.d0 / distance - 0.5d0 * a)
+
   else
      stop 'kinetic energy diverges at r=0'
   end if
+
 end function kinetic
 
 double precision function e_loc(a,r)
   implicit none
   double precision, intent(in) :: a, r(3)
+
   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
+  b(:)   = r(:) * ar_inv
+
 end subroutine drift
diff --git a/hydrogen.py b/hydrogen.py
index 7b6ea6e..b3746c5 100644
--- a/hydrogen.py
+++ b/hydrogen.py
@@ -11,7 +11,8 @@ def psi(a, r):
 def kinetic(a,r):
     distance = np.sqrt(np.dot(r,r))
     assert (distance > 0.)
-    return -0.5 * (a**2 - (2.*a)/distance)
+
+    return a * (1./distance - 0.5 * a)
 
 def e_loc(a,r):
     return kinetic(a,r) + potential(r)
diff --git a/qmc_metropolis.f90 b/qmc_metropolis.f90
index 41fbd9b..a3db63f 100644
--- a/qmc_metropolis.f90
+++ b/qmc_metropolis.f90
@@ -1,53 +1,65 @@
-subroutine metropolis_montecarlo(a,nmax,tau,energy,accep)
+subroutine metropolis_montecarlo(a,nmax,dt,energy,accep)
   implicit none
   double precision, intent(in)  :: a
   integer*8       , intent(in)  :: nmax 
-  double precision, intent(in)  :: tau
+  double precision, intent(in)  :: dt
   double precision, intent(out) :: energy
   double precision, intent(out) :: accep
 
-  integer*8 :: istep
-
   double precision :: r_old(3), r_new(3), psi_old, psi_new
   double precision :: v, ratio
   integer*8        :: n_accep
+  integer*8        :: istep
+
   double precision, external :: e_loc, psi, gaussian
 
-  energy = 0.d0
+  energy  = 0.d0
   n_accep = 0_8
+
   call random_number(r_old)
-  r_old(:) = tau * (2.d0*r_old(:) - 1.d0)
+  r_old(:) = dt * (2.d0*r_old(:) - 1.d0)
   psi_old = psi(a,r_old)
+
   do istep = 1,nmax
+     energy = energy + e_loc(a,r_old)
+
      call random_number(r_new)
-     r_new(:) = r_old(:) + tau * (2.d0*r_new(:) - 1.d0)
+     r_new(:) = r_old(:) + dt*(2.d0*r_new(:) - 1.d0)
+
      psi_new = psi(a,r_new)
+
      ratio = (psi_new / psi_old)**2
      call random_number(v)
+
      if (v <= ratio) then
+
+        n_accep = n_accep + 1_8
+
         r_old(:) = r_new(:)
         psi_old = psi_new
-        n_accep = n_accep + 1_8
+
      endif
-     energy = energy + e_loc(a,r_old)
+
   end do
+
   energy = energy / dble(nmax)
   accep  = dble(n_accep) / dble(nmax)
+
 end subroutine metropolis_montecarlo
 
 program qmc
   implicit none
   double precision, parameter :: a = 0.9d0
-  double precision, parameter :: tau = 1.3d0
+  double precision, parameter :: dt = 1.3d0
   integer*8       , parameter :: nmax = 100000
   integer         , parameter :: nruns = 30
 
-  integer :: irun
+  integer          :: irun
   double precision :: X(nruns), Y(nruns)
   double precision :: ave, err
 
   do irun=1,nruns
-     call metropolis_montecarlo(a,nmax,tau,X(irun),Y(irun))
+     call metropolis_montecarlo(a,nmax,dt,X(irun),Y(irun))
   enddo
 
   call ave_error(X,nruns,ave,err)
@@ -55,4 +67,5 @@ program qmc
 
   call ave_error(Y,nruns,ave,err)
   print *, 'A = ', ave, '+/-', err
+
 end program qmc
diff --git a/qmc_metropolis.py b/qmc_metropolis.py
index 7f16615..8581d30 100644
--- a/qmc_metropolis.py
+++ b/qmc_metropolis.py
@@ -1,28 +1,35 @@
 from hydrogen  import *
 from qmc_stats import *
 
-def MonteCarlo(a,nmax,tau):
-    energy = 0.
+def MonteCarlo(a,nmax,dt):
+    energy  = 0.
     N_accep = 0
-    r_old = np.random.uniform(-tau, tau, (3))
+
+    r_old = np.random.uniform(-dt, dt, (3))
     psi_old = psi(a,r_old)
+
     for istep in range(nmax):
-        r_new = r_old + np.random.uniform(-tau,tau,(3))
-        psi_new = psi(a,r_new)
-        ratio = (psi_new / psi_old)**2
-        v = np.random.uniform()
-        if v <= ratio:
-            N_accep += 1
-            r_old = r_new
-            psi_old = psi_new
         energy += e_loc(a,r_old)
+
+        r_new = r_old + np.random.uniform(-dt,dt,(3))
+        psi_new = psi(a,r_new)
+
+        ratio = (psi_new / psi_old)**2
+
+        if np.random.uniform() <= ratio:
+            N_accep += 1
+
+            r_old   = r_new
+            psi_old = psi_new
+
     return energy/nmax, N_accep/nmax
 
 # Run simulation
-a = 0.9
+a    = 0.9
 nmax = 100000
-tau = 1.3
-X0 = [ MonteCarlo(a,nmax,tau) for i in range(30)]
+dt   = 1.3
+
+X0 = [ MonteCarlo(a,nmax,dt) for i in range(30)]
 
 # Energy
 X = [ x for (x, _) in X0 ]
diff --git a/qmc_stats.f90 b/qmc_stats.f90
index 3af18d0..6cfef38 100644
--- a/qmc_stats.f90
+++ b/qmc_stats.f90
@@ -1,26 +1,25 @@
 subroutine ave_error(x,n,ave,err)
   implicit none
-  integer, intent(in)           :: n 
-  double precision, intent(in)  :: x(n) 
-  double precision, intent(out) :: ave, err
-  ! TODO
-end subroutine ave_error
 
-subroutine ave_error(x,n,ave,err)
-  implicit none
   integer, intent(in)           :: n 
   double precision, intent(in)  :: x(n) 
   double precision, intent(out) :: ave, err
-  double precision :: variance
+
+  double precision              :: variance
+
   if (n < 1) then
      stop 'n<1 in ave_error'
+
   else if (n == 1) then
      ave = x(1)
      err = 0.d0
+
   else
-     ave = sum(x(:)) / dble(n)
-     variance = sum( (x(:) - ave)**2 ) / dble(n-1)
-     err = dsqrt(variance/dble(n))
+     ave      = sum(x(:)) / dble(n)
+
+     variance = sum((x(:) - ave)**2) / dble(n-1)
+     err      = dsqrt(variance/dble(n))
+
   endif
 end subroutine ave_error
 
@@ -33,6 +32,7 @@ subroutine random_gauss(z,n)
   integer :: i
 
   call random_number(u)
+
   if (iand(n,1) == 0) then
      ! n is even
      do i=1,n,2
@@ -40,6 +40,7 @@ subroutine random_gauss(z,n)
         z(i+1) = z(i) * dsin( two_pi*u(i+1) )
         z(i)   = z(i) * dcos( two_pi*u(i+1) )
      end do
+
   else
      ! n is odd
      do i=1,n-1,2
@@ -47,7 +48,10 @@ subroutine random_gauss(z,n)
         z(i+1) = z(i) * dsin( two_pi*u(i+1) )
         z(i)   = z(i) * dcos( two_pi*u(i+1) )
      end do
+
      z(n)   = dsqrt(-2.d0*dlog(u(n))) 
      z(n)   = z(n) * dcos( two_pi*u(n+1) )
+
   end if
+
 end subroutine random_gauss
diff --git a/qmc_stats.py b/qmc_stats.py
index 5b75731..e21ebb9 100644
--- a/qmc_stats.py
+++ b/qmc_stats.py
@@ -2,10 +2,14 @@ from math import sqrt
 def ave_error(arr):
     M = len(arr)
     assert(M>0)
+
     if M == 1:
-        return (arr[0], 0.)
+        average = arr[0]
+        error   = 0.
+
     else:
         average = sum(arr)/M
         variance = 1./(M-1) * sum( [ (x - average)**2 for x in arr ] )
         error = sqrt(variance/M)
-        return (average, error)
+
+    return (average, error)
diff --git a/qmc_uniform.f90 b/qmc_uniform.f90
index a67fc70..ec95768 100644
--- a/qmc_uniform.f90
+++ b/qmc_uniform.f90
@@ -4,38 +4,46 @@ subroutine uniform_montecarlo(a,nmax,energy)
   integer*8       , intent(in)  :: nmax 
   double precision, intent(out) :: energy
 
-  integer*8 :: istep
-
+  integer*8        :: istep
   double precision :: norm, r(3), w
 
   double precision, external :: e_loc, psi
 
   energy = 0.d0
   norm   = 0.d0
+
   do istep = 1,nmax
+
      call random_number(r)
      r(:) = -5.d0 + 10.d0*r(:)
+
      w = psi(a,r)
      w = w*w
-     norm = norm + w
+
      energy = energy + w * e_loc(a,r)
+     norm   = norm   + w
+
   end do
+
   energy = energy / norm
+
 end subroutine uniform_montecarlo
 
 program qmc
   implicit none
-  double precision, parameter :: a = 0.9
-  integer*8       , parameter :: nmax = 100000
+  double precision, parameter :: a     = 0.9
+  integer*8       , parameter :: nmax  = 100000
   integer         , parameter :: nruns = 30
 
-  integer :: irun
+  integer          :: irun
   double precision :: X(nruns)
   double precision :: ave, err
 
   do irun=1,nruns
-     call uniform_montecarlo(a,nmax,X(irun))
+     call uniform_montecarlo(a, nmax, X(irun))
   enddo
-  call ave_error(X,nruns,ave,err)
+
+  call ave_error(X, nruns, ave, err)
+
   print *, 'E = ', ave, '+/-', err
 end program qmc
diff --git a/qmc_uniform.py b/qmc_uniform.py
index 516d196..66a8616 100644
--- a/qmc_uniform.py
+++ b/qmc_uniform.py
@@ -4,16 +4,22 @@ from qmc_stats import *
 def MonteCarlo(a, nmax):
      energy = 0.
      normalization = 0.
+
      for istep in range(nmax):
           r = np.random.uniform(-5., 5., (3))
+
           w = psi(a,r)
           w = w*w
-          normalization += w
-          energy += w * e_loc(a,r)
-     return energy/normalization
 
-a = 0.9
+          energy        += w * e_loc(a,r)
+          normalization += w
+
+     return energy / normalization
+
+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 9f53dfc..80d0b03 100644
--- a/variance_hydrogen.f90
+++ b/variance_hydrogen.f90
@@ -1,9 +1,11 @@
 program variance_hydrogen
   implicit none
+
+  double precision :: x(50), w, delta, energy, energy2
+  double precision :: dx, r(3), a(6), norm, e_tmp, s2
+  integer          :: i, k, l, j
+
   double precision, external :: e_loc, psi
-  double precision :: x(50), w, delta, energy, dx, r(3), a(6), norm, s2
-  double precision :: e, energy2
-  integer :: i, k, l, j
 
   a = (/ 0.1d0, 0.2d0, 0.5d0, 1.d0, 1.5d0, 2.d0 /)
 
@@ -17,27 +19,38 @@ program variance_hydrogen
   r(:) = 0.d0
 
   do j=1,size(a)
-     energy = 0.d0
+     energy  = 0.d0
      energy2 = 0.d0
-     norm = 0.d0
+     norm    = 0.d0
+
      do i=1,size(x)
         r(1) = x(i)
+
         do k=1,size(x)
            r(2) = x(k)
+
            do l=1,size(x)
               r(3) = x(l)
+
               w = psi(a(j),r)
               w = w * w * delta
-              e = e_loc(a(j), r)
-              energy  = energy  + w * e
-              energy2 = energy2 + w * e * e
-              norm   = norm   + w 
+
+              e_tmp = e_loc(a(j), r)
+
+              energy  = energy  + w * e_tmp
+              energy2 = energy2 + w * e_tmp * e_tmp
+              norm   = norm     + w 
            end do
+
         end do
+
      end do
+
      energy  = energy  / norm
      energy2 = energy2 / norm
+
      s2 = energy2 - energy*energy
+
      print *, 'a = ', a(j), ' E = ', energy, ' s2 = ', s2
   end do
 
diff --git a/variance_hydrogen.py b/variance_hydrogen.py
index 40f9a27..abaf3bd 100644
--- a/variance_hydrogen.py
+++ b/variance_hydrogen.py
@@ -2,27 +2,35 @@ import numpy as np
 from hydrogen import e_loc, psi
 
 interval = np.linspace(-5,5,num=50)
+
 delta = (interval[1]-interval[0])**3
 
 r = np.array([0.,0.,0.])
 
 for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]:
-    E = 0.
-    E2 = 0.
+    E    = 0.
+    E2   = 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 = psi(a,r)
                 w = w * w * delta
-                El = e_loc(a, r)
-                E  += w * El
-                E2 += w * El*El
+
+                e_tmp = e_loc(a,r)
+                E    += w * e_tmp
+                E2   += w * e_tmp * e_tmp
                 norm += w 
-    E = E / norm
+
+    E  = E  / norm
     E2 = E2 / norm
-    s2 = E2 - E*E
-    print(f"a = {a} \t E = {E:10.8f}  \t  \sigma^2 = {s2:10.8f}")
+
+    s2 = E2 - E**2
+    print(f"a = {a} \t E = {E:10.8f} \t \sigma^2 = {s2:10.8f}")
diff --git a/vmc.f90 b/vmc.f90
index 42ee411..6a3dde8 100644
--- a/vmc.f90
+++ b/vmc.f90
@@ -1,14 +1,14 @@
-subroutine variational_montecarlo(a,tau,nmax,energy)
+subroutine variational_montecarlo(a,dt,nmax,energy)
   implicit none
-  double precision, intent(in)  :: a, tau
+  double precision, intent(in)  :: a, dt
   integer*8       , intent(in)  :: nmax 
   double precision, intent(out) :: energy
 
   integer*8 :: istep
-  double precision :: norm, r_old(3), r_new(3), d_old(3), sq_tau, chi(3)
+  double precision :: norm, r_old(3), r_new(3), d_old(3), sq_dt, chi(3)
   double precision, external :: e_loc
 
-  sq_tau = dsqrt(tau)
+  sq_dt = dsqrt(dt)
   
   ! Initialization
   energy = 0.d0
@@ -18,7 +18,7 @@ subroutine variational_montecarlo(a,tau,nmax,energy)
   do istep = 1,nmax
      call drift(a,r_old,d_old)
      call random_gauss(chi,3)
-     r_new(:) = r_old(:) + tau * d_old(:) + chi(:)*sq_tau
+     r_new(:) = r_old(:) + dt * d_old(:) + chi(:)*sq_dt
      norm = norm + 1.d0
      energy = energy + e_loc(a,r_new)
      r_old(:) = r_new(:)
@@ -29,7 +29,7 @@ end subroutine variational_montecarlo
 program qmc
   implicit none
   double precision, parameter :: a = 0.9
-  double precision, parameter :: tau = 0.2
+  double precision, parameter :: dt = 0.2
   integer*8       , parameter :: nmax = 100000
   integer         , parameter :: nruns = 30
 
@@ -38,7 +38,7 @@ program qmc
   double precision :: ave, err
 
   do irun=1,nruns
-     call variational_montecarlo(a,tau,nmax,X(irun))
+     call variational_montecarlo(a,dt,nmax,X(irun))
   enddo
   call ave_error(X,nruns,ave,err)
   print *, 'E = ', ave, '+/-', err
diff --git a/vmc.py b/vmc.py
index ca422be..b34321e 100644
--- a/vmc.py
+++ b/vmc.py
@@ -1,8 +1,8 @@
 from hydrogen  import *
 from qmc_stats import *
 
-def MonteCarlo(a,tau,nmax):
-    sq_tau = np.sqrt(tau)
+def MonteCarlo(a,dt,nmax):
+    sq_dt = np.sqrt(dt)
 
     # Initialization
     E = 0.
@@ -12,7 +12,7 @@ def MonteCarlo(a,tau,nmax):
     for istep in range(nmax):
         d_old = drift(a,r_old)
         chi = np.random.normal(loc=0., scale=1.0, size=(3))
-        r_new = r_old + tau * d_old + chi*sq_tau
+        r_new = r_old + dt * d_old + chi*sq_dt
         N += 1.
         E += e_loc(a,r_new)
         r_old = r_new
@@ -21,7 +21,7 @@ def MonteCarlo(a,tau,nmax):
 
 a = 0.9
 nmax = 100000
-tau = 0.2
-X = [MonteCarlo(a,tau,nmax) for i in range(30)]
+dt = 0.2
+X = [MonteCarlo(a,dt,nmax) for i in range(30)]
 E, deltaE = ave_error(X)
 print(f"E = {E} +/- {deltaE}")
diff --git a/vmc_metropolis.f90 b/vmc_metropolis.f90
index 3202e1a..499f09b 100644
--- a/vmc_metropolis.f90
+++ b/vmc_metropolis.f90
@@ -1,69 +1,95 @@
-subroutine variational_montecarlo(a,tau,nmax,energy,accep_rate)
+subroutine variational_montecarlo(a,dt,nmax,energy,accep)
   implicit none
-  double precision, intent(in)  :: a, tau
+  double precision, intent(in)  :: a, dt
   integer*8       , intent(in)  :: nmax 
-  double precision, intent(out) :: energy, accep_rate
+  double precision, intent(out) :: energy, accep
 
-  integer*8 :: istep
-  double precision :: sq_tau, chi(3), d2_old, prod, u
+  integer*8        :: istep
+  integer*8        :: n_accep
+  double precision :: sq_dt, chi(3), d2_old, prod, u
   double precision :: psi_old, psi_new, d2_new, argexpo, q
   double precision :: r_old(3), r_new(3)
   double precision :: d_old(3), d_new(3)
+
   double precision, external :: e_loc, psi
 
-  sq_tau = dsqrt(tau)
+  sq_dt = dsqrt(dt)
 
   ! Initialization
-  energy = 0.d0
-  accep_rate = 0.d0
+  energy  = 0.d0
+  n_accep = 0_8
+
   call random_gauss(r_old,3)
+
   call drift(a,r_old,d_old)
-  d2_old = d_old(1)*d_old(1) + d_old(2)*d_old(2) + d_old(3)*d_old(3)
+  d2_old  = d_old(1)*d_old(1) + &
+            d_old(2)*d_old(2) + &
+            d_old(3)*d_old(3)
+
   psi_old = psi(a,r_old)
 
   do istep = 1,nmax
+     energy = energy + e_loc(a,r_old)
+
      call random_gauss(chi,3)
-     r_new(:) = r_old(:) + tau * d_old(:) + chi(:)*sq_tau
+     r_new(:) = r_old(:) + dt*d_old(:) + chi(:)*sq_dt
+
      call drift(a,r_new,d_new)
-     d2_new = d_new(1)*d_new(1) + d_new(2)*d_new(2) + d_new(3)*d_new(3)
+     d2_new = d_new(1)*d_new(1) + &
+              d_new(2)*d_new(2) + &
+              d_new(3)*d_new(3)
+
      psi_new = psi(a,r_new)
+
      ! Metropolis
      prod = (d_new(1) + d_old(1))*(r_new(1) - r_old(1)) + &
             (d_new(2) + d_old(2))*(r_new(2) - r_old(2)) + &
             (d_new(3) + d_old(3))*(r_new(3) - r_old(3))
-     argexpo = 0.5d0 * (d2_new - d2_old)*tau + prod
+
+     argexpo = 0.5d0 * (d2_new - d2_old)*dt + prod
+
      q = psi_new / psi_old
      q = dexp(-argexpo) * q*q
+
      call random_number(u)
-     if (u<q) then
-        accep_rate = accep_rate + 1.d0
+
+     if (u <= q) then
+
+        n_accep = n_accep + 1_8
+
         r_old(:) = r_new(:)
         d_old(:) = d_new(:)
-        d2_old = d2_new
-        psi_old = psi_new
+        d2_old   = d2_new
+        psi_old  = psi_new
+
      end if
-     energy = energy + e_loc(a,r_old)
+
   end do
+
   energy = energy / dble(nmax)
-  accep_rate = dble(accep_rate) / dble(nmax)
+  accep  = dble(n_accep) / dble(nmax)
+
 end subroutine variational_montecarlo
 
 program qmc
   implicit none
-  double precision, parameter :: a = 0.9
-  double precision, parameter :: tau = 1.0
-  integer*8       , parameter :: nmax = 100000
+  double precision, parameter :: a     = 0.9
+  double precision, parameter :: dt    = 1.0
+  integer*8       , parameter :: nmax  = 100000
   integer         , parameter :: nruns = 30
 
-  integer :: irun
+  integer          :: irun
   double precision :: X(nruns), accep(nruns)
   double precision :: ave, err
 
   do irun=1,nruns
-     call variational_montecarlo(a,tau,nmax,X(irun),accep(irun))
+     call variational_montecarlo(a,dt,nmax,X(irun),accep(irun))
   enddo
+
   call ave_error(X,nruns,ave,err)
   print *, 'E = ', ave, '+/-', err
+
   call ave_error(accep,nruns,ave,err)
   print *, 'A = ', ave, '+/-', err
+
 end program qmc
diff --git a/vmc_metropolis.py b/vmc_metropolis.py
index 716b828..4181fe4 100644
--- a/vmc_metropolis.py
+++ b/vmc_metropolis.py
@@ -1,40 +1,51 @@
 from hydrogen  import *
 from qmc_stats import *
 
-def MonteCarlo(a,nmax,tau):
-    E = 0.
-    accep_rate = 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)
+def MonteCarlo(a,nmax,dt):
+    sq_dt = np.sqrt(dt)
+
+    energy  = 0.
+    N_accep = 0
+
+    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):
         chi = np.random.normal(loc=0., scale=1.0, size=(3))
-        r_new = r_old + tau * d_old + sq_tau * chi
-        d_new = drift(a,r_new)
-        d2_new = np.dot(d_new,d_new)
+
+        energy += e_loc(a,r_old)
+
+        r_new   = r_old + dt * d_old + sq_dt * chi
+        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
+        prod    = np.dot((d_new + d_old), (r_new - r_old))
+        argexpo = 0.5 * (d2_new - d2_old)*dt + prod
+
         q = psi_new / psi_old
         q = np.exp(-argexpo) * q*q
-        if np.random.uniform() < q:
-            accep_rate += 1.
-            r_old = r_new
-            d_old = d_new
-            d2_old = d2_new
+
+        if np.random.uniform() <= q:
+            N_accep += 1
+
+            r_old   = r_new
+            d_old   = d_new
+            d2_old  = d2_new
             psi_old = psi_new
-        E += e_loc(a,r_old)
-    return E/nmax, accep_rate/nmax
+
+    return energy/nmax, accep_rate/nmax
 
 
 # Run simulation
-a = 0.9
+a    = 0.9
 nmax = 100000
-tau = 1.3
-X0 = [ MonteCarlo(a,nmax,tau) for i in range(30)]
+dt   = 1.3
+
+X0 = [ MonteCarlo(a,nmax,dt) for i in range(30)]
 
 # Energy
 X = [ x for (x, _) in X0 ]