From 46ca82707e79e8f3b7619b2dba62a24845c66b75 Mon Sep 17 00:00:00 2001 From: scemama Date: Mon, 25 Jan 2021 23:23:42 +0000 Subject: [PATCH] deploy: 336cc911c4d1f2f3801995e3ff7bbfde9bd17092 --- index.html | 497 +++++++++++++++++++++++++++++++---------------------- 1 file changed, 291 insertions(+), 206 deletions(-) diff --git a/index.html b/index.html index 2887e5e..363e4ac 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,128 +257,135 @@ for the JavaScript code in this tag.

Table of Contents

-
-

1 Introduction

+
+

1 Introduction

We propose different exercises to understand quantum Monte Carlo (QMC) @@ -413,8 +420,8 @@ coordinates, etc).

-
-

2 Numerical evaluation of the energy

+
+

2 Numerical evaluation of the energy

In this section we consider the Hydrogen atom with the following @@ -487,18 +494,36 @@ 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.

+ +
+
    +
  • When computing a square root in \(\mathbb{R}\), always make sure +that the argument of the square root is non-negative.
    • +
  • When you divide, always make sure that you will not divide by zero
    • +
+ +

+If a floating-point exception can occur, you should make a test +to catch the error. +

+
-
-

2.1.1 Exercise 1

+

+#+endnote +

+
+ +
+

2.1.1 Exercise 1

@@ -509,8 +534,8 @@ Find the theoretical value of \(a\) for which \(\Psi\) is an eigenfunction of \(

-
-

2.1.2 Exercise 2

+
+

2.1.2 Exercise 2

@@ -529,8 +554,8 @@ and returns the potential.

-
-
2.1.2.1 Python
+
+
2.1.2.1 Python
import numpy as np
@@ -542,21 +567,23 @@ and returns the potential.
-
-
2.1.2.2 Python   solution
+
+
2.1.2.2 Python   solution
import numpy as np
 
 def potential(r):
-    return -1. / np.sqrt(np.dot(r,r))
+    distance = np.sqrt(np.dot(r,r))
+    assert (distance > 0)
+    return -1. / distance
-
-
2.1.2.3 Fortran
+
+
2.1.2.3 Fortran
double precision function potential(r)
@@ -569,14 +596,20 @@ and returns the potential.
-
-
2.1.2.4 Fortran   solution
+
+
2.1.2.4 Fortran   solution
double precision function potential(r)
   implicit none
   double precision, intent(in) :: r(3)
-  potential = -1.d0 / dsqrt( r(1)*r(1) + r(2)*r(2) + r(3)*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) )
+  else
+     stop 'potential at r=0.d0 diverges'
+  end if
 end function potential
@@ -584,8 +617,8 @@ and returns the potential.
-
-

2.1.3 Exercise 3

+
+

2.1.3 Exercise 3

@@ -597,9 +630,8 @@ input arguments, and returns a scalar.

- -
-
2.1.3.1 Python
+
+
2.1.3.1 Python
def psi(a, r):
@@ -609,8 +641,8 @@ input arguments, and returns a scalar.
-
-
2.1.3.2 Python   solution
+
+
2.1.3.2 Python   solution
def psi(a, r):
@@ -620,8 +652,8 @@ input arguments, and returns a scalar.
-
-
2.1.3.3 Fortran
+
+
2.1.3.3 Fortran
double precision function psi(a, r)
@@ -634,8 +666,8 @@ input arguments, and returns a scalar.
-
-
2.1.3.4 Fortran   solution
+
+
2.1.3.4 Fortran   solution
double precision function psi(a, r)
@@ -649,8 +681,8 @@ input arguments, and returns a scalar.
-
-

2.1.4 Exercise 4

+
+

2.1.4 Exercise 4

@@ -708,8 +740,8 @@ So the local kinetic energy is

-
-
2.1.4.1 Python
+
+
2.1.4.1 Python
def kinetic(a,r):
@@ -719,19 +751,21 @@ So the local kinetic energy is
-
-
2.1.4.2 Python   solution
+
+
2.1.4.2 Python   solution
def kinetic(a,r):
-    return -0.5 * (a**2 - (2.*a)/np.sqrt(np.dot(r,r)))
+    distance = np.sqrt(np.dot(r,r))
+    assert (distance > 0.)
+    return -0.5 * (a**2 - (2.*a)/distance)
-
-
2.1.4.3 Fortran
+
+
2.1.4.3 Fortran
double precision function kinetic(a,r)
@@ -744,15 +778,20 @@ So the local kinetic energy is
-
-
2.1.4.4 Fortran   solution
+
+
2.1.4.4 Fortran   solution
double precision function kinetic(a,r)
   implicit none
   double precision, intent(in) :: a, r(3)
