This website contains the QMC tutorial of the 2021 LTTC winter school @@ -510,8 +514,8 @@ coordinates, etc).
For a given system with Hamiltonian \(\hat{H}\) and wave function \(\Psi\), we define the local energy as @@ -594,8 +598,8 @@ energy computed over these configurations:
In this section, we consider the hydrogen atom with the following @@ -624,8 +628,8 @@ To do that, we will compute the local energy and check whether it is constant.
You will now program all quantities needed to compute the local energy of the H atom for the given wave function. @@ -652,8 +656,8 @@ to catch the error.
@@ -698,8 +702,8 @@ and returns the potential.
Python @@ -740,8 +744,8 @@ and returns the potential.
@@ -776,8 +780,8 @@ input arguments, and returns a scalar.
Python @@ -804,8 +808,8 @@ input arguments, and returns a scalar.
@@ -886,8 +890,8 @@ Therefore, the local kinetic energy is
Python @@ -928,8 +932,8 @@ Therefore, the local kinetic energy is
@@ -988,8 +992,8 @@ are calling is yours.
Python @@ -1020,8 +1024,8 @@ are calling is yours.
@@ -1031,8 +1035,8 @@ Find the theoretical value of \(a\) for which \(\Psi\) is an eigenfunction of \(
The program you will write in this section will be written in @@ -1084,8 +1088,8 @@ In Fortran, you will need to compile all the source files together:
@@ -1179,8 +1183,8 @@ plot './data' index 0 using 1:2 with lines title 'a=0.1', \
Python @@ -1257,8 +1261,8 @@ plt.savefig("plot_py.png")
If the space is discretized in small volume elements \(\mathbf{r}_i\) @@ -1288,8 +1292,8 @@ The energy is biased because:
@@ -1360,8 +1364,8 @@ To compile the Fortran and run it:
Python @@ -1478,8 +1482,8 @@ a = 2.0000000000000000 E = -8.0869806678448772E-002
The variance of the local energy is a functional of \(\Psi\) @@ -1506,8 +1510,8 @@ energy can be used as a measure of the quality of a wave function.
@@ -1518,8 +1522,8 @@ Prove that :
\(\bar{E} = \langle E \rangle\) is a constant, so \(\langle \bar{E} @@ -1538,8 +1542,8 @@ Prove that :
@@ -1615,8 +1619,8 @@ To compile and run:
Python @@ -1755,8 +1759,8 @@ a = 2.0000000000000000 E = -8.0869806678448772E-002 s2 = 1.8068814
Numerical integration with deterministic methods is very efficient @@ -1772,8 +1776,8 @@ interval.
To compute the statistical error, you need to perform \(M\) @@ -1813,8 +1817,8 @@ And the confidence interval is given by
@@ -1854,8 +1858,8 @@ input array.
Python @@ -1916,8 +1920,8 @@ input array.
We will now perform our first Monte Carlo calculation to compute the @@ -1978,8 +1982,8 @@ compute the statistical error.
@@ -2081,8 +2085,8 @@ well as the index of the current step.
Python @@ -2188,8 +2192,8 @@ E = -0.48084122147238995 +/- 2.4983775878329355E-003
We will now use the square of the wave function to sample random @@ -2308,8 +2312,8 @@ All samples should be kept, from both accepted and rejected moves.
If the box is infinitely small, the ratio will be very close @@ -2344,8 +2348,8 @@ the same variable later on to store a time step.
@@ -2454,8 +2458,8 @@ Can you observe a reduction in the statistical error?
Python @@ -2602,8 +2606,8 @@ A = 0.50762633333333318 +/- 3.4601756760043725E-004
One can use more efficient numerical schemes to move the electrons by choosing a smarter expression for the transition probability. @@ -2724,8 +2728,8 @@ The algorithm of the previous exercise is only slighlty modified as:
To obtain Gaussian-distributed random numbers, you can apply the
@@ -2789,8 +2793,8 @@ In Python, you can use the
-
@@ -2832,8 +2836,8 @@ Write a function to compute the drift vector \(\frac{\nabla \Psi(\mathbf{r})}{\P
Python
@@ -2866,8 +2870,8 @@ Write a function to compute the drift vector \(\frac{\nabla \Psi(\mathbf{r})}{\P
@@ -2963,8 +2967,8 @@ Modify the previous program to introduce the drift-diffusion scheme.
Python
@@ -3152,8 +3156,8 @@ A = 0.62037333333333333 +/- 4.8970160591451110E-004
As we have seen, Variational Monte Carlo is a powerful method to
@@ -3170,8 +3174,8 @@ finding a near-exact numerical solution to the Schrödinger equation.
Consider the time-dependent Schrödinger equation:
@@ -3239,8 +3243,8 @@ system.
The diffusion equation of particles is given by
@@ -3320,8 +3324,8 @@ Therefore, in both cases, you are dealing with a "Bosonic" ground state.
In a molecular system, the potential is far from being constant
@@ -3419,8 +3423,8 @@ energies computed with the trial wave function.
\[
@@ -3481,8 +3485,8 @@ Defining \(\Pi(\mathbf{r},t) = \psi(\mathbf{r},\tau)
Instead of having a variable number of particles to simulate the
@@ -3571,13 +3575,13 @@ the DMC algorithm. However, its use reduces significantly the time-step error.
-
@@ -3675,6 +3679,222 @@ time \(\tau_{\text{max}}\) =100 a.u.
+Python
+
+Fortran
+
Change your PDMC code for one of the following:
@@ -3703,8 +3923,8 @@ And compute the ground state energy.
-3.4.2 Exercise 1
+3.4.2 Exercise 1
3.4.2.1 Solution solution2
+3.4.2.1 Solution solution
3.4.3 Exercise 2
+3.4.3 Exercise 2
3.4.3.1 Solution solution2
+3.4.3.1 Solution solution
4 Diffusion Monte Carlo
+4 Diffusion Monte Carlo
4.1 Schrödinger equation in imaginary time
+4.1 Schrödinger equation in imaginary time
4.2 Relation to diffusion
+4.2 Relation to diffusion
4.3 Importance sampling
+4.3 Importance sampling
4.3.1 Appendix : Details of the Derivation
+4.3.1 Appendix : Details of the Derivation
4.4 Pure Diffusion Monte Carlo
+4.4 Pure Diffusion Monte Carlo
4.5 Hydrogen atom
+4.5 Hydrogen atom
4.5.1 Exercise
+4.5.1 Exercise
gfortran hydrogen.f90 qmc_stats.f90 pdmc.f90 -o pdmc
./pdmc
+
+4.5.1.1 Solution solution
+#!/usr/bin/env python3
+
+from hydrogen import *
+from qmc_stats import *
+
+def MonteCarlo(a, nmax, dt, tau, Eref):
+ sq_dt = np.sqrt(dt)
+
+ energy = 0.
+ N_accep = 0
+ normalization = 0.
+
+ w = 1.
+ tau_current = 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):
+ el = e_loc(a,r_old)
+ w *= np.exp(-dt*(el - Eref))
+
+ normalization += w
+ energy += w * el
+
+ tau_current += dt
+
+ # Reset when tau is reached
+ if tau_current > tau:
+ w = 1.
+ tau_current = 0.
+
+ chi = np.random.normal(loc=0., scale=1.0, size=(3))
+
+ 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)*dt + prod
+
+ q = psi_new / psi_old
+ q = np.exp(-argexpo) * q*q
+
+ if np.random.uniform() <= q:
+ N_accep += 1
+ r_old = r_new
+ d_old = d_new
+ d2_old = d2_new
+ psi_old = psi_new
+
+ return energy/normalization, N_accep/nmax
+
+
+# Run simulation
+a = 1.2
+nmax = 100000
+dt = 0.05
+tau = 100.
+E_ref = -0.5
+
+X0 = [ MonteCarlo(a, nmax, dt, tau, E_ref) for i in range(30)]
+
+# Energy
+X = [ x for (x, _) in X0 ]
+E, deltaE = ave_error(X)
+print(f"E = {E} +/- {deltaE}")
+
+# Acceptance rate
+X = [ x for (_, x) in X0 ]
+A, deltaA = ave_error(X)
+print(f"A = {A} +/- {deltaA}")
+
+subroutine pdmc(a, dt, nmax, energy, accep, tau, E_ref)
+ implicit none
+ double precision, intent(in) :: a, dt, tau
+ integer*8 , intent(in) :: nmax
+ double precision, intent(out) :: energy, accep
+ double precision, intent(in) :: E_ref
+
+ integer*8 :: istep
+ 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 :: e, w, norm, tau_current
+
+ double precision, external :: e_loc, psi
+
+ sq_dt = dsqrt(dt)
+
+ ! Initialization
+ energy = 0.d0
+ n_accep = 0_8
+ norm = 0.d0
+
+ w = 1.d0
+ tau_current = 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
+ e = e_loc(a,r_old)
+ w = w * dexp(-dt*(e - E_ref))
+
+ norm = norm + w
+ energy = energy + w*e
+
+ tau_current = tau_current + dt
+
+ ! Reset when tau is reached
+ if (tau_current > tau) then
+ w = 1.d0
+ tau_current = 0.d0
+ endif
+
+ call random_gauss(chi,3)
+ 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)
+
+ 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)*dt + prod
+
+ q = psi_new / psi_old
+ q = dexp(-argexpo) * q*q
+
+ call random_number(u)
+
+ 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
+
+ end do
+
+ energy = energy / norm
+ accep = dble(n_accep) / dble(nmax)
+
+end subroutine pdmc
+
+program qmc
+ implicit none
+ double precision, parameter :: a = 1.2d0
+ double precision, parameter :: dt = 0.05d0
+ double precision, parameter :: E_ref = -0.5d0
+ double precision, parameter :: tau = 100.d0
+ integer*8 , parameter :: nmax = 100000
+ integer , parameter :: nruns = 30
+
+ integer :: irun
+ double precision :: X(nruns), accep(nruns)
+ double precision :: ave, err
+
+ do irun=1,nruns
+ call pdmc(a, dt, nmax, X(irun), accep(irun), tau, E_ref)
+ 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
+
+
+E = -0.49963953547336709 +/- 6.8755513992017491E-004
+A = 0.98963533333333342 +/- 6.3052128284666221E-005
+
5 Project
+5 Project
6 Acknowledgments
+6 Acknowledgments