[0/1] Last things to do[0/1] Last things to doWe propose different exercises to understand quantum Monte Carlo (QMC) @@ -436,8 +436,8 @@ coordinates, etc).
In this section we consider the Hydrogen atom with the following @@ -510,8 +510,8 @@ E & = & \frac{\langle \Psi| \hat{H} | \Psi\rangle}{\langle \Psi |\Psi \rangle} \end{eqnarray*}
Write all the functions of this section in a single file : @@ -534,8 +534,8 @@ to catch the error.
@@ -546,8 +546,8 @@ Find the theoretical value of \(a\) for which \(\Psi\) is an eigenfunction of \(
@@ -567,7 +567,7 @@ and returns the potential.
import numpy as np @@ -579,7 +579,7 @@ and returns the potential.
import numpy as np @@ -593,7 +593,7 @@ and returns the potential.
double precision function potential(r) @@ -606,7 +606,7 @@ and returns the potential.
double precision function potential(r) @@ -627,8 +627,8 @@ and returns the potential.
@@ -641,7 +641,7 @@ input arguments, and returns a scalar.
def psi(a, r): @@ -651,7 +651,7 @@ input arguments, and returns a scalar.
def psi(a, r): @@ -661,7 +661,7 @@ input arguments, and returns a scalar.
double precision function psi(a, r) @@ -674,7 +674,7 @@ input arguments, and returns a scalar.
double precision function psi(a, r) @@ -689,8 +689,8 @@ input arguments, and returns a scalar.
@@ -749,7 +749,7 @@ So the local kinetic energy is
def kinetic(a,r): @@ -759,7 +759,7 @@ So the local kinetic energy is
def kinetic(a,r): @@ -771,7 +771,7 @@ So the local kinetic energy is
double precision function kinetic(a,r) @@ -784,7 +784,7 @@ So the local kinetic energy is
double precision function kinetic(a,r) @@ -805,8 +805,8 @@ So the local kinetic energy is
@@ -827,7 +827,7 @@ local kinetic energy.
def e_loc(a,r): @@ -837,7 +837,7 @@ local kinetic energy.
def e_loc(a,r): @@ -847,7 +847,7 @@ local kinetic energy.
double precision function e_loc(a,r) @@ -860,7 +860,7 @@ local kinetic energy.
double precision function e_loc(a,r) @@ -877,8 +877,8 @@ local kinetic energy.
@@ -889,8 +889,8 @@ choose a grid which does not contain the origin.
@@ -905,7 +905,7 @@ Gnuplot to plot the files.
import numpy as np @@ -926,7 +926,7 @@ plt.savefig("plot_py.png")
import numpy as np @@ -955,7 +955,7 @@ plt.savefig("plot_py.png")
program plot @@ -1005,7 +1005,7 @@ plot './data' index 0 using 1:2 with lines title 'a=0.1', \
program plot @@ -1077,8 +1077,8 @@ plot './data' index 0 using 1:2 with lines title 'a=0.1', \
If the space is discretized in small volume elements \(\mathbf{r}_i\) @@ -1108,8 +1108,8 @@ The energy is biased because:
@@ -1122,7 +1122,7 @@ Compute a numerical estimate of the energy in a grid of
import numpy as np @@ -1142,7 +1142,7 @@ Compute a numerical estimate of the energy in a grid of
import numpy as np @@ -1174,7 +1174,7 @@ Compute a numerical estimate of the energy in a grid of
program energy_hydrogen @@ -1211,7 +1211,7 @@ To compile the Fortran and run it:
program energy_hydrogen @@ -1280,8 +1280,8 @@ a = 2.0000000000000000 E = -8.0869806678448772E-002
The variance of the local energy is a functional of \(\Psi\) @@ -1308,8 +1308,8 @@ energy can be used as a measure of the quality of a wave function.
@@ -1321,8 +1321,8 @@ Prove that :
@@ -1337,7 +1337,7 @@ in a grid of \(50\times50\times50\) points in the range
import numpy as np @@ -1356,7 +1356,7 @@ in a grid of \(50\times50\times50\) points in the range
import numpy as np @@ -1392,7 +1392,7 @@ in a grid of \(50\times50\times50\) points in the range
program variance_hydrogen @@ -1434,7 +1434,7 @@ To compile and run:
program variance_hydrogen @@ -1511,8 +1511,8 @@ a = 2.0000000000000000 E = -8.0869806678448772E-002 s2 = 1.806881
Numerical integration with deterministic methods is very efficient @@ -1528,8 +1528,8 @@ interval.
To compute the statistical error, you need to perform \(M\) @@ -1569,8 +1569,8 @@ And the confidence interval is given by
@@ -1582,7 +1582,7 @@ input array.
from math import sqrt @@ -1594,7 +1594,7 @@ input array.
from math import sqrt @@ -1613,7 +1613,7 @@ input array.
subroutine ave_error(x,n,ave,err) @@ -1628,7 +1628,7 @@ input array.
subroutine ave_error(x,n,ave,err) @@ -1656,8 +1656,8 @@ input array.
We will now do our first Monte Carlo calculation to compute the @@ -1691,8 +1691,8 @@ compute the statistical error.
@@ -1706,7 +1706,7 @@ compute the average energy and the associated error bar.
@@ -1732,7 +1732,7 @@ the Python solution
+
from hydrogen import * @@ -1759,7 +1759,7 @@ the Fortran
+
@@ -1819,7 +1819,7 @@ well as the index of the current step.
@@ -1892,8 +1892,8 @@ E = -0.49588321986667677 +/- 7.1758863546737969E-004
We will now use the square of the wave function to sample random @@ -1981,8 +1981,8 @@ step such that the acceptance rate is close to 0.5 is a good compromise.
@@ -2003,7 +2003,7 @@ Can you observe a reduction in the statistical error?
from hydrogen import * @@ -2033,7 +2033,7 @@ Can you observe a reduction in the statistical error?
from hydrogen import * @@ -2076,7 +2076,7 @@ Can you observe a reduction in the statistical error?
subroutine metropolis_montecarlo(a,nmax,tau,energy,accep) @@ -2129,7 +2129,7 @@ Can you observe a reduction in the statistical error?
subroutine metropolis_montecarlo(a,nmax,tau,energy,accep) @@ -2209,8 +2209,8 @@ A = 0.51713866666666664 +/- 3.7072551835783688E-004
To obtain Gaussian-distributed random numbers, you can apply the
@@ -2267,8 +2267,8 @@ In Python, you can use the
-
One can use more efficient numerical schemes to move the electrons,
@@ -2367,8 +2367,8 @@ The transition probability becomes:
@@ -2379,7 +2379,7 @@ Write a function to compute the drift vector \(\frac{\nabla \Psi(\mathbf{r})}{\P
@@ -2445,7 +2445,7 @@ Modify the previous program to introduce the drifted diffusion scheme.
Consider the time-dependent Schrödinger equation:
@@ -2713,7 +2713,7 @@ Now, let's replace the time variable \(t\) by an imaginary time variable
\[
- -\frac{\partial \psi(\mathbf{r}, t)}{\partial \tau} = \hat{H} \psi(\mathbf{r}, t)
+ -\frac{\partial \psi(\mathbf{r}, \tau)}{\partial \tau} = \hat{H} \psi(\mathbf{r}, \tau)
\]
+The diffusion equation of particles is given by
+
+\[
+ \frac{\partial \phi(\mathbf{r},t)}{\partial t} = D\, \Delta \phi(\mathbf{r},t).
+ \]
+
+The rate of reaction \(v\) is the speed at which a chemical reaction
+takes place. In a solution, the rate is given as a function of the
+concentration \([A]\) by
+
+\[
+ v = \frac{d[A]}{dt},
+ \]
+
+where the concentration \([A]\) is proportional to the number of particles.
+
+These two equations allow us to interpret the Schrödinger equation
+in imaginary time as the combination of:
+
+The diffusion equation can be simulated by a Brownian motion:
+\[ \mathbf{r}_{n+1} = \mathbf{r}_{n} + \sqrt{\tau}\, \chi \]
+where \(\chi\) is a Gaussian random variable, and the rate equation
+can be simulated by creating or destroying particles over time.
+
+Diffusion Monte Carlo (DMC) consists in obtaining the ground state of a
+system by simulating the Schrödinger equation in imaginary time, by
+the combination of a diffusion process and a rate process.
+
+In a molecular system, the potential is far from being constant,
+and diverges at inter-particle coalescence points. Hence, when the
+rate equation is simulated, it results in very large fluctuations
+in the numbers of particles, making the calculations impossible in
+practice.
+Fortunately, if we multiply the Schrödinger equation by a chosen
+trial wave function \(\Psi_T(\mathbf{r})\) (Hartree-Fock, Kohn-Sham
+determinant, CI wave function, etc), one obtains
+
+\[
+ -\frac{\partial \psi(\mathbf{r},\tau)}{\partial \tau} \Psi_T(\mathbf{r}) =
+ \left[ -\frac{1}{2} \Delta \psi(\mathbf{r},\tau) + V(\mathbf{r}) \psi(\mathbf{r},\tau) \right] \Psi_T(\mathbf{r})
+\]
+
+\[
+-\frac{\partial \big[ \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) \big]}{\partial \tau}
+= -\frac{1}{2} \Big( \Delta \big[
+\psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) \big] -
+\psi(\mathbf{r},\tau) \Delta \Psi_T(\mathbf{r}) - 2
+\nabla \psi(\mathbf{r},\tau) \nabla \Psi_T(\mathbf{r}) \Big) + V(\mathbf{r}) \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r})
+\]
+
+\[
+-\frac{\partial \big[ \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) \big]}{\partial \tau}
+= -\frac{1}{2} \Delta \big[\psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) \big] +
+ \frac{1}{2} \psi(\mathbf{r},\tau) \Delta \Psi_T(\mathbf{r}) +
+\Psi_T(\mathbf{r})\nabla \psi(\mathbf{r},\tau) \frac{\nabla \Psi_T(\mathbf{r})}{\Psi_T(\mathbf{r})} + V(\mathbf{r}) \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r})
+\]
+
+\[
+-\frac{\partial \big[ \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) \big]}{\partial \tau}
+= -\frac{1}{2} \Delta \big[\psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) \big] +
+ \psi(\mathbf{r},\tau) \Delta \Psi_T(\mathbf{r}) +
+\Psi_T(\mathbf{r})\nabla \psi(\mathbf{r},\tau) \frac{\nabla \Psi_T(\mathbf{r})}{\Psi_T(\mathbf{r})} + E_L(\mathbf{r}) \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r})
+\]
+\[
+-\frac{\partial \big[ \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) \big]}{\partial \tau}
+= -\frac{1}{2} \Delta \big[\psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) \big] +
+\nabla \left[ \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r})
+\frac{\nabla \Psi_T(\mathbf{r})}{\Psi_T(\mathbf{r})}
+\right] + E_L(\mathbf{r}) \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r})
+\]
+
+Defining \(\Pi(\mathbf{r},t) = \psi(\mathbf{r},\tau)
+\Psi_T(\mathbf{r})\),
+
+\[
+-\frac{\partial \Pi(\mathbf{r},\tau)}{\partial \tau}
+= -\frac{1}{2} \Delta \Pi(\mathbf{r},\tau) +
+\nabla \left[ \Pi(\mathbf{r},\tau) \frac{\nabla \Psi_T(\mathbf{r})}{\Psi_T(\mathbf{r})}
+\right] + E_L(\mathbf{r}) \Pi(\mathbf{r},\tau)
+\]
+
+The new "potential" is the local energy, which has smaller fluctuations
+as \(\Psi_T\) tends to the exact wave function. The new "kinetic energy"
+can be simulated by the drifted diffusion scheme presented in the
+previous section (VMC).
@@ -2742,11 +2872,12 @@ In the limit \(\tau \rightarrow 0\), you should recover the exact
energy of H for any value of \(a\).
-Python
-
We will now consider the H2 molecule in a minimal basis composed of the
@@ -2917,8 +3048,8 @@ the nuclei.
3.5 Generalized Metropolis algorithm
+3.5 Generalized Metropolis algorithm
3.5.1 Exercise 1
+3.5.1 Exercise 1
-
+
def drift(a,r):
@@ -2389,7 +2389,7 @@ Write a function to compute the drift vector \(\frac{\nabla \Psi(\mathbf{r})}{\P
+
def drift(a,r):
@@ -2400,7 +2400,7 @@ Write a function to compute the drift vector \(\frac{\nabla \Psi(\mathbf{r})}{\P
+
subroutine drift(a,r,b)
@@ -2414,7 +2414,7 @@ Write a function to compute the drift vector \(\frac{\nabla \Psi(\mathbf{r})}{\P
+
subroutine drift(a,r,b)
@@ -2432,8 +2432,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
-
+
from hydrogen import *
@@ -2474,7 +2474,7 @@ Modify the previous program to introduce the drifted diffusion scheme.
+
from hydrogen import *
@@ -2529,7 +2529,7 @@ Modify the previous program to introduce the drifted diffusion scheme.
+
subroutine variational_montecarlo(a,tau,nmax,energy,accep_rate)
@@ -2579,7 +2579,7 @@ Modify the previous program to introduce the drifted diffusion scheme.
+
subroutine variational_montecarlo(a,tau,nmax,energy,accep_rate)
@@ -2672,8 +2672,8 @@ A = 0.78861366666666666 +/- 3.7855335138754813E-004
4 TODO Diffusion Monte Carlo
+4 Diffusion Monte Carlo
+
+
+4.1 TODO Hydrogen atom
+4.1 TODO Hydrogen atom
4.1.1 Exercise
+4.1.1 Exercise
+
+from hydrogen import *
from qmc_stats import *
@@ -2798,10 +2929,10 @@ energy of H for any value of \(a\).
-
-
+subroutine variational_montecarlo(a,tau,nmax,energy,accep_rate)
implicit none
@@ -2895,8 +3026,8 @@ A = 0.78861366666666655 +/- 3.5096729498002445E-004
4.2 TODO Dihydrogen
+4.2 TODO Dihydrogen
5 TODO
+[0/1] Last things to do5 TODO
[0/1] Last things to do
[ ] Prepare 4 questions for the exam: multiple-choice questions
@@ -2933,7 +3064,7 @@ the H\(_2\) molecule at $R$=1.4010 bohr. Answer: 0.17406 a.u.