-  kinetic = -0.5d0 * (a*a - (2.d0*a) / &
-       dsqrt( r(1)*r(1) + r(2)*r(2) + r(3)*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)
+  else
+     stop 'kinetic energy diverges at r=0'
+  end if
 end function kinetic
@@ -760,8 +799,8 @@ So the local kinetic energy is
-
-

2.1.5 Exercise 5

+
+

2.1.5 Exercise 5

@@ -781,8 +820,8 @@ local kinetic energy.

-
-
2.1.5.1 Python
+
+
2.1.5.1 Python
def e_loc(a,r):
@@ -792,8 +831,8 @@ local kinetic energy.
-
-
2.1.5.2 Python   solution
+
+
2.1.5.2 Python   solution
def e_loc(a,r):
@@ -803,8 +842,8 @@ local kinetic energy.
-
-
2.1.5.3 Fortran
+
+
2.1.5.3 Fortran
double precision function e_loc(a,r)
@@ -817,8 +856,8 @@ local kinetic energy.
-
-
2.1.5.4 Fortran   solution
+
+
2.1.5.4 Fortran   solution
double precision function e_loc(a,r)
@@ -834,13 +873,20 @@ local kinetic energy.
-
-

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

+
+

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

+
+

+The potential and the kinetic energy both diverge at \(r=0\), so we +choose a grid which does not contain the origin. +

+ +
-
-

2.2.1 Exercise

+
+

2.2.1 Exercise

@@ -854,8 +900,8 @@ Gnuplot to plot the files.

-
-
2.2.1.1 Python
+
+
2.2.1.1 Python
import numpy as np
@@ -876,8 +922,8 @@ plt.savefig("plot_py.png")
-
-
2.2.1.2 Python   solution
+
+
2.2.1.2 Python   solution
import numpy as np
@@ -906,8 +952,8 @@ plt.savefig("plot_py.png")
-
-
2.2.1.3 Fortran
+
+
2.2.1.3 Fortran
program plot
@@ -957,8 +1003,8 @@ plot './data' index 0 using 1:2 with lines title 'a=0.1', \
-
-
2.2.1.4 Fortran   solution
+
+
2.2.1.4 Fortran   solution
program plot
@@ -1029,8 +1075,8 @@ plot './data' index 0 using 1:2 with lines title 'a=0.1', \
-
-

2.3 Numerical estimation of the energy

+
+

2.3 Numerical estimation of the energy

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

-
-

2.3.1 Exercise

+
+

2.3.1 Exercise

@@ -1073,8 +1119,8 @@ Compute a numerical estimate of the energy in a grid of

-
-
2.3.1.1 Python
+
+
2.3.1.1 Python
import numpy as np
@@ -1094,8 +1140,8 @@ Compute a numerical estimate of the energy in a grid of
-
-
2.3.1.2 Python   solution
+
+
2.3.1.2 Python   solution
import numpy as np
@@ -1127,8 +1173,8 @@ Compute a numerical estimate of the energy in a grid of
-
-
2.3.1.3 Fortran
+
+
2.3.1.3 Fortran
program energy_hydrogen
@@ -1165,8 +1211,8 @@ To compile the Fortran and run it:
-
-
2.3.1.4 Fortran   solution
+
+
2.3.1.4 Fortran   solution
program energy_hydrogen
@@ -1234,8 +1280,8 @@ a =    2.0000000000000000          E =   -8.0869806678448772E-002
-
-

2.4 Variance of the local energy

+
+

2.4 Variance of the local energy

The variance of the local energy is a functional of \(\Psi\) @@ -1262,8 +1308,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)

@@ -1275,8 +1321,8 @@ Prove that :

-
-

2.4.2 Exercise

+
+

2.4.2 Exercise

@@ -1290,8 +1336,8 @@ in a grid of \(50\times50\times50\) points in the range

-
-
2.4.2.1 Python
+
+
2.4.2.1 Python
import numpy as np
@@ -1310,8 +1356,8 @@ in a grid of \(50\times50\times50\) points in the range
-
-
2.4.2.2 Python   solution
+
+
2.4.2.2 Python   solution
import numpy as np
@@ -1347,8 +1393,8 @@ in a grid of \(50\times50\times50\) points in the range
-
-
2.4.2.3 Fortran
+
+
2.4.2.3 Fortran
program variance_hydrogen
@@ -1390,8 +1436,8 @@ To compile and run:
-
-
2.4.2.4 Fortran   solution
+
+
2.4.2.4 Fortran   solution
program variance_hydrogen
@@ -1467,8 +1513,8 @@ a =    2.0000000000000000       E =   -8.0869806678448772E-002  s2 =    1.806881
-
-

