From d50db6bf0cd9cf3529cf6fa2dc2d99aba7f55145 Mon Sep 17 00:00:00 2001 From: scemama Date: Thu, 4 Feb 2021 16:02:43 +0000 Subject: [PATCH] deploy: efb48d936f29665cce00393a9a2bbde02b230f72 --- index.html | 1000 ++++++++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 846 insertions(+), 154 deletions(-) diff --git a/index.html b/index.html index 6166797..edeac2e 100644 --- a/index.html +++ b/index.html @@ -3,7 +3,7 @@ "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd"> - + Quantum Monte Carlo @@ -329,124 +329,148 @@ for the JavaScript code in this tag.

Table of Contents

+ + + +
  • 3. Variational Monte Carlo -
    • -
  • 3.2. Uniform sampling in the box +
  • 3.1. Computation of the statistical error -
    • -
  • 3.3. Metropolis sampling with \(\Psi^2\) +
  • 3.1.1. Exercise
    • -
  • 3.4. Generalized Metropolis algorithm + +
    • +
  • 3.2. Uniform sampling in the box -
    • - - -
  • 4. Diffusion Monte Carlo +
  • 3.2.1. Exercise +
    • + + +
  • 3.3. Metropolis sampling with \(\Psi^2\) -
    • -
  • 4.4. Pure Diffusion Monte Carlo
    • -
  • 4.5. Hydrogen atom +
  • 3.3.1. Optimal step size
    • +
  • 3.3.2. Exercise
    • -
  • 5. Project
    • -
  • 6. Acknowledgments
    • +
  • 3.4. Generalized Metropolis algorithm + +
    • + + +
  • 4. Diffusion Monte Carlo + +
    • +
  • 5. Project
    • +
  • 6. Acknowledgments
    • -
    -

    1 Introduction

    +
    +

    1 Introduction

    This website contains the QMC tutorial of the 2021 LTTC winter school @@ -486,8 +510,8 @@ coordinates, etc).

    -
    -

    1.1 Energy and local energy

    +
    +

    1.1 Energy and local energy

    For a given system with Hamiltonian \(\hat{H}\) and wave function \(\Psi\), we define the local energy as @@ -570,8 +594,8 @@ energy computed over these configurations:

    -
    -

    2 Numerical evaluation of the energy of the hydrogen atom

    +
    +

    2 Numerical evaluation of the energy of the hydrogen atom

    In this section, we consider the hydrogen atom with the following @@ -600,8 +624,8 @@ To do that, we will compute the local energy and check whether it is constant.

    -
    -

    2.1 Local energy

    +
    +

    2.1 Local energy

    You will now program all quantities needed to compute the local energy of the H atom for the given wave function. @@ -628,8 +652,8 @@ to catch the error.

    -
    -

    2.1.1 Exercise 1

    +
    +

    2.1.1 Exercise 1

    @@ -674,8 +698,8 @@ and returns the potential.

    -
    -
    2.1.1.1 Solution   solution2
    +
    +
    2.1.1.1 Solution   solution2

    Python @@ -716,8 +740,8 @@ and returns the potential.

    -
    -

    2.1.2 Exercise 2

    +
    +

    2.1.2 Exercise 2

    @@ -752,8 +776,8 @@ input arguments, and returns a scalar.

    -
    -
    2.1.2.1 Solution   solution2
    +
    +
    2.1.2.1 Solution   solution2

    Python @@ -780,8 +804,8 @@ input arguments, and returns a scalar.

    -
    -

    2.1.3 Exercise 3

    +
    +

    2.1.3 Exercise 3

    @@ -862,8 +886,8 @@ Therefore, the local kinetic energy is

    -
    -
    2.1.3.1 Solution   solution2
    +
    +
    2.1.3.1 Solution   solution2

    Python @@ -904,8 +928,8 @@ Therefore, the local kinetic energy is

    -
    -

    2.1.4 Exercise 4

    +
    +

    2.1.4 Exercise 4

    @@ -964,8 +988,8 @@ are calling is yours.

    -
    -
    2.1.4.1 Solution   solution2
    +
    +
    2.1.4.1 Solution   solution2

    Python @@ -996,8 +1020,8 @@ are calling is yours.

    -
    -

    2.1.5 Exercise 5

    +
    +

    2.1.5 Exercise 5

    @@ -1007,8 +1031,8 @@ Find the theoretical value of \(a\) for which \(\Psi\) is an eigenfunction of \(

    -
    -
    2.1.5.1 Solution   solution2
    +
    +
    2.1.5.1 Solution   solution2
    \begin{eqnarray*} E &=& \frac{\hat{H} \Psi}{\Psi} = - \frac{1}{2} \frac{\Delta \Psi}{\Psi} - @@ -1028,8 +1052,8 @@ equal to -0.5 atomic units.
    -
    -

    2.2 Plot of the local energy along the \(x\) axis

    +
    +

    2.2 Plot of the local energy along the \(x\) axis

    The program you will write in this section will be written in @@ -1060,8 +1084,8 @@ In Fortran, you will need to compile all the source files together:

    -
    -

    2.2.1 Exercise

    +
    +

    2.2.1 Exercise

    @@ -1155,8 +1179,8 @@ plot './data' index 0 using 1:2 with lines title 'a=0.1', \

    -
    -
    2.2.1.1 Solution   solution2
    +
    +
    2.2.1.1 Solution   solution2

    Python @@ -1233,8 +1257,8 @@ plt.savefig("plot_py.png")

    -
    -

    2.3 Numerical estimation of the energy

    +
    +

    2.3 Numerical estimation of the energy

    If the space is discretized in small volume elements \(\mathbf{r}_i\) @@ -1264,8 +1288,8 @@ The energy is biased because:

    -
    -

    2.3.1 Exercise

    +
    +

    2.3.1 Exercise

    @@ -1336,8 +1360,8 @@ To compile the Fortran and run it:

    -
    -
    2.3.1.1 Solution   solution2
    +
    +
    2.3.1.1 Solution   solution2

    Python @@ -1454,8 +1478,8 @@ a = 2.0000000000000000 E = -8.0869806678448772E-002

    -
    -

    2.4 Variance of the local energy

    +
    +

    2.4 Variance of the local energy

    The variance of the local energy is a functional of \(\Psi\) @@ -1482,8 +1506,8 @@ energy can be used as a measure of the quality of a wave function.

    -
    -

    2.4.1 Exercise (optional)

    +
    +

    2.4.1 Exercise (optional)

    @@ -1494,8 +1518,8 @@ Prove that :

    -
    -
    2.4.1.1 DONE Solution   solution2
    +
    +
    2.4.1.1 DONE Solution   solution2

    \(\bar{E} = \langle E \rangle\) is a constant, so \(\langle \bar{E} @@ -1514,8 +1538,8 @@ Prove that :

    -
    -

    2.4.2 Exercise

    +
    +

    2.4.2 Exercise

    @@ -1587,6 +1611,143 @@ To compile and run: 

    gfortran hydrogen.f90 variance_hydrogen.f90 -o variance_hydrogen
 ./variance_hydrogen
+
    +
    +
    + +
    +
    2.4.2.1 Solution   solution2
    +
    +

    +Python +

    +
    +
    #!/usr/bin/env python3
+
+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.
+    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_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}")
+
+
    +
    + +
    +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 +

    + +
    +
    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)
