From e6ca7dd0a513d83edca750ecd84e52202a3d932b Mon Sep 17 00:00:00 2001 From: scemama Date: Thu, 21 Jan 2021 22:25:25 +0000 Subject: [PATCH] deploy: a2373b198b49eb68afa426682bab03083782a5c6 --- index.html | 553 ++++++++++++++++++++++++++++++++++++++--------------- 1 file changed, 398 insertions(+), 155 deletions(-) diff --git a/index.html b/index.html index 4caf686..39d9b3e 100644 --- a/index.html +++ b/index.html @@ -3,7 +3,7 @@ "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd"> - + Quantum Monte Carlo @@ -257,74 +257,113 @@ for the JavaScript code in this tag.

Table of Contents

- -
-

1 Introduction

+
+

1 Introduction

We propose different exercises to understand quantum Monte Carlo (QMC) @@ -357,7 +396,7 @@ is defined everywhere, continuous and infinitely differentiable.

In Fortran, when you use a double precision constant, don't forget -to put d0 as a suffix (for example 2.0d0), or it will be +to put d0 as a suffix (for example 2.0d0), or it will be interpreted as a single precision value

@@ -366,8 +405,8 @@ interpreted as a single precision value
-
-

2 Numerical evaluation of the energy

+
+

2 Numerical evaluation of the energy

In this section we consider the Hydrogen atom with the following @@ -381,7 +420,8 @@ wave function:

-We will first verify that \(\Psi\) is an eigenfunction of the Hamiltonian +We will first verify that, for a given value of \(a\), \(\Psi\) is an +eigenfunction of the Hamiltonian

@@ -391,8 +431,7 @@ We will first verify that \(\Psi\) is an eigenfunction of the Hamiltonian

-when \(a=1\), by checking that \(\hat{H}\Psi(\mathbf{r}) = E\Psi(\mathbf{r})\) for -all \(\mathbf{r}\). We will check that the local energy, defined as +To do that, we will check if the local energy, defined as

@@ -402,8 +441,7 @@ all \(\mathbf{r}\). We will check that the local energy, defined as

-is constant. We will also see that when \(a \ne 1\) the local energy -is not constant, so \(\hat{H} \Psi \ne E \Psi\). +is constant.

@@ -413,7 +451,7 @@ with respect to a probability density function \(p(x)\) is given by

-\[ \langle f \rangle_p = \int_{-\infty}^\infty p(x)\, f(x)\,dx \]. +\[ \langle f \rangle_p = \int_{-\infty}^\infty p(x)\, f(x)\,dx. \]

@@ -422,7 +460,7 @@ and integrates to one:

-\[ \int_{-\infty}^\infty p(x)\,dx = 1 \]. +\[ \int_{-\infty}^\infty p(x)\,dx = 1. \]

@@ -441,16 +479,33 @@ E & = & \frac{\langle \Psi| \hat{H} | \Psi\rangle}{\langle \Psi |\Psi \rangle} \end{eqnarray*}
-
-

2.1 Local energy

+
+

2.1 Local energy

+

+Write all the functions of this section in a single file : +hydrogen.py if you use Python, or hydrogen.f90 is you use +Fortran. +

-
-

2.1.1 Exercise 1

+
+

2.1.1 Exercise 1

+Find the theoretical value of \(a\) for which \(\Psi\) is an eigenfunction of \(\hat{H}\). +

+ +
+
+
+ +
+

2.1.2 Exercise 2

+
+
+

Write a function which computes the potential at \(\mathbf{r}\). The function accepts a 3-dimensional vector r as input arguments and returns the potential. @@ -464,10 +519,24 @@ and returns the potential. V(\mathbf{r}) = -\frac{1}{\sqrt{x^2 + y^2 + z^2}} \]

+
-

-Python -

+
+
2.1.2.1 Python
+
+
+
import numpy as np
+
+def potential(r):
+    # TODO
+
+
+
+
+ +
+
2.1.2.2 Python   solution
+
import numpy as np
 
