diff --git a/qmc_gaussian.f90 b/qmc_gaussian.f90 new file mode 100644 index 0000000..6adfc3f --- /dev/null +++ b/qmc_gaussian.f90 @@ -0,0 +1,48 @@ +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 * (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*8 , 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*8 , 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 diff --git a/qmc_gaussian.py b/qmc_gaussian.py new file mode 100644 index 0000000..129c79f --- /dev/null +++ b/qmc_gaussian.py @@ -0,0 +1,23 @@ +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) + +def MonteCarlo(a,nmax): + E = 0. + N = 0. + for istep in range(nmax): + r = np.random.normal(loc=0., scale=1.0, size=(3)) + w = psi(a,r) + w = w*w / gaussian(r) + 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}") diff --git a/vmc.f90 b/vmc.f90 new file mode 100644 index 0000000..42ee411 --- /dev/null +++ b/vmc.f90 @@ -0,0 +1,45 @@ +subroutine variational_montecarlo(a,tau,nmax,energy) + implicit none + double precision, intent(in) :: a, tau + 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, external :: e_loc + + sq_tau = dsqrt(tau) + + ! Initialization + energy = 0.d0 + norm = 0.d0 + call random_gauss(r_old,3) + + 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 + norm = norm + 1.d0 + energy = energy + e_loc(a,r_new) + r_old(:) = r_new(:) + end do + energy = energy / norm +end subroutine variational_montecarlo + +program qmc + implicit none + double precision, parameter :: a = 0.9 + double precision, parameter :: tau = 0.2 + integer*8 , parameter :: nmax = 100000 + integer , parameter :: nruns = 30 + + integer :: irun + double precision :: X(nruns) + double precision :: ave, err + + do irun=1,nruns + call variational_montecarlo(a,tau,nmax,X(irun)) + enddo + call ave_error(X,nruns,ave,err) + print *, 'E = ', ave, '+/-', err +end program qmc diff --git a/vmc.py b/vmc.py new file mode 100644 index 0000000..ca422be --- /dev/null +++ b/vmc.py @@ -0,0 +1,27 @@ +from hydrogen import * +from qmc_stats import * + +def MonteCarlo(a,tau,nmax): + sq_tau = np.sqrt(tau) + + # Initialization + E = 0. + N = 0. + r_old = np.random.normal(loc=0., scale=1.0, size=(3)) + + 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 + N += 1. + E += e_loc(a,r_new) + r_old = r_new + return E/N + + +a = 0.9 +nmax = 100000 +tau = 0.2 +X = [MonteCarlo(a,tau,nmax) for i in range(30)] +E, deltaE = ave_error(X) +print(f"E = {E} +/- {deltaE}") diff --git a/vmc_metropolis.f90 b/vmc_metropolis.f90 new file mode 100644 index 0000000..3d843e2 --- /dev/null +++ b/vmc_metropolis.f90 @@ -0,0 +1,71 @@ +subroutine variational_montecarlo(a,tau,nmax,energy,accep_rate) + implicit none + double precision, intent(in) :: a, tau + integer*8 , intent(in) :: nmax + 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 + + sq_tau = dsqrt(tau) + + ! 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(chi,3) + r_new(:) = r_old(:) + tau * d_old(:) + chi(:)*sq_tau + 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) + 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