PDMC codes OK

energy_hydrogen.f90

@@ -12,7 +12,9 @@ program energy_hydrogen
   end do
 
   do j=1,size(a)
+
      ! TODO
+
      print *, 'a = ', a(j), '    E = ', energy
   end do
 
@@ -38,19 +40,27 @@ program energy_hydrogen
   do j=1,size(a)
      energy = 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
+
               energy = energy + w * e_loc(a(j), r)
               norm   = norm   + w 
            end do
+
         end do
+
      end do
+
      energy = energy / norm
      print *, 'a = ', a(j), '    E = ', energy
   end do

energy_hydrogen.py

@@ -9,15 +9,19 @@ r = np.array([0.,0.,0.])
 for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]:
     E    = 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    += w * e_loc(a,r)
                 norm += w 
+
     E = E / norm
     print(f"a = {a} \t E = {E}")

hydrogen.f90

@@ -1,45 +1,62 @@
 double precision function potential(r)
   implicit none
   double precision, intent(in) :: r(3)
+
   double precision             :: distance
+
   distance = dsqrt( r(1)*r(1) + r(2)*r(2) + r(3)*r(3) )
+
   if (distance > 0.d0) then
-     potential = -1.d0 / dsqrt( r(1)*r(1) + r(2)*r(2) + r(3)*r(3) )
+     potential = -1.d0 / distance
   else
      stop 'potential at r=0.d0 diverges'
   end if
+
 end function potential
 
 double precision function psi(a, r)
   implicit none
   double precision, intent(in) :: a, r(3)
+
   psi = dexp(-a * dsqrt( r(1)*r(1) + r(2)*r(2) + r(3)*r(3) ))
 end function psi
 
 double precision function kinetic(a,r)
   implicit none
   double precision, intent(in) :: a, r(3)
+
   double precision             :: distance
+
   distance = dsqrt( r(1)*r(1) + r(2)*r(2) + r(3)*r(3) ) 
+
   if (distance > 0.d0) then
-     kinetic = -0.5d0 * (a*a - (2.d0*a) / distance)
+
+     kinetic =  a * (1.d0 / distance - 0.5d0 * a)
+
   else
      stop 'kinetic energy diverges at r=0'
   end if
+
 end function kinetic
 
 double precision function e_loc(a,r)
   implicit none
   double precision, intent(in) :: a, r(3)
+
   double precision, external   :: kinetic, potential
+
   e_loc = kinetic(a,r) + potential(r)
+
 end function e_loc
 
 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

hydrogen.py

@@ -11,7 +11,8 @@ def psi(a, r):
 def kinetic(a,r):
     distance = np.sqrt(np.dot(r,r))
     assert (distance > 0.)
-    return -0.5 * (a**2 - (2.*a)/distance)
+
+    return a * (1./distance - 0.5 * a)
 
 def e_loc(a,r):
     return kinetic(a,r) + potential(r)

qmc_metropolis.f90

@@ -1,44 +1,56 @@
-subroutine metropolis_montecarlo(a,nmax,tau,energy,accep)
+subroutine metropolis_montecarlo(a,nmax,dt,energy,accep)
   implicit none
   double precision, intent(in)  :: a
   integer*8       , intent(in)  :: nmax 
-  double precision, intent(in)  :: tau
+  double precision, intent(in)  :: dt
   double precision, intent(out) :: energy
   double precision, intent(out) :: accep
 
-  integer*8 :: istep
-
   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(:) = tau * (2.d0*r_old(:) - 1.d0)
+  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(:) + tau * (2.d0*r_new(:) - 1.d0)
+     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
-        n_accep = n_accep + 1_8
+
      endif
-     energy = energy + e_loc(a,r_old)
+
   end do
+
   energy = energy / dble(nmax)
   accep  = dble(n_accep) / dble(nmax)
+
 end subroutine metropolis_montecarlo
 
 program qmc
   implicit none
   double precision, parameter :: a = 0.9d0
-  double precision, parameter :: tau = 1.3d0
+  double precision, parameter :: dt = 1.3d0
   integer*8       , parameter :: nmax = 100000
   integer         , parameter :: nruns = 30
 
@@ -47,7 +59,7 @@ program qmc
   double precision :: ave, err
 
   do irun=1,nruns
-     call metropolis_montecarlo(a,nmax,tau,X(irun),Y(irun))
+     call metropolis_montecarlo(a,nmax,dt,X(irun),Y(irun))
   enddo
 
   call ave_error(X,nruns,ave,err)
@@ -55,4 +67,5 @@ program qmc
 
   call ave_error(Y,nruns,ave,err)
   print *, 'A = ', ave, '+/-', err
+
 end program qmc

qmc_metropolis.py

@@ -1,28 +1,35 @@
 from hydrogen  import *
 from qmc_stats import *
 
-def MonteCarlo(a,nmax,tau):
+def MonteCarlo(a,nmax,dt):
     energy  = 0.
     N_accep = 0
