Uniform sampling

QMC.org

@@ -1,9 +1,7 @@
 #+TITLE: Quantum Monte Carlo
 #+AUTHOR: Anthony Scemama, Claudia Filippi
 #+SETUPFILE: https://fniessen.github.io/org-html-themes/org/theme-readtheorg.setup
-
 #+STARTUP: latexpreview
-#+STARTUP: indent
 
 
 * Introduction
@@ -62,7 +60,7 @@ can be chosen.
 :header-args:f90: :tangle hydrogen.f90
 :END:
 *** Write a function which computes the potential at $\mathbf{r}$
-   The function accepts q 3-dimensional vector =r= as input arguments
+The function accepts a 3-dimensional vector =r= as input arguments
 and returns the potential.
 
 $\mathbf{r}=\sqrt{x^2 + y^2 + z^2})$, so
@@ -70,13 +68,14 @@ can be chosen.
 V(x,y,z) = -\frac{1}{\sqrt{x^2 + y^2 + z^2})$
 $$
 
-   #+BEGIN_SRC python 
+#+BEGIN_SRC python :results none
 import numpy as np
 
 def potential(r):
     return -1. / np.sqrt(np.dot(r,r))
 #+END_SRC
 
+
 #+BEGIN_SRC f90 
 double precision function potential(r)
   implicit none
@@ -89,7 +88,7 @@ end function potential
 The function accepts a scalar =a= and a 3-dimensional vector =r= as
 input arguments, and returns a scalar.
 
-   #+BEGIN_SRC python 
+#+BEGIN_SRC python :results none
 def psi(a, r):
     return np.exp(-a*np.sqrt(np.dot(r,r)))
 #+END_SRC
@@ -140,7 +139,7 @@ end function psi
 -\frac{1}{2} \frac{\Delta \Psi}{\Psi} (x,y,z) = -\frac{1}{2}\left(a^2 - \frac{2a}{\mathbf{|r|}} \right) 
 $$
      
-   #+BEGIN_SRC python 
+#+BEGIN_SRC python :results none
 def kinetic(a,r):
     return -0.5 * (a**2 - (2.*a)/np.sqrt(np.dot(r,r)))
 #+END_SRC
@@ -162,7 +161,7 @@ end function kinetic
 E_L(x,y,z) = -\frac{1}{2} \frac{\Delta \Psi}{\Psi} (x,y,z) + V(x,y,z)
 $$
 
-   #+BEGIN_SRC python
+#+BEGIN_SRC python :results none
 def e_loc(a,r):
     return kinetic(a,r) + potential(r)
 #+END_SRC
@@ -181,14 +180,11 @@ end function e_loc
 :header-args:python: :tangle plot_hydrogen.py
 :header-args:f90: :tangle plot_hydrogen.f90
 :END:
-:LOGBOOK:
-CLOCK: [2021-01-03 Sun 17:48]
-:END:
 
 For multiple values of $a$ (0.1, 0.2, 0.5, 1., 1.5, 2.), plot the
 local energy along the $x$ axis.
 
-#+begin_src python 
+#+begin_src python :results output
 import numpy as np
 import matplotlib.pyplot as plt
 
@@ -212,6 +208,8 @@ plt.legend()
 plt.savefig("plot_py.png")
 #+end_src
 
+#+RESULTS:
+
 [[./plot_py.png]]
 
 
@@ -286,15 +284,14 @@ plot './data' index 0 using 1:2 with lines title 'a=0.1', \
   & = & \frac{\int \left[\Psi(\mathbf{r})\right]^2\, \frac{\hat{H} \Psi(\mathbf{r})}{\Psi(\mathbf{r})}\,d\mathbf{r}}{\int \left[\Psi(\mathbf{r}) \right]^2 d\mathbf{r}}
 \end{eqnarray}
 
-   If the space is discretized in small volume elements
-   $\delta x\, \delta y\, \delta z$, this last equation corresponds
-   to a weighted average of the local energy, where the weights are
-   the values of the square of the wave function at $(x,y,z)$
-   multiplied by the volume element:
+If the space is discretized in small volume elements $\delta
+\mathbf{r}$, this last equation corresponds to a weighted average of
+the local energy, where the weights are the values of the square of
+the wave function at $\mathbf{r}$ multiplied by the volume element:
      
 $$
 E \approx \frac{\sum_i w_i E_L(\mathbf{r}_i)}{\sum_i w_i}, \;\;