+     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_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
+
    +
    + +
    +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     
+
    @@ -1594,8 +1755,8 @@ To compile and run:
    -
    -

    3 Variational Monte Carlo

    +
    +

    3 Variational Monte Carlo

    Numerical integration with deterministic methods is very efficient @@ -1611,8 +1772,8 @@ interval.

    -
    -

    3.1 Computation of the statistical error

    +
    +

    3.1 Computation of the statistical error

    To compute the statistical error, you need to perform \(M\) @@ -1652,8 +1813,8 @@ And the confidence interval is given by

    -
    -

    3.1.1 Exercise

    +
    +

    3.1.1 Exercise

    @@ -1692,11 +1853,71 @@ input array.

    + +
    +
    3.1.1.1 Solution   solution2
    +
    +

    +Python +

    +
    +
    #!/usr/bin/env python3
+
+from math import sqrt
+def ave_error(arr):
+    M = len(arr)
+    assert(M>0)
+
+    if M == 1:
+        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)
+
    +
    + +

    +Fortran +

    +
    +
    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
+
+  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))
+
+  endif
+end subroutine ave_error
+
    +
    +
    +
    -
    -

    3.2 Uniform sampling in the box

    +
    +

    3.2 Uniform sampling in the box

    We will now perform our first Monte Carlo calculation to compute the @@ -1757,8 +1978,8 @@ compute the statistical error.

    -
    -

    3.2.1 Exercise

    +
    +

    3.2.1 Exercise

    @@ -1856,14 +2077,119 @@ well as the index of the current step. 

    gfortran hydrogen.f90 qmc_stats.f90 qmc_uniform.f90 -o qmc_uniform
 ./qmc_uniform