-    r_old = np.random.uniform(-tau, tau, (3))
+
+    r_old = np.random.uniform(-dt, dt, (3))
     psi_old = psi(a,r_old)
+
     for istep in range(nmax):
-        r_new = r_old + np.random.uniform(-tau,tau,(3))
+        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
-        v = np.random.uniform()
+
-        if v <= ratio:
+        if np.random.uniform() <= ratio:
             N_accep += 1
+
             r_old   = r_new
             psi_old = psi_new
-        energy += e_loc(a,r_old)
+
     return energy/nmax, N_accep/nmax
 
 # Run simulation
 a    = 0.9
 nmax = 100000
-tau = 1.3
+dt   = 1.3
-X0 = [ MonteCarlo(a,nmax,tau) for i in range(30)]
+
+X0 = [ MonteCarlo(a,nmax,dt) for i in range(30)]
 
 # Energy
 X = [ x for (x, _) in X0 ]

qmc_stats.f90

@@ -1,26 +1,25 @@
 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
-  ! TODO
-end subroutine ave_error
 
-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
 
@@ -33,6 +32,7 @@ subroutine random_gauss(z,n)
   integer :: i
 
   call random_number(u)
+
   if (iand(n,1) == 0) then
      ! n is even
      do i=1,n,2
@@ -40,6 +40,7 @@ subroutine random_gauss(z,n)
         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
@@ -47,7 +48,10 @@ subroutine random_gauss(z,n)
         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

qmc_stats.py

@@ -2,10 +2,14 @@ from math import sqrt
 def ave_error(arr):
     M = len(arr)
     assert(M>0)
+
     if M == 1:
-        return (arr[0], 0.)
+        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)

qmc_uniform.f90

@@ -5,22 +5,28 @@ subroutine uniform_montecarlo(a,nmax,energy)
   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)
+     norm   = norm   + w
+
   end do
+
   energy = energy / norm
+
 end subroutine uniform_montecarlo
 
 program qmc
@@ -36,6 +42,8 @@ program qmc
   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

@@ -4,16 +4,22 @@ 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
-          normalization += w
+
           energy        += w * e_loc(a,r)
+          normalization += w
+
      return energy / normalization
 
 a    = 0.9
 nmax = 100000
+
 X = [MonteCarlo(a,nmax) for i in range(30)]
 E, deltaE = ave_error(X)
+
 print(f"E = {E} +/- {deltaE}")

variance_hydrogen.f90

@@ -1,10 +1,12 @@
 program variance_hydrogen
   implicit none
-  double precision, external :: e_loc, psi
+
-  double precision :: x(50), w, delta, energy, dx, r(3), a(6), norm, s2
+  double precision :: x(50), w, delta, energy, energy2
-  double precision :: e, 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)
@@ -20,24 +22,35 @@ program variance_hydrogen
      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 = e_loc(a(j), r)
+
-              energy  = energy  + w * e
+              e_tmp = e_loc(a(j), r)
-              energy2 = energy2 + w * e * e
+
+              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
 

variance_hydrogen.py

@@ -2,6 +2,7 @@ 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.])
@@ -10,19 +11,26 @@ 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
-                El = e_loc(a, r)
+
-                E  += w * El
+                e_tmp = e_loc(a,r)
-                E2 += w * El*El
+                E    += w * e_tmp
+                E2   += w * e_tmp * e_tmp
                 norm += w 
+
     E  = E  / norm
     E2 = E2 / norm
-    s2 = E2 - E*E
+
+    s2 = E2 - E**2
     print(f"a = {a} \t E = {E:10.8f} \t \sigma^2 = {s2:10.8f}")

vmc.f90

@@ -1,14 +1,14 @@
-subroutine variational_montecarlo(a,tau,nmax,energy)
+subroutine variational_montecarlo(a,dt,nmax,energy)
   implicit none
-  double precision, intent(in)  :: a, tau
+  double precision, intent(in)  :: a, dt
   integer*8       , intent(in)  :: nmax 
   double precision, intent(out) :: energy
 
   integer*8 :: istep
-  double precision :: norm, r_old(3), r_new(3), d_old(3), sq_tau, chi(3)
+  double precision :: norm, r_old(3), r_new(3), d_old(3), sq_dt, chi(3)
   double precision, external :: e_loc
 
-  sq_tau = dsqrt(tau)
+  sq_dt = dsqrt(dt)
   
   ! Initialization
   energy = 0.d0
@@ -18,7 +18,7 @@ subroutine variational_montecarlo(a,tau,nmax,energy)
   do istep = 1,nmax
      call drift(a,r_old,d_old)
      call random_gauss(chi,3)
-     r_new(:) = r_old(:) + tau * d_old(:) + chi(:)*sq_tau
+     r_new(:) = r_old(:) + dt * d_old(:) + chi(:)*sq_dt
      norm = norm + 1.d0
      energy = energy + e_loc(a,r_new)
      r_old(:) = r_new(:)