-   w_i = \left[\Psi(\mathbf{r}_i)\right]^2 \delta x\, \delta y\, \delta z
+w_i = \left[\Psi(\mathbf{r}_i)\right]^2 \delta \mathbf{r}
 $$
      
 We now compute an numerical estimate of the energy in a grid of
@@ -340,6 +337,7 @@ for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]:
     : a = 1.5 	 E = -0.39242967082602226
     : a = 2.0 	 E = -0.08086980667844901
 
+
 #+begin_src f90 
 program energy_hydrogen
   implicit none
@@ -403,7 +401,7 @@ gfortran hydrogen.f90 energy_hydrogen.f90 -o energy_hydrogen
 :header-args:f90: :tangle variance_hydrogen.f90
 :END:
 
-   The variance of the local energy measures the intensity of the
+The variance of the local energy measures the magnitude of the
 fluctuations of the local energy around the average. If the local
 energy is constant (i.e. $\Psi$ is an eigenfunction of $\hat{H}$)
 the variance is zero.
@@ -413,7 +411,7 @@ gfortran hydrogen.f90 energy_hydrogen.f90 -o energy_hydrogen
 E_L(\mathbf{r}) - E \right]^2 \, d\mathbf{r}}{\int \left[\Psi(\mathbf{r}) \right]^2 d\mathbf{r}}
 $$
 
-   Compute an numerical estimate of the variance of the local energy
+Compute a numerical estimate of the variance of the local energy
 in a grid of $50\times50\times50$ points in the range  $(-5,-5,-5) \le \mathbf{r} \le (5,5,5)$.
      
 #+BEGIN_SRC python :results output :exports both
@@ -545,13 +543,19 @@ gfortran hydrogen.f90 variance_hydrogen.f90 -o variance_hydrogen
 
 Instead of computing the average energy as a numerical integration
 on a grid, we will do a Monte Carlo sampling, which is an extremely
-  efficient method to compute integrals in large dimensions.
+efficient method to compute integrals when the number of dimensions is
+large.
 
 Moreover, a Monte Carlo sampling will alow us to remove the bias due
 to the discretization of space, and compute a statistical confidence
 interval.
   
+
 ** Computation of the statistical error
+:PROPERTIES:
+:header-args:python: :tangle qmc_stats.py
+:header-args:f90: :tangle qmc_stats.f90
+:END:
 
 To compute the statistical error, you need to perform $M$
 independent Monte Carlo calculations. You will obtain $M$ different
@@ -579,7 +583,8 @@ $$
 Write a function returning the average and statistical error of an
 input array.
 
-#+BEGIN_SRC python :results output
+#+BEGIN_SRC python
+from math import sqrt
 def ave_error(arr):
     M = len(arr)
     assert (M>1)
@@ -588,12 +593,33 @@ def ave_error(arr):
     return (average, sqrt(variance/M))
 #+END_SRC
 
-#+RESULTS:
+#+BEGIN_SRC f90
+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
+     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
+#+END_SRC
    
 ** Uniform sampling in the box
+:PROPERTIES:
+:header-args:python: :tangle qmc_uniform.py
+:header-args:f90: :tangle qmc_uniform.f90
+:END:
 
-   Write a function to perform a Monte Carlo calculation of the
-   average energy. At every Monte Carlo step,
+In this section we write a function to perform a Monte Carlo
+calculation of the average energy.
+At every Monte Carlo step:
 
 - Draw 3 uniform random numbers in the interval $(-5,-5,-5) \le
   (x,y,z) \le (5,5,5)$
@@ -604,12 +630,16 @@ def ave_error(arr):
   Once all the steps have been computed, return the average energy
   computed on the Monte Carlo calculation.
 
-   Then, write a loop to perform 30 Monte Carlo runs, and compute the
-   average energy and the associated statistical error.
+  In the main program, write a loop to perform 30 Monte Carlo runs,
+  and compute the average energy and the associated statistical error.
 
   Compute the energy of the wave function with $a=0.9$.
    
-   #+BEGIN_SRC python 
+
+  #+BEGIN_SRC python :results output
+from hydrogen  import *
+from qmc_stats import *
+
 def MonteCarlo(a, nmax):
     E = 0.
     N = 0.
@@ -620,9 +650,7 @@ def MonteCarlo(a, nmax):
         N += w
         E += w * e_loc(a,r)
     return E/N
-   #+END_SRC
 
