Changing a=0.9 to a=1.2 to avoid population explosion in DMC

QMC.org

@@ -1185,7 +1185,7 @@ end subroutine ave_error
 *** Exercise
 
     #+begin_exercise
-    Parameterize the wave function with $a=0.9$.  Perform 30
+    Parameterize the wave function with $a=1.2$.  Perform 30
     independent Monte Carlo runs, each with 100 000 Monte Carlo
     steps. Store the final energies of each run and use this array to
     compute the average energy and the associated error bar.
@@ -1206,7 +1206,7 @@ from qmc_stats import *
 def MonteCarlo(a, nmax):
      # TODO
 
-a    = 0.9
+a    = 1.2
 nmax = 100000
 
 #TODO
@@ -1244,7 +1244,7 @@ end subroutine uniform_montecarlo
 
 program qmc
   implicit none
-  double precision, parameter :: a = 0.9
+  double precision, parameter :: a = 1.2d0
   integer*8       , parameter :: nmax = 100000
   integer         , parameter :: nruns = 30
 
@@ -1287,7 +1287,7 @@ def MonteCarlo(a, nmax):
 
      return energy / normalization
 
-a    = 0.9
+a    = 1.2
 nmax = 100000
 
 X = [MonteCarlo(a,nmax) for i in range(30)]
@@ -1297,7 +1297,7 @@ print(f"E = {E} +/- {deltaE}")
      #+END_SRC
 
      #+RESULTS:
-     : E = -0.4956255109300764 +/- 0.0007082875482711226
+     : E = -0.4773221805255284 +/- 0.0022489426510479975
 
  *Fortran*
      #+begin_note
@@ -1341,7 +1341,7 @@ end subroutine uniform_montecarlo
 
 program qmc
   implicit none
-  double precision, parameter :: a     = 0.9
+  double precision, parameter :: a     = 1.2d0
   integer*8       , parameter :: nmax  = 100000
   integer         , parameter :: nruns = 30
 
@@ -1365,7 +1365,7 @@ gfortran hydrogen.f90 qmc_stats.f90 qmc_uniform.f90 -o qmc_uniform
      #+end_src
 
      #+RESULTS:
-     :  E =  -0.49518773675598715      +/-   5.2391494923686175E-004
+     :  E =  -0.48084122147238995      +/-   2.4983775878329355E-003
 
 ** Metropolis sampling with $\Psi^2$
    :PROPERTIES:
@@ -1420,10 +1420,14 @@ gfortran hydrogen.f90 qmc_stats.f90 qmc_uniform.f90 -o qmc_uniform
    which, for our choice of transition probability, becomes
    
    $$
-   A(\mathbf{r}_{n}\rightarrow\mathbf{r}_{n+1}) = \min\left(1,\frac{P(\mathbf{r}_{n+1})}{P(\mathbf{r}_{n})}\right)= \min\left(1,\frac{\Psi(\mathbf{r}_{n+1})^2}{\Psi(\mathbf{r}_{n})^2}\right)\,.
+   A(\mathbf{r}_{n}\rightarrow\mathbf{r}_{n+1}) = \min\left(1,\frac{P(\mathbf{r}_{n+1})}{P(\mathbf{r}_{n})}\right)= \min\left(1,\frac{|\Psi(\mathbf{r}_{n+1})|^2}{|\Psi(\mathbf{r}_{n})|^2}\right)\,.
    $$
    
-   Explain why the transition probability cancels out in the expression of $A$. Also note that we do not need to compute the norm of the wave function!
+   #+begin_exercise
+   Explain why the transition probability cancels out in the
+   expression of $A$.
+   #+end_exercise
+   Also note that we do not need to compute the norm of the wave function!
    
    The algorithm is summarized as follows:
    
@@ -1439,27 +1443,32 @@ gfortran hydrogen.f90 qmc_stats.f90 qmc_uniform.f90 -o qmc_uniform
     A common error is to remove the rejected samples from the
     calculation of the average. *Don't do it!*
 
-    All samples should be kept, from both accepted and rejected moves.
+    All samples should be kept, from both accepted /and/ rejected moves.
    #+end_note
    
-   If the box is infinitely small, the ratio will be very close
-   to one and all the steps will be accepted. However, the moves will be 
-   very correlated and you will visit the configurational space very slowly.
-
-   On the other hand, if you propose too large moves, the number of
-   accepted steps will decrease because the ratios might become
-   small. If the number of accepted steps is close to zero, then the
-   space is not well sampled either.
-
-   The size of the move should be adjusted so that it is as large as
-   possible, keeping the number of accepted steps not too small. To
-   achieve that, we define the acceptance rate as the number of
-   accepted steps over the total number of steps. Adjusting the time
-   step such that the acceptance rate is close to 0.5 is a good 
-   compromise for the current problem.
    