@@ -29,7 +29,7 @@ end subroutine variational_montecarlo
 program qmc
   implicit none
   double precision, parameter :: a = 0.9
-  double precision, parameter :: tau = 0.2
+  double precision, parameter :: dt = 0.2
   integer*8       , parameter :: nmax = 100000
   integer         , parameter :: nruns = 30
 
@@ -38,7 +38,7 @@ program qmc
   double precision :: ave, err
 
   do irun=1,nruns
-     call variational_montecarlo(a,tau,nmax,X(irun))
+     call variational_montecarlo(a,dt,nmax,X(irun))
   enddo
   call ave_error(X,nruns,ave,err)
   print *, 'E = ', ave, '+/-', err

vmc.py

@@ -1,8 +1,8 @@
 from hydrogen  import *
 from qmc_stats import *
 
-def MonteCarlo(a,tau,nmax):
+def MonteCarlo(a,dt,nmax):
-    sq_tau = np.sqrt(tau)
+    sq_dt = np.sqrt(dt)
 
     # Initialization
     E = 0.
@@ -12,7 +12,7 @@ def MonteCarlo(a,tau,nmax):
     for istep in range(nmax):
         d_old = drift(a,r_old)
         chi = np.random.normal(loc=0., scale=1.0, size=(3))
-        r_new = r_old + tau * d_old + chi*sq_tau
+        r_new = r_old + dt * d_old + chi*sq_dt
         N += 1.
         E += e_loc(a,r_new)
         r_old = r_new
@@ -21,7 +21,7 @@ def MonteCarlo(a,tau,nmax):
 
 a = 0.9
 nmax = 100000
-tau = 0.2
+dt = 0.2
-X = [MonteCarlo(a,tau,nmax) for i in range(30)]
+X = [MonteCarlo(a,dt,nmax) for i in range(30)]
 E, deltaE = ave_error(X)
 print(f"E = {E} +/- {deltaE}")

vmc_metropolis.f90

@@ -1,57 +1,80 @@
-subroutine variational_montecarlo(a,tau,nmax,energy,accep_rate)
+subroutine variational_montecarlo(a,dt,nmax,energy,accep)
   implicit none
-  double precision, intent(in)  :: a, tau
+  double precision, intent(in)  :: a, dt
   integer*8       , intent(in)  :: nmax 
-  double precision, intent(out) :: energy, accep_rate
+  double precision, intent(out) :: energy, accep
 
   integer*8        :: istep
-  double precision :: sq_tau, chi(3), d2_old, prod, u
+  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_tau = dsqrt(tau)
+  sq_dt = dsqrt(dt)
 
   ! Initialization
   energy  = 0.d0
-  accep_rate = 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)
+  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(:) + tau * d_old(:) + chi(:)*sq_tau
+     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)
+     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)*tau + prod
+
+     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
+
-        accep_rate = accep_rate + 1.d0
+     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
-     energy = energy + e_loc(a,r_old)
+
   end do
+
   energy = energy / dble(nmax)
-  accep_rate = dble(accep_rate) / dble(nmax)
+  accep  = dble(n_accep) / dble(nmax)
+
 end subroutine variational_montecarlo
 
 program qmc
   implicit none
   double precision, parameter :: a     = 0.9
-  double precision, parameter :: tau = 1.0
+  double precision, parameter :: dt    = 1.0
   integer*8       , parameter :: nmax  = 100000
   integer         , parameter :: nruns = 30
 
@@ -60,10 +83,13 @@ program qmc
   double precision :: ave, err
 
   do irun=1,nruns
-     call variational_montecarlo(a,tau,nmax,X(irun),accep(irun))
+     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

vmc_metropolis.py

@@ -1,40 +1,51 @@
 from hydrogen  import *
 from qmc_stats import *
 
-def MonteCarlo(a,nmax,tau):
+def MonteCarlo(a,nmax,dt):
-    E = 0.
+    sq_dt = np.sqrt(dt)
-    accep_rate = 0.
+
-    sq_tau = np.sqrt(tau)
+    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))
-        r_new = r_old + tau * d_old + sq_tau * chi
+
+        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)*tau + prod
+        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:
+
-            accep_rate += 1.
+        if np.random.uniform() <= q:
+            N_accep += 1
+
             r_old   = r_new
             d_old   = d_new
             d2_old  = d2_new
             psi_old = psi_new
-        E += e_loc(a,r_old)
+
-    return E/nmax, accep_rate/nmax
+    return energy/nmax, accep_rate/nmax
 
 
 # Run simulation
 a    = 0.9
 nmax = 100000
-tau = 1.3
+dt   = 1.3
-X0 = [ MonteCarlo(a,nmax,tau) for i in range(30)]
+
+X0 = [ MonteCarlo(a,nmax,dt) for i in range(30)]
 
 # Energy
 X = [ x for (x, _) in X0 ]