@@ -475,11 +544,26 @@ and returns the potential.
     return -1. / np.sqrt(np.dot(r,r))
+
+
+
+
2.1.2.3 Fortran
+
+
+
double precision function potential(r)
+  implicit none
+  double precision, intent(in) :: r(3)
+  ! TODO
+end function potential
+
+
+
+
-

-Fortran -

+
+
2.1.2.4 Fortran   solution
+
double precision function potential(r)
   implicit none
@@ -490,10 +574,11 @@ and returns the potential.
+
-
-

2.1.2 Exercise 2

-
+
+

2.1.3 Exercise 3

+

Write a function which computes the wave function at \(\mathbf{r}\). @@ -502,20 +587,48 @@ input arguments, and returns a scalar.

+
-

-Python -

+
+
2.1.3.1 Python
+
+
+
def psi(a, r):
+    # TODO
+
+
+
+
+ +
+
2.1.3.2 Python   solution
+
def psi(a, r):
     return np.exp(-a*np.sqrt(np.dot(r,r)))
+
+
-

-Fortran -

+
+
2.1.3.3 Fortran
+
+
+
double precision function psi(a, r)
+  implicit none
+  double precision, intent(in) :: a, r(3)
+  ! TODO
+end function psi
+
+
+
+
+ +
+
2.1.3.4 Fortran   solution
+
double precision function psi(a, r)
   implicit none
@@ -526,10 +639,11 @@ input arguments, and returns a scalar.
+
-
-

2.1.3 Exercise 3

-
+
+

2.1.4 Exercise 4

+

Write a function which computes the local kinetic energy at \(\mathbf{r}\). @@ -584,19 +698,47 @@ So the local kinetic energy is -\frac{1}{2} \frac{\Delta \Psi}{\Psi} (\mathbf{r}) = -\frac{1}{2}\left(a^2 - \frac{2a}{\mathbf{|r|}} \right) \]

+
-

-Python -

+
+
2.1.4.1 Python
+
+
+
def kinetic(a,r):
+    # TODO
+
+
+
+
+ +
+
2.1.4.2 Python   solution
+
def kinetic(a,r):
     return -0.5 * (a**2 - (2.*a)/np.sqrt(np.dot(r,r)))
+
+
-

-Fortran -

+
+
2.1.4.3 Fortran
+
+
+
double precision function kinetic(a,r)
+  implicit none
+  double precision, intent(in) :: a, r(3)
+  ! TODO
+end function kinetic
+
+
+
+
+ +
+
2.1.4.4 Fortran   solution
+
double precision function kinetic(a,r)
   implicit none
@@ -608,15 +750,17 @@ So the local kinetic energy is
+
-
-

2.1.4 Exercise 4

-
+
+

2.1.5 Exercise 5

+

-Write a function which computes the local energy at \(\mathbf{r}\). -The function accepts x,y,z as input arguments and returns the -local energy. +Write a function which computes the local energy at \(\mathbf{r}\), +using the previously defined functions. +The function accepts a and r as input arguments and returns the +local kinetic energy.

@@ -626,20 +770,48 @@ local energy. E_L(\mathbf{r}) = -\frac{1}{2} \frac{\Delta \Psi}{\Psi} (\mathbf{r}) + V(\mathbf{r}) \]

+
-

-Python -

+
+
2.1.5.1 Python
+
+
+
def e_loc(a,r):
+    #TODO
+
+
+
+
+ +
+
2.1.5.2 Python   solution
+
def e_loc(a,r):
     return kinetic(a,r) + potential(r)
+
+
-

-Fortran -

+
+
2.1.5.3 Fortran
+
+
+
double precision function e_loc(a,r)
+  implicit none
+  double precision, intent(in) :: a, r(3)
+  ! TODO
+end function e_loc
+
+
+
+
+ +
+
2.1.5.4 Fortran   solution
+
double precision function e_loc(a,r)
   implicit none
