Gaussian sampling

QMC.org

@@ -550,7 +550,6 @@ 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
@@ -635,7 +634,6 @@ At every Monte Carlo step:
 
      Compute the energy of the wave function with $a=0.9$.
 
-
    #+BEGIN_SRC python :results output
 from hydrogen  import *
 from qmc_stats import *
@@ -714,13 +712,57 @@ gfortran hydrogen.f90 qmc_stats.f90 qmc_uniform.f90 -o qmc_uniform
    :  E =  -0.49588321986667677      +/-   7.1758863546737969E-004
 
 ** Gaussian sampling 
+   :PROPERTIES:
+   :header-args:python: :tangle qmc_gaussian.py
+   :header-args:f90: :tangle qmc_gaussian.f90
+   :END:
 
    We will now improve the sampling and allow to sample in the whole
    3D space, correcting the bias related to the sampling in the box.
 
    Instead of drawing uniform random numbers, we will draw Gaussian
-random numbers centered on 0 and with a variance of 1. Now the
+   random numbers centered on 0 and with a variance of 1.
-equation for the energy is changed into
+
+   To obtain Gaussian-distributed random numbers, you can apply the
+   [[https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform][Box Muller transform]] to uniform random numbers:
+
+   \begin{eqnarray*}
+   z_1 &=& \sqrt{-2 \ln u_1} \cos(2 \pi u_2) \\
+   z_2 &=& \sqrt{-2 \ln u_1} \sin(2 \pi u_2) 
+   \end{eqnarray*}
+
+   #+BEGIN_SRC f90 :tangle qmc_stats.f90
+subroutine random_gauss(z,n)
+  implicit none
+  integer, intent(in) :: n
+  double precision, intent(out) :: z(n)
+  double precision :: u(n+1)
+  double precision, parameter :: two_pi = 2.d0*dacos(-1.d0)
+  integer :: i
+
+  call random_number(u)
+  if (iand(n,1) == 0) then
+     ! n is even
+     do i=1,n,2
+        z(i)   = dsqrt(-2.d0*dlog(u(i))) 
+        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
+        z(i)   = dsqrt(-2.d0*dlog(u(i))) 
+        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
+
+
+   Now the equation for the energy is changed into
    
    \[
    E = \frac{\int P(\mathbf{r}) \frac{\left[\Psi(\mathbf{r})\right]^2}{P(\mathbf{r})}\, \frac{\hat{H} \Psi(\mathbf{r})}{\Psi(\mathbf{r})}\,d\mathbf{r}}{\int P(\mathbf{r}) \frac{\left[\Psi(\mathbf{r}) \right]^2}{P(\mathbf{r})} d\mathbf{r}}
@@ -738,15 +780,14 @@ 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}_i)} \delta \mathbf{r}
    $$
 
-#+BEGIN_SRC python
+   #+BEGIN_SRC python :results output
+from hydrogen  import *
+from qmc_stats import *
+
 norm_gauss = 1./(2.*np.pi)**(1.5)
 def gaussian(r):
     return norm_gauss * np.exp(-np.dot(r,r)*0.5)
-#+END_SRC
 
-#+RESULTS:
-
-#+BEGIN_SRC python
 def MonteCarlo(a,nmax):
     E = 0.
     N = 0.
@@ -757,11 +798,7 @@ def MonteCarlo(a,nmax):
         N += w
         E += w * e_loc(a,r)
     return E/N
-#+END_SRC
 
-#+RESULTS:
-
-#+BEGIN_SRC python :results output
 a = 0.9
 nmax = 100000
 X = [MonteCarlo(a,nmax) for i in range(30)]
@@ -770,8 +807,67 @@ print(f"E = {E} +/- {deltaE}")
    #+END_SRC
 
    #+RESULTS:
-: E = -0.4952488228427792 +/- 0.00011913174676540714
+   : E = -0.49507506093129827 +/- 0.00014164037765553668
 
+   
+   #+BEGIN_SRC f90
+double precision function gaussian(r)
+  implicit none
+  double precision, intent(in) :: r(3)
+  double precision, parameter :: norm_gauss = 1.d0/(2.d0*dacos(-1.d0))**(1.5d0)
+  gaussian = norm_gauss * dexp( -0.5d0 * dsqrt(r(1)*r(1) + r(2)*r(2) + r(3)*r(3) ))
+end function gaussian
+
+
+subroutine gaussian_montecarlo(a,nmax,energy)
+  implicit none
+  double precision, intent(in)  :: a
+  integer         , intent(in)  :: nmax 
+  double precision, intent(out) :: energy
+
+  integer*8 :: istep
+
+  double precision :: norm, r(3), w
+
+  double precision, external :: e_loc, psi, gaussian
+
+  energy = 0.d0
+  norm   = 0.d0
+  do istep = 1,nmax
+     call random_gauss(r,3)
+     w = psi(a,r) 
+     w = w*w / gaussian(r)
+     norm = norm + w
+     energy = energy + w * e_loc(a,r)
+  end do
+  energy = energy / norm
+end subroutine gaussian_montecarlo
+
+program qmc
+  implicit none
+  double precision, parameter :: a = 0.9
+  integer         , parameter :: nmax = 100000
+  integer         , parameter :: nruns = 30
+
+  integer :: irun
+  double precision :: X(nruns)
+  double precision :: ave, err
+
+  do irun=1,nruns
+     call gaussian_montecarlo(a,nmax,X(irun))
+  enddo
+  call ave_error(X,nruns,ave,err)
+  print *, 'E = ', ave, '+/-', err
+end program qmc
+   #+END_SRC
+
+   #+begin_src sh :results output :exports both
+gfortran hydrogen.f90 qmc_stats.f90 qmc_gaussian.f90 -o qmc_gaussian
+./qmc_gaussian
+   #+end_src
+
+   #+RESULTS:
+   :  E =  -0.49606057056767766      +/-   1.3918807547836872E-004
 ** Sampling with $\Psi^2$
 
    We will now use the square of the wave function to make the sampling:

qmc_stats.f90

@@ -13,3 +13,31 @@ subroutine ave_error(x,n,ave,err)
      err = dsqrt(variance/dble(n))
   endif
 end subroutine ave_error
+
+subroutine random_gauss(z,n)
+  implicit none
+  integer, intent(in) :: n
+  double precision, intent(out) :: z(n)
+  double precision :: u(n+1)
+  double precision, parameter :: two_pi = 2.d0*dacos(-1.d0)
+  integer :: i
+
+  call random_number(u)
+  if (iand(n,1) == 0) then
+     ! n is even
+     do i=1,n,2
+        z(i)   = dsqrt(-2.d0*dlog(u(i))) 
+        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
+        z(i)   = dsqrt(-2.d0*dlog(u(i))) 
+        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