3 TODO Variational Monte Carlo

+
+

3 Variational Monte Carlo

Numerical integration with deterministic methods is very efficient @@ -1484,14 +1530,14 @@ interval.

-
-

3.1 TODO Computation of the statistical error

+
+

3.1 Computation of the statistical error

To compute the statistical error, you need to perform \(M\) independent Monte Carlo calculations. You will obtain \(M\) different estimates of the energy, which are expected to have a Gaussian -distribution by the central limit theorem. +distribution according to the Central Limit Theorem.

@@ -1525,8 +1571,8 @@ And the confidence interval is given by

-
-

3.1.1 Exercise

+
+

3.1.1 Exercise

@@ -1534,25 +1580,61 @@ Write a function returning the average and statistical error of an input array.

+
-

-Python -

+
+
3.1.1.1 Python
+
+
+
from math import sqrt
+def ave_error(arr):
+    #TODO
+    return (average, error)
+
+
+
+
+ +
+
3.1.1.2 Python   solution
+
from math import sqrt
 def ave_error(arr):
     M = len(arr)
-    assert (M>1)
-    average = sum(arr)/M
-    variance = 1./(M-1) * sum( [ (x - average)**2 for x in arr ] )
-    return (average, sqrt(variance/M))
+    assert(M>0)
+    if M == 1:
+        return (arr[0], 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)
+
+
-

-Fortran -

+
+
3.1.1.3 Fortran
+
+
+
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
+
+
+
+
+ +
+
3.1.1.4 Fortran   solution
+
subroutine ave_error(x,n,ave,err)
   implicit none
@@ -1560,7 +1642,9 @@ input array.
   double precision, intent(in)  :: x(n) 
   double precision, intent(out) :: ave, err
   double precision :: variance
-  if (n == 1) then
+  if (n < 1) then
+     stop 'n<1 in ave_error'
+  else if (n == 1) then
      ave = x(1)
      err = 0.d0
   else
@@ -1574,9 +1658,10 @@ input array.
+
-
-

3.2 TODO 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 @@ -1610,8 +1695,8 @@ statistical error.

-
-

3.2.1 Exercise

+
+

3.2.1 Exercise

@@ -1721,8 +1806,8 @@ E = -0.49588321986667677 +/- 7.1758863546737969E-004

-
-

3.3 TODO 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 @@ -1809,8 +1894,8 @@ step such that the acceptance rate is close to 0.5 is a good compromise.

-
-

3.3.1 Exercise

+
+

3.3.1 Exercise

@@ -1940,8 +2025,8 @@ A = 0.51737800000000000 +/- 4.1827406733181444E-004

-
-

3.4 TODO Gaussian random number generator

+
+

3.4 TODO Gaussian random number generator

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

-
-

3.5 TODO Generalized Metropolis algorithm

+
+

3.5 TODO Generalized Metropolis algorithm

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

-
-

3.5.1 Exercise 1

+
+

3.5.1 Exercise 1

@@ -2132,8 +2217,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

@@ -2285,17 +2370,17 @@ A = 0.78861366666666655 +/- 3.5096729498002445E-004

-
-

4 TODO Diffusion Monte Carlo

+
+

4 TODO Diffusion Monte Carlo

-
-

4.1 TODO Hydrogen atom

+
+

4.1 TODO Hydrogen atom

-
-

4.1.1 Exercise

+
+

4.1.1 Exercise

@@ -2453,8 +2538,8 @@ A = 0.78861366666666655 +/- 3.5096729498002445E-004

-
-

4.2 TODO Dihydrogen

+
+

4.2 TODO Dihydrogen

We will now consider the H2 molecule in a minimal basis composed of the @@ -2475,8 +2560,8 @@ the nuclei.

-
-

5 TODO [0/1] Last things to do

+
+

5 TODO [0/1] Last things to do

  • [ ] Prepare 4 questions for the exam: multiple-choice questions @@ -2491,7 +2576,7 @@ the H\(_2\) molecule at $R$=1.4010 bohr. Answer: 0.17406 a.u.

Author: Anthony Scemama, Claudia Filippi

-

Created: 2021-01-25 Mon 22:54

+

Created: 2021-01-25 Mon 23:23

Validate