-   NOTE: below, we use the symbol $\delta t$ to denote $\delta L$ since we will use
-   the same variable later on to store a time step.
+*** Optimal step size
+    
+    If the box is infinitely small, the ratio will be very close
+    to one and all the steps will be accepted. However, the moves will be 
+    very correlated and you will visit the configurational space very slowly.
+
+    On the other hand, if you propose too large moves, the number of
+    accepted steps will decrease because the ratios might become
+    small. If the number of accepted steps is close to zero, then the
+    space is not well sampled either.
+
+    The size of the move should be adjusted so that it is as large as
+    possible, keeping the number of accepted steps not too small. To
+    achieve that, we define the acceptance rate as the number of
+    accepted steps over the total number of steps. Adjusting the time
+    step such that the acceptance rate is close to 0.5 is a good 
+    compromise for the current problem.
+   
+   #+begin_note
+    Below, we use the symbol $\delta t$ to denote $\delta L$ since we will use
+    the same variable later on to store a time step.
+   #+end_note
    
    
 *** Exercise
@@ -1489,7 +1498,7 @@ def MonteCarlo(a,nmax,dt):
 
 
 # Run simulation
-a    = 0.9
+a    = 1.2
 nmax = 100000
 dt   = #TODO
 
@@ -1506,6 +1515,8 @@ A, deltaA = ave_error(X)
 print(f"A = {A} +/- {deltaA}")
      #+END_SRC
 
+     #+RESULTS:
+
  *Fortran*
      #+BEGIN_SRC f90 :tangle none
 subroutine metropolis_montecarlo(a,nmax,dt,energy,accep)
@@ -1529,7 +1540,7 @@ end subroutine metropolis_montecarlo
 
 program qmc
   implicit none
-  double precision, parameter :: a     = 0.9d0
+  double precision, parameter :: a     = 1.2d0
   double precision, parameter :: dt    = ! TODO
   integer*8       , parameter :: nmax  = 100000
   integer         , parameter :: nruns = 30
@@ -1588,9 +1599,9 @@ def MonteCarlo(a,nmax,dt):
     return energy/nmax, N_accep/nmax
 
 # Run simulation
-a    = 0.9
+a    = 1.2
 nmax = 100000
-dt   = 1.3
+dt   = 1.0
 
 X0 = [ MonteCarlo(a,nmax,dt) for i in range(30)]
 
@@ -1606,8 +1617,8 @@ print(f"A = {A} +/- {deltaA}")
      #+END_SRC
 
      #+RESULTS:
-     : E = -0.4950720838131573 +/- 0.00019089638602238043
-     : A = 0.5172960000000001 +/- 0.0003443446549306529
+     : E = -0.4802595860693983 +/- 0.0005124420418289207
+     : A = 0.5074913333333334 +/- 0.000350889422714878
 
  *Fortran*
      #+BEGIN_SRC f90 :exports code
@@ -1662,8 +1673,8 @@ end subroutine metropolis_montecarlo
 
 program qmc
   implicit none
-  double precision, parameter :: a = 0.9d0
-  double precision, parameter :: dt = 1.3d0
+  double precision, parameter :: a = 1.2d0
+  double precision, parameter :: dt = 1.0d0
   integer*8       , parameter :: nmax = 100000
   integer         , parameter :: nruns = 30
 
@@ -1689,8 +1700,8 @@ gfortran hydrogen.f90 qmc_stats.f90 qmc_metropolis.f90 -o qmc_metropolis
 ./qmc_metropolis
      #+end_src
      #+RESULTS:
-     :  E =  -0.49503130891988767      +/-   1.7160104275040037E-004
-     :  A =   0.51695266666666673      +/-   4.0445505648997396E-004
+     :  E =  -0.47948142754168033      +/-   4.8410177863919307E-004
+     :  A =   0.50762633333333318      +/-   3.4601756760043725E-004
 
 ** Gaussian random number generator
    
@@ -1899,7 +1910,7 @@ def MonteCarlo(a,nmax,dt):
    # TODO
 
 # Run simulation
-a    = 0.9
+a    = 1.2
 nmax = 100000
 dt   = # TODO
 
@@ -1939,7 +1950,7 @@ end subroutine variational_montecarlo
 
 program qmc
   implicit none
-  double precision, parameter :: a     = 0.9
+  double precision, parameter :: a     = 1.2d0
   double precision, parameter :: dt    = ! TODO
   integer*8       , parameter :: nmax  = 100000
   integer         , parameter :: nruns = 30