@@ -652,27 +824,31 @@ local energy.
+
-
-

2.2 Plot of the local energy along the \(x\) axis

+
+

2.2 Plot of the local energy along the \(x\) axis

- -
-

2.2.1 Exercise

+
+

2.2.1 Exercise

For multiple values of \(a\) (0.1, 0.2, 0.5, 1., 1.5, 2.), plot the -local energy along the \(x\) axis. +local energy along the \(x\) axis. In Python, you can use matplotlib +for example. In Fortran, it is convenient to write in a text file +the values of \(x\) and \(E_L(\mathbf{r})\) for each point, and use +Gnuplot to plot the files.

+
-

-Python -

+
+
2.2.1.1 Python
+
import numpy as np
 import matplotlib.pyplot as plt
@@ -680,20 +856,36 @@ local energy along the \(x\) axis.
 from hydrogen import e_loc
 
 x=np.linspace(-5,5)
-
-def make_array(a):
-  y=np.array([ e_loc(a, np.array([t,0.,0.]) ) for t in x])
-  return y
-
 plt.figure(figsize=(10,5))
+
+# TODO
+
+plt.tight_layout()
+plt.legend()
+plt.savefig("plot_py.png")
+
+
+
+
+ +
+
2.2.1.2 Python   solution
+
+
+
import numpy as np
+import matplotlib.pyplot as plt
+
+from hydrogen import e_loc
+
+x=np.linspace(-5,5)
+plt.figure(figsize=(10,5))
+
 for a in [0.1, 0.2, 0.5, 1., 1.5, 2.]:
-  y = make_array(a)
+  y=np.array([ e_loc(a, np.array([t,0.,0.]) ) for t in x])
   plt.plot(x,y,label=f"a={a}")
 
 plt.tight_layout()
-
 plt.legend()
-
 plt.savefig("plot_py.png")
@@ -703,12 +895,63 @@ plt.savefig("plot_py.png")

plot_py.png

+
+
+
+
2.2.1.3 Fortran
+
+
+
program plot
+  implicit none
+  double precision, external :: e_loc
 
+  double precision :: x(50), dx
+  integer :: i, j
+
+  dx = 10.d0/(size(x)-1)
+  do i=1,size(x)
+     x(i) = -5.d0 + (i-1)*dx
+  end do
+
+  ! TODO
+
+end program plot
+
+

-Fortran +To compile and run:

+ +
+
gfortran hydrogen.f90 plot_hydrogen.f90 -o plot_hydrogen
+./plot_hydrogen > data
+
+
+ +

+To plot the data using gnuplot: +

+ +
+
set grid
+set xrange [-5:5]
+set yrange [-2:1]
+plot './data' index 0 using 1:2 with lines title 'a=0.1', \
+     './data' index 1 using 1:2 with lines title 'a=0.2', \
+     './data' index 2 using 1:2 with lines title 'a=0.5', \
+     './data' index 3 using 1:2 with lines title 'a=1.0', \
+     './data' index 4 using 1:2 with lines title 'a=1.5', \
+     './data' index 5 using 1:2 with lines title 'a=2.0'
+
+
+
+
+ +
+
2.2.1.4 Fortran   solution
+
program plot
   implicit none
@@ -776,9 +1019,10 @@ plot './data' index 0 using 1:2 with lines title 'a=0.1', \
+
-
-

2.3 Numerical estimation of the energy

+
+

2.3 TODO Numerical estimation of the energy

If the space is discretized in small volume elements \(\mathbf{r}_i\) @@ -808,8 +1052,8 @@ The energy is biased because:

-
-

2.3.1 Exercise

+
+

2.3.1 Exercise

@@ -919,8 +1163,8 @@ a = 2.0000000000000000 E = -8.0869806678448772E-002

-
-

2.4 Variance of the local energy

+
+

2.4 TODO Variance of the local energy

The variance of the local energy is a functional of \(\Psi\) @@ -947,8 +1191,8 @@ energy can be used as a measure of the quality of a wave function.