-   #+BEGIN_SRC python 
 a = 0.9
 nmax = 100000
 X = [MonteCarlo(a,nmax) for i in range(30)]
@@ -631,7 +659,59 @@ print(f"E = {E} +/- {deltaE}")
   #+END_SRC
 
   #+RESULTS:
-   : E = -0.4952626284319677 +/- 0.0006877988969872546
+  : E = -0.4956255109300764 +/- 0.0007082875482711226
+   
+#+BEGIN_SRC f90
+subroutine uniform_montecarlo(a,nmax,energy)
+  implicit none
+  double precision, intent(in)  :: a
+  integer         , intent(in)  :: nmax 
+  double precision, intent(out) :: energy
+  
+  integer*8 :: istep
+  
+  double precision :: norm, r(3), w
+
+  double precision, external :: e_loc, psi
+
+  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)
+  end do
+  energy = energy / norm
+end subroutine uniform_montecarlo
+
+program qmc
+  implicit none
+  double precision, parameter :: a = 0.9
+  integer         , parameter :: nmax = 100000
+  integer         , parameter :: nruns = 30
+
+  integer :: irun
+  double precision :: X(nruns)
+  double precision :: ave, err
+
+  do irun=1,nruns
+     call uniform_montecarlo(a,nmax,X(irun))
+  enddo
+  call ave_error(X,nruns,ave,err)
+  print *, 'E = ', ave, '+/-', err
+end program qmc
+#+END_SRC
+
+#+begin_src sh :results output :exports both
+gfortran hydrogen.f90 qmc_stats.f90 qmc_uniform.f90 -o qmc_uniform
+./qmc_uniform
+#+end_src
+
+#+RESULTS:
+:  E =  -0.49588321986667677      +/-   7.1758863546737969E-004
 
 ** Gaussian sampling 
 
@@ -655,7 +735,7 @@ print(f"E = {E} +/- {deltaE}")
 
 $$
 E \approx \frac{\sum_i w_i E_L(\mathbf{r}_i)}{\sum_i w_i}, \;\;
-   w_i = \frac{\left[\Psi(\mathbf{r}_i)\right]^2}{P(\mathbf{r})} \delta x\, \delta y\, \delta z
+w_i = \frac{\left[\Psi(\mathbf{r}_i)\right]^2}{P(\mathbf{r}_i)} \delta \mathbf{r}
 $$
 
 #+BEGIN_SRC python

qmc_stats.f90

@@ -0,0 +1,15 @@
+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
+     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

qmc_stats.py

@@ -0,0 +1,7 @@
+from math import sqrt
+def ave_error(arr):
+    M = len(arr)
+    assert (M>1)
+    average = sum(arr)/M
+    variance = 1./(M-1) * sum( [ (x - average)**2 for x in arr ] )
+    return (average, sqrt(variance/M))

qmc_uniform.f90

@@ -0,0 +1,41 @@
+subroutine uniform_montecarlo(a,nmax,energy)
+  implicit none
+  double precision, intent(in)  :: a
+  integer         , intent(in)  :: nmax 
+  double precision, intent(out) :: energy
+  
+  integer*8 :: istep
+  
+  double precision :: norm, r(3), w
+
+  double precision, external :: e_loc, psi
+
+  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)
+  end do
+  energy = energy / norm
+end subroutine uniform_montecarlo
+
+program qmc
+  implicit none
+  double precision, parameter :: a = 0.9
+  integer         , parameter :: nmax = 100000
+  integer         , parameter :: nruns = 30
+
+  integer :: irun
+  double precision :: X(nruns)
+  double precision :: ave, err
+
+  do irun=1,nruns
+     call uniform_montecarlo(a,nmax,X(irun))
+  enddo
+  call ave_error(X,nruns,ave,err)
+  print *, 'E = ', ave, '+/-', err
+end program qmc

qmc_uniform.py

@@ -0,0 +1,19 @@
+from hydrogen  import *
+from qmc_stats import *
+
+def MonteCarlo(a, nmax):
+    E = 0.
+    N = 0.
+    for istep in range(nmax):
+        r = np.random.uniform(-5., 5., (3))
+        w = psi(a,r)
+        w = w*w
+        N += w
+        E += w * e_loc(a,r)
+    return E/N
+
+a = 0.9
+nmax = 100000
+X = [MonteCarlo(a,nmax) for i in range(30)]
+E, deltaE = ave_error(X)
+print(f"E = {E} +/- {deltaE}")