VMC OK

QMC.org

@@ -1248,27 +1248,36 @@ gfortran hydrogen.f90 qmc_stats.f90 vmc.f90 -o vmc
       the move: set $\mathbf{r}_{n+1} = \mathbf{r}_{n}$, but *don't
       remove the sample from the average!*
 
+    The /acceptance rate/ is the ratio of the number of accepted step
+    over the total number of steps. The time step should be adapted so
+    that the acceptance rate is around 0.5 for a good efficiency of
+    the simulation.
 
 **** TODO Exercise 
 
      #+begin_exercise
      Modify the previous program to introduce the accept/reject step.
      You should recover the unbiased result.
+     Adjust the time-step so that the acceptance rate is 0.5.
      #+end_exercise
    
  *Python*
       #+BEGIN_SRC python :results output
+from hydrogen  import *
+from qmc_stats import *
+
 def MonteCarlo(a,tau,nmax):
     E = 0.
     N = 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)
     for istep in range(nmax):
-        eta = np.random.normal(loc=0., scale=1.0, size=(3))
+        chi = np.random.normal(loc=0., scale=1.0, size=(3))
-        r_new = r_old + tau * d_old + sq_tau * eta
+        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)
@@ -1278,75 +1287,111 @@ def MonteCarlo(a,tau,nmax):
         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
             psi_old = psi_new
         N += 1.
         E += e_loc(a,r_old)
-    return E/N
+    return E/N, accep_rate/N
 
 
+a = 0.9
 nmax = 100000
-tau = 0.1
+tau = 1.0
 X = [MonteCarlo(a,tau,nmax) for i in range(30)]
-E, deltaE = ave_error(X)
+E, deltaE = ave_error([x[0] for x in X])
-print(f"E = {E} +/- {deltaE}")
+A, deltaA = ave_error([x[1] for x in X])
+print(f"E = {E} +/- {deltaE}  {A} +/- {deltaA}")
       #+END_SRC
 
       #+RESULTS:
-      : E = -0.4951783346213532 +/- 0.00022067316984271938
+      : E = -0.49387078389332206 +/- 0.0033326460286729792  0.4983000000000001 +/- 0.006825097363627021
    
  *Fortran*
       #+BEGIN_SRC f90
-subroutine variational_montecarlo(a,nmax,energy)
+subroutine variational_montecarlo(a,tau,nmax,energy,accep_rate)
   implicit none
-  double precision, intent(in)  :: a
+  double precision, intent(in)  :: a, tau
-  integer         , intent(in)  :: nmax 
+  integer*8       , intent(in)  :: nmax 
-  double precision, intent(out) :: energy
+  double precision, intent(out) :: energy, accep_rate
 
   integer*8 :: istep
+  double precision :: norm, sq_tau, 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
 
-  double precision :: norm, r(3), w
+  sq_tau = dsqrt(tau)
-
-  double precision, external :: e_loc, psi, gaussian
   
+  ! 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)
+
   do istep = 1,nmax
-     call random_gauss(r,3)
+     call random_gauss(chi,3)
-     w = psi(a,r) 
+     r_new(:) = r_old(:) + tau * d_old(:) + chi(:)*sq_tau
-     w = w*w / gaussian(r)
+     call drift(a,r_new,d_new)
-     norm = norm + w
+     d2_new = d_new(1)*d_new(1) + d_new(2)*d_new(2) + d_new(3)*d_new(3)
-     energy = energy + w * e_loc(a,r)
+     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
 
 program qmc
   implicit none
   double precision, parameter :: a = 0.9
-  integer         , parameter :: nmax = 100000
+  double precision, parameter :: tau = 1.0
+  integer*8       , parameter :: nmax = 100000
   integer         , parameter :: nruns = 30
 
   integer :: irun
-  double precision :: X(nruns)
+  double precision :: X(nruns), accep(nruns)
   double precision :: ave, err
 
   do irun=1,nruns
-     call gaussian_montecarlo(a,nmax,X(irun))
+     call variational_montecarlo(a,tau,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
 
       #+begin_src sh :results output :exports both
-gfortran hydrogen.f90 qmc_stats.f90 vmc.f90 -o vmc
+gfortran hydrogen.f90 qmc_stats.f90 vmc_metropolis.f90 -o vmc_metropolis
-./vmc
+./vmc_metropolis
       #+end_src
 
       #+RESULTS:
+      :  E =  -0.49499990423528023      +/-   1.5958250761863871E-004
+      :  A =   0.78861366666666655      +/-   3.5096729498002445E-004
      
   
 * TODO Diffusion Monte Carlo