+
    +
    +
    + +
    +
    3.2.1.1 Solution   solution2
    +
    +

    +Python +

    +
    +
    #!/usr/bin/env python3
+
+from hydrogen  import *
+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
+
+          energy        += w * e_loc(a,r)
+          normalization += w
+
+     return energy / normalization
+
+a    = 1.2
+nmax = 100000
+
+X = [MonteCarlo(a,nmax) for i in range(30)]
+E, deltaE = ave_error(X)
+
+print(f"E = {E} +/- {deltaE}")
+
    +
    + +
    +E = -0.48363807880008725 +/- 0.002330876047368999
+
+
    + +

    +Fortran +

    +
    +
    subroutine uniform_montecarlo(a,nmax,energy)
+  implicit none
+  double precision, intent(in)  :: a
+  integer*8       , intent(in)  :: nmax 
+  double precision, intent(out) :: energy
+
+  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
+
+     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     = 1.2d0
+  integer*8       , parameter :: nmax  = 100000
+  integer         , parameter :: nruns = 30
+
+  integer          :: irun
+  double precision :: X(nruns)
+  double precision :: ave, err
+
+  do irun=1,nruns
+     call uniform_montecarlo(a, nmax, X(irun))
+  enddo
+
+  call ave_error(X, nruns, ave, err)
+
+  print *, 'E = ', ave, '+/-', err
+end program qmc
+
    +
    + +
    +E =  -0.48084122147238995      +/-   2.4983775878329355E-003
+
    -
    -

    3.3 Metropolis sampling with \(\Psi^2\)

    +
    +

    3.3 Metropolis sampling with \(\Psi^2\)

    We will now use the square of the wave function to sample random @@ -1982,8 +2308,8 @@ All samples should be kept, from both accepted and rejected moves.

    -
    -

    3.3.1 Optimal step size

    +
    +

    3.3.1 Optimal step size

    If the box is infinitely small, the ratio will be very close @@ -2018,8 +2344,8 @@ the same variable later on to store a time step.

    -
    -

    3.3.2 Exercise

    +
    +

    3.3.2 Exercise

    @@ -2124,14 +2450,160 @@ Can you observe a reduction in the statistical error? 

    gfortran hydrogen.f90 qmc_stats.f90 qmc_metropolis.f90 -o qmc_metropolis
 ./qmc_metropolis
+
    +
    +
    + +
    +
    3.3.2.1 Solution   solution2
    +
    +

    +Python +

    +
    +
    #!/usr/bin/env python3
+
+from hydrogen  import *
+from qmc_stats import *
+
+def MonteCarlo(a,nmax,dt):
+    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
+
+        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    = 1.2
+nmax = 100000
+dt   = 1.0
+
+X0 = [ MonteCarlo(a,nmax,dt) 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}")
+
    +
    + +
    +E = -0.4802595860693983 +/- 0.0005124420418289207
+A = 0.5074913333333334 +/- 0.000350889422714878
+
+
    + +

    +Fortran +

    +
    +
    subroutine metropolis_montecarlo(a,nmax,dt,energy,accep)