@@ -2014,7 +2025,7 @@ def MonteCarlo(a,nmax,dt):
 
 
 # Run simulation
-a    = 0.9
+a    = 1.2
 nmax = 100000
 dt   = 1.3
 
@@ -2112,8 +2123,8 @@ end subroutine variational_montecarlo
 
 program qmc
   implicit none
-  double precision, parameter :: a     = 0.9
-  double precision, parameter :: dt    = 1.0
+  double precision, parameter :: a     = 1.2d0
+  double precision, parameter :: dt    = 1.0d0
   integer*8       , parameter :: nmax  = 100000
   integer         , parameter :: nruns = 30
 
@@ -2442,7 +2453,7 @@ def MonteCarlo(a, nmax, dt, Eref):
     # TODO
 
 # Run simulation
-a     = 0.9
+a     = 1.2
 nmax  = 100000
 dt    = 0.01
 E_ref = -0.5
@@ -2484,7 +2495,7 @@ end subroutine pdmc
 
 program qmc
   implicit none
-  double precision, parameter :: a     = 0.9
+  double precision, parameter :: a     = 1.2d0
   double precision, parameter :: dt    = 0.1d0
   double precision, parameter :: E_ref = -0.5d0
   double precision, parameter :: tau   = 10.d0
@@ -2576,7 +2587,7 @@ def MonteCarlo(a, nmax, dt, tau, Eref):
 
 
 # Run simulation
-a     = 0.9
+a     = 1.2
 nmax  = 100000
 dt    = 0.1
 tau   = 10.
@@ -2694,7 +2705,7 @@ end subroutine pdmc
 
 program qmc
   implicit none
-  double precision, parameter :: a     = 0.9
+  double precision, parameter :: a     = 1.2d0
   double precision, parameter :: dt    = 0.1d0
   double precision, parameter :: E_ref = -0.5d0
   double precision, parameter :: tau   = 10.d0
@@ -2814,7 +2825,7 @@ def MonteCarlo(a,nmax):
         E += w * e_loc(a,r)
     return E/N
 
-a = 0.9
+a = 1.2
 nmax = 100000
 X = [MonteCarlo(a,nmax) for i in range(30)]
 E, deltaE = ave_error(X)
@@ -2860,7 +2871,7 @@ end subroutine gaussian_montecarlo
 
 program qmc
   implicit none
-  double precision, parameter :: a = 0.9
+  double precision, parameter :: a = 1.2d0
   integer*8       , parameter :: nmax = 100000
   integer         , parameter :: nruns = 30
 
@@ -3003,7 +3014,7 @@ def MonteCarlo(a,dt,nmax):
     return E/N
 
 
-a = 0.9
+a = 1.2
 nmax = 100000
 dt = 0.2
 X = [MonteCarlo(a,dt,nmax) for i in range(30)]
@@ -3046,8 +3057,8 @@ end subroutine variational_montecarlo
 
 program qmc
   implicit none
-  double precision, parameter :: a = 0.9
-  double precision, parameter :: dt = 0.2
+  double precision, parameter :: a = 1.2d0
+  double precision, parameter :: dt = 0.2d0
   integer*8       , parameter :: nmax = 100000
   integer         , parameter :: nruns = 30
 
@@ -3093,4 +3104,6 @@ gfortran hydrogen.f90 qmc_stats.f90 vmc.f90 -o vmc
  | <2021-02-04 Thu 12:15-12:30> |     2.4 |
  |------------------------------+---------|
  | <2021-02-04 Thu 14:00-14:10> |     3.1 |
+ | <2021-02-04 Thu 14:10-14:30> |     3.2 |
+ | <2021-02-04 Thu 14:30-15:30> |     3.3 |
  |------------------------------+---------|

energy_hydrogen.py

@@ -1,3 +1,5 @@
+#!/usr/bin/env python3
+
 import numpy as np
 from hydrogen import e_loc, psi
 

hydrogen.f90

@@ -43,7 +43,8 @@ double precision function e_loc(a,r)
   implicit none
   double precision, intent(in) :: a, r(3)
 
-  double precision, external   :: kinetic, potential
+  double precision, external :: kinetic
+  double precision, external :: potential
 
   e_loc = kinetic(a,r) + potential(r)
 

hydrogen.py

@@ -1,3 +1,4 @@
+#!/usr/bin/env python3
 import numpy as np
 
 def potential(r):

plot_hydrogen.py

@@ -1,3 +1,5 @@
+#!/usr/bin/env python3
+
 import numpy as np
 import matplotlib.pyplot as plt
 