-
-

2.4.1 Exercise (optional)

+
+

2.4.1 Exercise (optional)

@@ -960,8 +1204,8 @@ Prove that :

-
-

2.4.2 Exercise

+
+

2.4.2 Exercise

@@ -1085,8 +1329,8 @@ a = 2.0000000000000000 E = -8.0869806678448772E-002 s2 = 1.806881

-
-

3 Variational Monte Carlo

+
+

3 TODO Variational Monte Carlo

Numerical integration with deterministic methods is very efficient @@ -1102,8 +1346,8 @@ interval.

-
-

3.1 Computation of the statistical error

+
+

3.1 TODO Computation of the statistical error

To compute the statistical error, you need to perform \(M\) @@ -1143,8 +1387,8 @@ And the confidence interval is given by

-
-

3.1.1 Exercise

+
+

3.1.1 Exercise

@@ -1193,8 +1437,8 @@ input array.

-
-

3.2 Uniform sampling in the box

+
+

3.2 TODO Uniform sampling in the box

We will now do our first Monte Carlo calculation to compute the @@ -1228,8 +1472,8 @@ statistical error.

-
-

3.2.1 Exercise

+
+

3.2.1 Exercise

@@ -1339,8 +1583,8 @@ E = -0.49588321986667677 +/- 7.1758863546737969E-004

-
-

3.3 Metropolis sampling with \(\Psi^2\)

+
+

3.3 TODO Metropolis sampling with \(\Psi^2\)

We will now use the square of the wave function to sample random @@ -1427,8 +1671,8 @@ step such that the acceptance rate is close to 0.5 is a good compromise.

-
-

3.3.1 Exercise

+
+

3.3.1 Exercise

@@ -1558,8 +1802,8 @@ A = 0.51737800000000000 +/- 4.1827406733181444E-004

-
-

3.4 Gaussian random number generator

+
+

3.4 TODO Gaussian random number generator

To obtain Gaussian-distributed random numbers, you can apply the @@ -1613,8 +1857,8 @@ following sections.

-
-

3.5 Generalized Metropolis algorithm

+
+

3.5 TODO Generalized Metropolis algorithm

One can use more efficient numerical schemes to move the electrons. @@ -1713,8 +1957,8 @@ The transition probability becomes:

-
-

3.5.1 Exercise 1

+
+

3.5.1 Exercise 1

@@ -1750,8 +1994,8 @@ Write a function to compute the drift vector \(\frac{\nabla \Psi(\mathbf{r})}{\P

-
-

3.5.2 Exercise 2

+
+

3.5.2 Exercise 2

@@ -1794,8 +2038,8 @@ Modify the previous program to introduce the drifted diffusion scheme. d_old = d_new d2_old = d2_new psi_old = psi_new - N += 1. - E += e_loc(a,r_old) + N += 1. + E += e_loc(a,r_old) return E/N, accep_rate/N @@ -1903,18 +2147,18 @@ A = 0.78861366666666655 +/- 3.5096729498002445E-004

-
-

4 TODO Diffusion Monte Carlo

+
+

4 TODO Diffusion Monte Carlo

-
-

4.1 Hydrogen atom

+
+

4.1 TODO Hydrogen atom

-
    -
  1. Exercise
    -
    +
    +

    4.1.1 Exercise

    +

    Modify the Metropolis VMC program to introduce the PDMC weight. @@ -2067,13 +2311,12 @@ A = 0.78861366666666655 +/- 3.5096729498002445E-004

    -
    2. -
+
-
-

4.2 Dihydrogen

+
+

4.2 TODO Dihydrogen

We will now consider the H2 molecule in a minimal basis composed of the @@ -2095,7 +2338,7 @@ the nuclei.

Author: Anthony Scemama, Claudia Filippi

-

Created: 2021-01-20 Wed 20:19

+

Created: 2021-01-21 Thu 22:25

Validate