+  implicit none
+  double precision, intent(in)  :: a
+  integer*8       , intent(in)  :: nmax 
+  double precision, intent(in)  :: dt
+  double precision, intent(out) :: energy
+  double precision, intent(out) :: accep
+
+  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
+  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)
+
+     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
+
+     endif
+
+  end do
+
+  energy = energy / dble(nmax)
+  accep  = dble(n_accep) / dble(nmax)
+
+end subroutine metropolis_montecarlo
+
+program qmc
+  implicit none
+  double precision, parameter :: a = 1.2d0
+  double precision, parameter :: dt = 1.0d0
+  integer*8       , parameter :: nmax = 100000
+  integer         , parameter :: nruns = 30
+
+  integer          :: irun
+  double precision :: X(nruns), Y(nruns)
+  double precision :: ave, err
+
+  do irun=1,nruns
+     call metropolis_montecarlo(a,nmax,dt,X(irun),Y(irun))
+  enddo
+
+  call ave_error(X,nruns,ave,err)
+  print *, 'E = ', ave, '+/-', err
+
+  call ave_error(Y,nruns,ave,err)
+  print *, 'A = ', ave, '+/-', err
+
+end program qmc
+
    +
    + +
    +E =  -0.47948142754168033      +/-   4.8410177863919307E-004
+A =   0.50762633333333318      +/-   3.4601756760043725E-004
+
    -
    -

    3.4 Generalized Metropolis algorithm

    +
    +

    3.4 Generalized Metropolis algorithm

    One can use more efficient numerical schemes to move the electrons by choosing a smarter expression for the transition probability. @@ -2252,8 +2724,8 @@ The algorithm of the previous exercise is only slighlty modified as:

    -
    -

    3.4.1 Gaussian random number generator

    +
    +

    3.4.1 Gaussian random number generator

    To obtain Gaussian-distributed random numbers, you can apply the @@ -2317,8 +2789,8 @@ In Python, you can use the -

    3.4.2 Exercise 1

    +
    +

    3.4.2 Exercise 1

    @@ -2359,10 +2831,43 @@ Write a function to compute the drift vector \(\frac{\nabla \Psi(\mathbf{r})}{\P

    + +
    +
    3.4.2.1 Solution   solution2
    +
    +

    +Python +

    +
    +
    def drift(a,r):
+   ar_inv = -a/np.sqrt(np.dot(r,r))
+   return r * ar_inv
+
    -
    -

    3.4.3 Exercise 2

    +

    +Fortran +

    +
    +
    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
+
    +
    +
    +
    +
    + +
    +

    3.4.3 Exercise 2

    @@ -2454,6 +2959,192 @@ Modify the previous program to introduce the drift-diffusion scheme. 

    gfortran hydrogen.f90 qmc_stats.f90 vmc_metropolis.f90 -o vmc_metropolis
 ./vmc_metropolis
+
    +
    +
    + +
    +
    3.4.3.1 Solution   solution2
    +
    +

    +Python +

    +
    +
    #!/usr/bin/env python3
+
+from hydrogen  import *
+from qmc_stats import *
+
+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))
+
+        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)*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/nmax, N_accep/nmax
+
+
+# Run simulation
+a    = 1.2
+nmax = 100000
+dt   = 1.0
+
+X0 = [ MonteCarlo(a,nmax,dt) 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}")
+
    +
    + +
    +E = -0.48034171558629885 +/- 0.0005286038561061781
+A = 0.6210380000000001 +/- 0.0005457375900937905
+
+
    + +

    +Fortran +

    +
    +
    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
+
+  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
+  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)
+
+  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
+
+     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 / dble(nmax)
+  accep  = dble(n_accep) / dble(nmax)
+
+end subroutine variational_montecarlo
+
+program qmc
+  implicit none
+  double precision, parameter :: a     = 1.2d0
+  double precision, parameter :: dt    = 1.0d0
+  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 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
+
    +
    + +
    +E =  -0.47940635575542706      +/-   5.5613594433433764E-004
+A =   0.62037333333333333      +/-   4.8970160591451110E-004
+
    @@ -2461,8 +3152,8 @@ Modify the previous program to introduce the drift-diffusion scheme.
    -
    -

    4 Diffusion Monte Carlo

    +
    +

    4 Diffusion Monte Carlo

    As we have seen, Variational Monte Carlo is a powerful method to @@ -2479,8 +3170,8 @@ finding a near-exact numerical solution to the Schrödinger equation.

    -