qmc_gaussian.f90

@@ -32,7 +32,7 @@ end subroutine gaussian_montecarlo
 
 program qmc
   implicit none
-  double precision, parameter :: a = 0.9
+  double precision, parameter :: a = 1.2d0
   integer*8       , parameter :: nmax = 100000
   integer         , parameter :: nruns = 30
 

qmc_gaussian.py

@@ -1,3 +1,5 @@
+#!/usr/bin/env python3
+
 from hydrogen  import *
 from qmc_stats import *
 
@@ -16,7 +18,7 @@ def MonteCarlo(a,nmax):
         E += w * e_loc(a,r)
     return E/N
 
-a = 0.9
+a = 1.2
 nmax = 100000
 X = [MonteCarlo(a,nmax) for i in range(30)]
 E, deltaE = ave_error(X)

qmc_metropolis.f90

@@ -49,8 +49,8 @@ end subroutine metropolis_montecarlo
 
 program qmc
   implicit none
-  double precision, parameter :: a = 0.9d0
-  double precision, parameter :: dt = 1.3d0
+  double precision, parameter :: a = 1.2d0
+  double precision, parameter :: dt = 1.0d0
   integer*8       , parameter :: nmax = 100000
   integer         , parameter :: nruns = 30
 

qmc_metropolis.py

@@ -1,3 +1,5 @@
+#!/usr/bin/env python3
+
 from hydrogen  import *
 from qmc_stats import *
 
@@ -25,9 +27,9 @@ def MonteCarlo(a,nmax,dt):
     return energy/nmax, N_accep/nmax
 
 # Run simulation
-a    = 0.9
+a    = 1.2
 nmax = 100000
-dt   = 1.3
+dt   = 1.0
 
 X0 = [ MonteCarlo(a,nmax,dt) for i in range(30)]
 

qmc_stats.py

@@ -1,3 +1,5 @@
+#!/usr/bin/env python3
+
 from math import sqrt
 def ave_error(arr):
     M = len(arr)

qmc_uniform.f90

@@ -31,7 +31,7 @@ end subroutine uniform_montecarlo
 
 program qmc
   implicit none
-  double precision, parameter :: a     = 0.9
+  double precision, parameter :: a     = 1.2d0
   integer*8       , parameter :: nmax  = 100000
   integer         , parameter :: nruns = 30
 

qmc_uniform.py

@@ -1,3 +1,5 @@
+#!/usr/bin/env python3
+
 from hydrogen  import *
 from qmc_stats import *
 
@@ -16,7 +18,7 @@ def MonteCarlo(a, nmax):
 
      return energy / normalization
 
-a    = 0.9
+a    = 1.2
 nmax = 100000
 
 X = [MonteCarlo(a,nmax) for i in range(30)]

variance_hydrogen.py

@@ -1,3 +1,5 @@
+#!/usr/bin/env python3
+
 import numpy as np
 from hydrogen import e_loc, psi
 

vmc.f90

@@ -28,8 +28,8 @@ end subroutine variational_montecarlo
 
 program qmc
   implicit none
-  double precision, parameter :: a = 0.9
-  double precision, parameter :: dt = 0.2
+  double precision, parameter :: a = 1.2d0
+  double precision, parameter :: dt = 0.2d0
   integer*8       , parameter :: nmax = 100000
   integer         , parameter :: nruns = 30
 

vmc.py

@@ -1,3 +1,5 @@
+#!/usr/bin/env python3
+
 from hydrogen  import *
 from qmc_stats import *
 
@@ -19,7 +21,7 @@ def MonteCarlo(a,dt,nmax):
     return E/N
 
 
-a = 0.9
+a = 1.2
 nmax = 100000
 dt = 0.2
 X = [MonteCarlo(a,dt,nmax) for i in range(30)]

vmc_metropolis.f90

@@ -73,8 +73,8 @@ end subroutine variational_montecarlo
 
 program qmc
   implicit none
-  double precision, parameter :: a     = 0.9
-  double precision, parameter :: dt    = 1.0
+  double precision, parameter :: a     = 1.2d0
+  double precision, parameter :: dt    = 1.0d0
   integer*8       , parameter :: nmax  = 100000
   integer         , parameter :: nruns = 30
 

vmc_metropolis.py

@@ -1,3 +1,5 @@
+#!/usr/bin/env python3
+
 from hydrogen  import *
 from qmc_stats import *
 
@@ -41,7 +43,7 @@ def MonteCarlo(a,nmax,dt):
 
 
 # Run simulation
-a    = 0.9
+a    = 1.2
 nmax = 100000
 dt   = 1.3