v1j4y/qmc-lttc
1
0
Fork 0
forked from scemama/qmc-lttc
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Working on VMC

Browse Source
This commit is contained in:
Anthony Scemama 2021-01-12 01:01:52 +01:00
parent 277e243545
commit b6d70add67
9 changed files with 453 additions and 231 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

582
QMC.org
View File

@ -90,9 +90,14 @@
:header-args:python: :tangle hydrogen.py :header-args:python: :tangle hydrogen.py
:header-args:f90: :tangle hydrogen.f90 :header-args:f90: :tangle hydrogen.f90
:END: :END:
*** Write a function which computes the potential at $\mathbf{r}$
*** Exercise 1
#+begin_exercise
Write a function which computes the potential at $\mathbf{r}$.
The function accepts a 3-dimensional vector =r= as input arguments The function accepts a 3-dimensional vector =r= as input arguments
and returns the potential. and returns the potential.
#+end_exercise
$\mathbf{r}=\left( \begin{array}{c} x \\ y\\ z\end{array} \right)$, so $\mathbf{r}=\left( \begin{array}{c} x \\ y\\ z\end{array} \right)$, so
$$ $$
@ -117,9 +122,12 @@ double precision function potential(r)
end function potential end function potential
#+END_SRC #+END_SRC
*** Write a function which computes the wave function at $\mathbf{r}$ *** Exercise 2
#+begin_exercise
Write a function which computes the wave function at $\mathbf{r}$.
The function accepts a scalar =a= and a 3-dimensional vector =r= as The function accepts a scalar =a= and a 3-dimensional vector =r= as
input arguments, and returns a scalar. input arguments, and returns a scalar.
#+end_exercise
*Python* *Python*
@ -137,9 +145,12 @@ double precision function psi(a, r)
end function psi end function psi
#+END_SRC #+END_SRC
*** Write a function which computes the local kinetic energy at $\mathbf{r}$ *** Exercise 3
#+begin_exercise
Write a function which computes the local kinetic energy at $\mathbf{r}$.
The function accepts =a= and =r= as input arguments and returns the The function accepts =a= and =r= as input arguments and returns the
local kinetic energy. local kinetic energy.
#+end_exercise
The local kinetic energy is defined as $$-\frac{1}{2}\frac{\Delta \Psi}{\Psi}.$$ The local kinetic energy is defined as $$-\frac{1}{2}\frac{\Delta \Psi}{\Psi}.$$
@ -187,9 +198,12 @@ double precision function kinetic(a,r)
end function kinetic end function kinetic
#+END_SRC #+END_SRC
*** Write a function which computes the local energy at $\mathbf{r}$ *** Exercise 4
#+begin_exercise
Write a function which computes the local energy at $\mathbf{r}$.
The function accepts =x,y,z= as input arguments and returns the The function accepts =x,y,z= as input arguments and returns the
local energy. local energy.
#+end_exercise
$$ $$
E_L(\mathbf{r}) = -\frac{1}{2} \frac{\Delta \Psi}{\Psi} (\mathbf{r}) + V(\mathbf{r}) E_L(\mathbf{r}) = -\frac{1}{2} \frac{\Delta \Psi}{\Psi} (\mathbf{r}) + V(\mathbf{r})
@ -218,11 +232,15 @@ end function e_loc
:header-args:f90: :tangle plot_hydrogen.f90 :header-args:f90: :tangle plot_hydrogen.f90
:END: :END:
For multiple values of $a$ (0.1, 0.2, 0.5, 1., 1.5, 2.), plot the
local energy along the $x$ axis.
*Python* *** Exercise
#+BEGIN_SRC python :results none #+begin_exercise
For multiple values of $a$ (0.1, 0.2, 0.5, 1., 1.5, 2.), plot the
local energy along the $x$ axis.
#+end_exercise
*Python*
#+BEGIN_SRC python :results none
import numpy as np import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
@ -244,16 +262,16 @@ plt.tight_layout()
plt.legend() plt.legend()
plt.savefig("plot_py.png") plt.savefig("plot_py.png")
#+end_src #+end_src
#+RESULTS: #+RESULTS:
[[./plot_py.png]] [[./plot_py.png]]
*Fortran* *Fortran*
#+begin_src f90 #+begin_src f90
program plot program plot
implicit none implicit none
double precision, external :: e_loc double precision, external :: e_loc
@ -282,20 +300,20 @@ program plot
end do end do
end program plot end program plot
#+end_src #+end_src
To compile and run: To compile and run:
#+begin_src sh :exports both #+begin_src sh :exports both
gfortran hydrogen.f90 plot_hydrogen.f90 -o plot_hydrogen gfortran hydrogen.f90 plot_hydrogen.f90 -o plot_hydrogen
./plot_hydrogen > data ./plot_hydrogen > data
#+end_src #+end_src
#+RESULTS: #+RESULTS:
To plot the data using gnuplot: To plot the data using gnuplot:
#+begin_src gnuplot :file plot.png :exports both #+begin_src gnuplot :file plot.png :exports both
set grid set grid
set xrange [-5:5] set xrange [-5:5]
set yrange [-2:1] set yrange [-2:1]
@ -305,12 +323,12 @@ plot './data' index 0 using 1:2 with lines title 'a=0.1', \
'./data' index 3 using 1:2 with lines title 'a=1.0', \ './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 4 using 1:2 with lines title 'a=1.5', \
'./data' index 5 using 1:2 with lines title 'a=2.0' './data' index 5 using 1:2 with lines title 'a=2.0'
#+end_src #+end_src
#+RESULTS: #+RESULTS:
[[file:plot.png]] [[file:plot.png]]
** Compute numerically the expectation value of the energy ** Numerical estimation of the energy
:PROPERTIES: :PROPERTIES:
:header-args:python: :tangle energy_hydrogen.py :header-args:python: :tangle energy_hydrogen.py
:header-args:f90: :tangle energy_hydrogen.f90 :header-args:f90: :tangle energy_hydrogen.f90
@ -327,18 +345,22 @@ plot './data' index 0 using 1:2 with lines title 'a=0.1', \
w_i = \left[\Psi(\mathbf{r}_i)\right]^2 \delta \mathbf{r} w_i = \left[\Psi(\mathbf{r}_i)\right]^2 \delta \mathbf{r}
$$ $$
In this section, we will compute a numerical estimate of the
energy in a grid of $50\times50\times50$ points in the range
$(-5,-5,-5) \le \mathbf{r} \le (5,5,5)$.
#+begin_note #+begin_note
The energy is biased because: The energy is biased because:
- The volume elements are not infinitely small (discretization error) - The volume elements are not infinitely small (discretization error)
- The energy is evaluated only inside the box (incompleteness of the space) - The energy is evaluated only inside the box (incompleteness of the space)
#+end_note #+end_note
*Python*
#+BEGIN_SRC python :results none *** Exercise
#+begin_exercise
Compute a numerical estimate of the energy in a grid of
$50\times50\times50$ points in the range $(-5,-5,-5) \le
\mathbf{r} \le (5,5,5)$.
#+end_exercise
*Python*
#+BEGIN_SRC python :results none
import numpy as np import numpy as np
from hydrogen import e_loc, psi from hydrogen import e_loc, psi
@ -363,19 +385,19 @@ for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]:
E = E / norm E = E / norm
print(f"a = {a} \t E = {E}") print(f"a = {a} \t E = {E}")
#+end_src #+end_src
#+RESULTS: #+RESULTS:
: a = 0.1 E = -0.24518438948809218 : a = 0.1 E = -0.24518438948809218
: a = 0.2 E = -0.26966057967803525 : a = 0.2 E = -0.26966057967803525
: a = 0.5 E = -0.3856357612517407 : a = 0.5 E = -0.3856357612517407
: a = 0.9 E = -0.49435709786716214 : a = 0.9 E = -0.49435709786716214
: a = 1.0 E = -0.5 : a = 1.0 E = -0.5
: a = 1.5 E = -0.39242967082602226 : a = 1.5 E = -0.39242967082602226
: a = 2.0 E = -0.08086980667844901 : a = 2.0 E = -0.08086980667844901
*Fortran* *Fortran*
#+begin_src f90 #+begin_src f90
program energy_hydrogen program energy_hydrogen
implicit none implicit none
double precision, external :: e_loc, psi double precision, external :: e_loc, psi
@ -414,24 +436,24 @@ program energy_hydrogen
end do end do
end program energy_hydrogen end program energy_hydrogen
#+end_src #+end_src
To compile the Fortran and run it: To compile the Fortran and run it:
#+begin_src sh :results output :exports both #+begin_src sh :results output :exports both
gfortran hydrogen.f90 energy_hydrogen.f90 -o energy_hydrogen gfortran hydrogen.f90 energy_hydrogen.f90 -o energy_hydrogen
./energy_hydrogen ./energy_hydrogen
#+end_src #+end_src
#+RESULTS: #+RESULTS:
: a = 0.10000000000000001 E = -0.24518438948809140 : a = 0.10000000000000001 E = -0.24518438948809140
: a = 0.20000000000000001 E = -0.26966057967803236 : a = 0.20000000000000001 E = -0.26966057967803236
: a = 0.50000000000000000 E = -0.38563576125173815 : a = 0.50000000000000000 E = -0.38563576125173815
: a = 1.0000000000000000 E = -0.50000000000000000 : a = 1.0000000000000000 E = -0.50000000000000000
: a = 1.5000000000000000 E = -0.39242967082602065 : a = 1.5000000000000000 E = -0.39242967082602065
: a = 2.0000000000000000 E = -8.0869806678448772E-002 : a = 2.0000000000000000 E = -8.0869806678448772E-002
** Compute the variance of the local energy ** Variance of the local energy
:PROPERTIES: :PROPERTIES:
:header-args:python: :tangle variance_hydrogen.py :header-args:python: :tangle variance_hydrogen.py
:header-args:f90: :tangle variance_hydrogen.f90 :header-args:f90: :tangle variance_hydrogen.f90
@ -450,8 +472,13 @@ gfortran hydrogen.f90 energy_hydrogen.f90 -o energy_hydrogen
$\hat{H}$) the variance is zero, so the variance of the local $\hat{H}$) the variance is zero, so the variance of the local
energy can be used as a measure of the quality of a wave function. energy can be used as a measure of the quality of a wave function.
*** Exercise
#+begin_exercise
Compute a numerical estimate of the variance of the local energy Compute a numerical estimate of the variance of the local energy
in a grid of $50\times50\times50$ points in the range $(-5,-5,-5) \le \mathbf{r} \le (5,5,5)$. in a grid of $50\times50\times50$ points in the range
$(-5,-5,-5)
\le \mathbf{r} \le (5,5,5)$ for different values of $a$.
#+end_exercise
*Python* *Python*
#+begin_src python :results none #+begin_src python :results none
@ -617,8 +644,11 @@ gfortran hydrogen.f90 variance_hydrogen.f90 -o variance_hydrogen
E \pm \delta E, \text{ where } \delta E = \frac{\sigma}{\sqrt{M}} E \pm \delta E, \text{ where } \delta E = \frac{\sigma}{\sqrt{M}}
$$ $$
*** Exercise
#+begin_exercise
Write a function returning the average and statistical error of an Write a function returning the average and statistical error of an
input array. input array.
#+end_exercise
*Python* *Python*
#+BEGIN_SRC python :results none #+BEGIN_SRC python :results none
@ -656,26 +686,38 @@ end subroutine ave_error
:header-args:f90: :tangle qmc_uniform.f90 :header-args:f90: :tangle qmc_uniform.f90
:END: :END:
In this section we write a function to perform a Monte Carlo We will now do our first Monte Carlo calculation to compute the
calculation of the average energy. energy of the hydrogen atom.
At every Monte Carlo step: At every Monte Carlo step:
- Draw 3 uniform random numbers in the interval $(-5,-5,-5) \le - Draw a random point $\mathbf{r}_i$ in the box $(-5,-5,-5) \le
(x,y,z) \le (5,5,5)$ (x,y,z) \le (5,5,5)$
- Compute $\Psi^2 \times E_L$ at this point and accumulate the - Compute $[\Psi(\mathbf{r}_i)]^2$ and accumulate the result in a
result in E variable =normalization=
- Compute $\Psi^2$ at this point and accumulate the result in N - Compute $[\Psi(\mathbf{r}_i)]^2 \times E_L(\mathbf{r}_i)$, and accumulate the
result in a variable =energy=
Once all the steps have been computed, return the average energy One Monte Carlo run will consist of $N$ Monte Carlo steps. Once all the
computed on the Monte Carlo calculation. steps have been computed, the run returns the average energy
$\bar{E}_k$ over the $N$ steps of the run.
In the main program, write a loop to perform 30 Monte Carlo runs, To compute the statistical error, perform $M$ runs. The final
and compute the average energy and the associated statistical error. estimate of the energy will be the average over the $\bar{E}_k$,
and the variance of the $\bar{E}_k$ will be used to compute the
statistical error.
Compute the energy of the wave function with $a=0.9$. *** Exercise
*Python* #+begin_exercise
#+BEGIN_SRC python :results output Parameterize the wave function with $a=0.9$. Perform 30
independent Monte Carlo runs, each with 100 000 Monte Carlo
steps. Store the final energies of each run and use this array to
compute the average energy and the associated error bar.
#+end_exercise
*Python*
#+BEGIN_SRC python :results output
from hydrogen import * from hydrogen import *
from qmc_stats import * from qmc_stats import *
@ -695,13 +737,13 @@ nmax = 100000
X = [MonteCarlo(a,nmax) for i in range(30)] X = [MonteCarlo(a,nmax) for i in range(30)]
E, deltaE = ave_error(X) E, deltaE = ave_error(X)
print(f"E = {E} +/- {deltaE}") print(f"E = {E} +/- {deltaE}")
#+END_SRC #+END_SRC
#+RESULTS: #+RESULTS:
: E = -0.4956255109300764 +/- 0.0007082875482711226 : E = -0.4956255109300764 +/- 0.0007082875482711226
*Fortran* *Fortran*
#+BEGIN_SRC f90 #+BEGIN_SRC f90
subroutine uniform_montecarlo(a,nmax,energy) subroutine uniform_montecarlo(a,nmax,energy)
implicit none implicit none
double precision, intent(in) :: a double precision, intent(in) :: a
@ -743,15 +785,15 @@ program qmc
call ave_error(X,nruns,ave,err) call ave_error(X,nruns,ave,err)
print *, 'E = ', ave, '+/-', err print *, 'E = ', ave, '+/-', err
end program qmc end program qmc
#+END_SRC #+END_SRC
#+begin_src sh :results output :exports both #+begin_src sh :results output :exports both
gfortran hydrogen.f90 qmc_stats.f90 qmc_uniform.f90 -o qmc_uniform gfortran hydrogen.f90 qmc_stats.f90 qmc_uniform.f90 -o qmc_uniform
./qmc_uniform ./qmc_uniform
#+end_src #+end_src
#+RESULTS: #+RESULTS:
: E = -0.49588321986667677 +/- 7.1758863546737969E-004 : E = -0.49588321986667677 +/- 7.1758863546737969E-004
** Gaussian sampling ** Gaussian sampling
:PROPERTIES: :PROPERTIES:
@ -773,6 +815,9 @@ gfortran hydrogen.f90 qmc_stats.f90 qmc_uniform.f90 -o qmc_uniform
z_2 &=& \sqrt{-2 \ln u_1} \sin(2 \pi u_2) z_2 &=& \sqrt{-2 \ln u_1} \sin(2 \pi u_2)
\end{eqnarray*} \end{eqnarray*}
Here is a Fortran implementation returning a Gaussian-distributed
n-dimensional vector $\mathbf{z}$;
*Fortran* *Fortran*
#+BEGIN_SRC f90 :tangle qmc_stats.f90 #+BEGIN_SRC f90 :tangle qmc_stats.f90
subroutine random_gauss(z,n) subroutine random_gauss(z,n)
@ -805,14 +850,16 @@ end subroutine random_gauss
#+END_SRC #+END_SRC
Now the equation for the energy is changed into Now the sampling probability can be inserted into the equation of the energy:
\[ \[
E = \frac{\int P(\mathbf{r}) \frac{\left[\Psi(\mathbf{r})\right]^2}{P(\mathbf{r})}\, \frac{\hat{H} \Psi(\mathbf{r})}{\Psi(\mathbf{r})}\,d\mathbf{r}}{\int P(\mathbf{r}) \frac{\left[\Psi(\mathbf{r}) \right]^2}{P(\mathbf{r})} d\mathbf{r}} E = \frac{\int P(\mathbf{r}) \frac{\left[\Psi(\mathbf{r})\right]^2}{P(\mathbf{r})}\, \frac{\hat{H} \Psi(\mathbf{r})}{\Psi(\mathbf{r})}\,d\mathbf{r}}{\int P(\mathbf{r}) \frac{\left[\Psi(\mathbf{r}) \right]^2}{P(\mathbf{r})} d\mathbf{r}}
\] \]
with
with the Gaussian probability
\[ \[
P(\mathbf{r}) = \frac{1}{(2 \pi)^{3/2}}\exp\left( -\frac{\mathbf{r}^2}{2} \right) P(\mathbf{r}) = \frac{1}{(2 \pi)^{3/2}}\exp\left( -\frac{\mathbf{r}^2}{2} \right).
\] \]
As the coordinates are drawn with probability $P(\mathbf{r})$, the As the coordinates are drawn with probability $P(\mathbf{r})$, the
@ -823,8 +870,17 @@ end subroutine random_gauss
w_i = \frac{\left[\Psi(\mathbf{r}_i)\right]^2}{P(\mathbf{r}_i)} \delta \mathbf{r} w_i = \frac{\left[\Psi(\mathbf{r}_i)\right]^2}{P(\mathbf{r}_i)} \delta \mathbf{r}
$$ $$
*Python*
#+BEGIN_SRC python :results output *** Exercise
#+begin_exercise
Modify the exercise of the previous section to sample with
Gaussian-distributed random numbers. Can you see an reduction in
the statistical error?
#+end_exercise
*Python*
#+BEGIN_SRC python :results output
from hydrogen import * from hydrogen import *
from qmc_stats import * from qmc_stats import *
@ -848,14 +904,14 @@ nmax = 100000
X = [MonteCarlo(a,nmax) for i in range(30)] X = [MonteCarlo(a,nmax) for i in range(30)]
E, deltaE = ave_error(X) E, deltaE = ave_error(X)
print(f"E = {E} +/- {deltaE}") print(f"E = {E} +/- {deltaE}")
#+END_SRC #+END_SRC
#+RESULTS: #+RESULTS:
: E = -0.49507506093129827 +/- 0.00014164037765553668 : E = -0.49507506093129827 +/- 0.00014164037765553668
*Fortran* *Fortran*
#+BEGIN_SRC f90 #+BEGIN_SRC f90
double precision function gaussian(r) double precision function gaussian(r)
implicit none implicit none
double precision, intent(in) :: r(3) double precision, intent(in) :: r(3)
@ -904,15 +960,15 @@ program qmc
call ave_error(X,nruns,ave,err) call ave_error(X,nruns,ave,err)
print *, 'E = ', ave, '+/-', err print *, 'E = ', ave, '+/-', err
end program qmc end program qmc
#+END_SRC #+END_SRC
#+begin_src sh :results output :exports both #+begin_src sh :results output :exports both
gfortran hydrogen.f90 qmc_stats.f90 qmc_gaussian.f90 -o qmc_gaussian gfortran hydrogen.f90 qmc_stats.f90 qmc_gaussian.f90 -o qmc_gaussian
./qmc_gaussian ./qmc_gaussian
#+end_src #+end_src
#+RESULTS: #+RESULTS:
: E = -0.49606057056767766 +/- 1.3918807547836872E-004 : E = -0.49606057056767766 +/- 1.3918807547836872E-004
** Sampling with $\Psi^2$ ** Sampling with $\Psi^2$
:PROPERTIES: :PROPERTIES:
@ -933,88 +989,105 @@ gfortran hydrogen.f90 qmc_stats.f90 qmc_gaussian.f90 -o qmc_gaussian
E \approx \frac{1}{M}\sum_{i=1}^M E_L(\mathbf{r}_i) E \approx \frac{1}{M}\sum_{i=1}^M E_L(\mathbf{r}_i)
$$ $$
To generate the probability density $\Psi^2$, we consider a
diffusion process characterized by a time-dependent density
$[\Psi(\mathbf{r},t)]^2$, which obeys the Fokker-Planck equation:
\[ *** Importance sampling
\frac{\partial \Psi^2}{\partial t} = \sum_i D
\frac{\partial}{\partial \mathbf{r}_i} \left(
\frac{\partial}{\partial \mathbf{r}_i} - F_i(\mathbf{r}) \right)
[\Psi(\mathbf{r},t)]^2.
\]
$D$ is the diffusion constant and $F_i$ is the i-th component of a To generate the probability density $\Psi^2$, we consider a
drift velocity caused by an external potential. For a stationary diffusion process characterized by a time-dependent density
density, \( \frac{\partial \Psi^2}{\partial t} = 0 \), so $[\Psi(\mathbf{r},t)]^2$, which obeys the Fokker-Planck equation:
\begin{eqnarray*} \[
0 & = & \sum_i D \frac{\partial \Psi^2}{\partial t} = \sum_i D
\frac{\partial}{\partial \mathbf{r}_i} \left( \frac{\partial}{\partial \mathbf{r}_i} \left(
\frac{\partial}{\partial \mathbf{r}_i} - F_i(\mathbf{r}) \right) \frac{\partial}{\partial \mathbf{r}_i} - F_i(\mathbf{r}) \right)
[\Psi(\mathbf{r})]^2 \\ [\Psi(\mathbf{r},t)]^2.
0 & = & \sum_i D \]
\frac{\partial}{\partial \mathbf{r}_i} \left(
\frac{\partial [\Psi(\mathbf{r})]^2}{\partial \mathbf{r}_i} -
F_i(\mathbf{r})\,[\Psi(\mathbf{r})]^2 \right) \\
0 & = &
\frac{\partial^2 \Psi^2}{\partial \mathbf{r}_i^2} -
\frac{\partial F_i }{\partial \mathbf{r}_i}[\Psi(\mathbf{r})]^2 -
\frac{\partial \Psi^2}{\partial \mathbf{r}_i} F_i(\mathbf{r})
\end{eqnarray*}
we search for a drift function which satisfies $D$ is the diffusion constant and $F_i$ is the i-th component of a
drift velocity caused by an external potential. For a stationary
density, \( \frac{\partial \Psi^2}{\partial t} = 0 \), so
\[ \begin{eqnarray*}
\frac{\partial^2 \Psi^2}{\partial \mathbf{r}_i^2} = 0 & = & \sum_i D
[\Psi(\mathbf{r})]^2 \frac{\partial F_i }{\partial \mathbf{r}_i} + \frac{\partial}{\partial \mathbf{r}_i} \left(
\frac{\partial \Psi^2}{\partial \mathbf{r}_i} F_i(\mathbf{r}) \frac{\partial}{\partial \mathbf{r}_i} - F_i(\mathbf{r}) \right)
\] [\Psi(\mathbf{r})]^2 \\
0 & = & \sum_i D
\frac{\partial}{\partial \mathbf{r}_i} \left(
\frac{\partial [\Psi(\mathbf{r})]^2}{\partial \mathbf{r}_i} -
F_i(\mathbf{r})\,[\Psi(\mathbf{r})]^2 \right) \\
0 & = &
\frac{\partial^2 \Psi^2}{\partial \mathbf{r}_i^2} -
\frac{\partial F_i }{\partial \mathbf{r}_i}[\Psi(\mathbf{r})]^2 -
\frac{\partial \Psi^2}{\partial \mathbf{r}_i} F_i(\mathbf{r})
\end{eqnarray*}
to obtain a second derivative on the left, we need the drift to be we search for a drift function which satisfies
of the form
\[
F_i(\mathbf{r}) = g(\mathbf{r}) \frac{\partial \Psi^2}{\partial \mathbf{r}_i}
\]
\[ \[
\frac{\partial^2 \Psi^2}{\partial \mathbf{r}_i^2} = \frac{\partial^2 \Psi^2}{\partial \mathbf{r}_i^2} =
[\Psi(\mathbf{r})]^2 \frac{\partial [\Psi(\mathbf{r})]^2 \frac{\partial F_i }{\partial \mathbf{r}_i} +
g(\mathbf{r})}{\partial \mathbf{r}_i}\frac{\partial \Psi^2}{\partial \mathbf{r}_i} + \frac{\partial \Psi^2}{\partial \mathbf{r}_i} F_i(\mathbf{r})
[\Psi(\mathbf{r})]^2 g(\mathbf{r}) \frac{\partial^2 \]
\Psi^2}{\partial \mathbf{r}_i^2} +
\frac{\partial \Psi^2}{\partial \mathbf{r}_i}
g(\mathbf{r}) \frac{\partial \Psi^2}{\partial \mathbf{r}_i}
\]
$g = 1 / \Psi^2$ satisfies this equation, so to obtain a second derivative on the left, we need the drift to be
of the form
\[
F_i(\mathbf{r}) = g(\mathbf{r}) \frac{\partial \Psi^2}{\partial \mathbf{r}_i}
\]
\[ \[
F(\mathbf{r}) = \frac{\nabla [\Psi(\mathbf{r})]^2}{[\Psi(\mathbf{r})]^2} = 2 \frac{\nabla \frac{\partial^2 \Psi^2}{\partial \mathbf{r}_i^2} =
\Psi(\mathbf{r})}{\Psi(\mathbf{r})} = 2 \nabla \left( \log \Psi(\mathbf{r}) \right) [\Psi(\mathbf{r})]^2 \frac{\partial
\] g(\mathbf{r})}{\partial \mathbf{r}_i}\frac{\partial \Psi^2}{\partial \mathbf{r}_i} +
[\Psi(\mathbf{r})]^2 g(\mathbf{r}) \frac{\partial^2
\Psi^2}{\partial \mathbf{r}_i^2} +
\frac{\partial \Psi^2}{\partial \mathbf{r}_i}
g(\mathbf{r}) \frac{\partial \Psi^2}{\partial \mathbf{r}_i}
\]
following drifted diffusion scheme: $g = 1 / \Psi^2$ satisfies this equation, so
\[ \[
\mathbf{r}_{n+1} = \mathbf{r}_{n} + \tau \frac{\nabla F(\mathbf{r}) = \frac{\nabla [\Psi(\mathbf{r})]^2}{[\Psi(\mathbf{r})]^2} = 2 \frac{\nabla
\Psi(\mathbf{r})}{\Psi(\mathbf{r})} + \eta \sqrt{\tau} \Psi(\mathbf{r})}{\Psi(\mathbf{r})} = 2 \nabla \left( \log \Psi(\mathbf{r}) \right)
\] \]
where $\eta$ is a normally-distributed Gaussian random number. In statistical mechanics, Fokker-Planck trajectories are generated
by a Langevin equation:
\[
\frac{\partial \mathbf{r}(t)}{\partial t} = 2D \frac{\nabla
\Psi(\mathbf{r}(t))}{\Psi} + \eta
\]
First, write a function to compute the drift vector $\frac{\nabla \Psi(\mathbf{r})}{\Psi(\mathbf{r})}$. where $\eta$ is a normally-distributed fluctuating random force.
*Python* Discretizing this differential equation gives the following drifted
#+BEGIN_SRC python diffusion scheme:
\[
\mathbf{r}_{n+1} = \mathbf{r}_{n} + \tau\, 2D \frac{\nabla
\Psi(\mathbf{r})}{\Psi(\mathbf{r})} + \chi
\]
where $\chi$ is a Gaussian random variable with zero mean and
variance $\tau\,2D$.
**** Exercise 1
#+begin_exercise
Write a function to compute the drift vector $\frac{\nabla \Psi(\mathbf{r})}{\Psi(\mathbf{r})}$.
#+end_exercise
*Python*
#+BEGIN_SRC python :tangle hydrogen.py
def drift(a,r): def drift(a,r):
ar_inv = -a/np.sqrt(np.dot(r,r)) ar_inv = -a/np.sqrt(np.dot(r,r))
return r * ar_inv return r * ar_inv
#+END_SRC #+END_SRC
*Fortran* *Fortran*
#+BEGIN_SRC f90 #+BEGIN_SRC f90 :tangle hydrogen.f90
subroutine drift(a,r,b) subroutine drift(a,r,b)
implicit none implicit none
double precision, intent(in) :: a, r(3) double precision, intent(in) :: a, r(3)
@ -1023,54 +1096,50 @@ subroutine drift(a,r,b)
ar_inv = -a / dsqrt(r(1)*r(1) + r(2)*r(2) + r(3)*r(3)) ar_inv = -a / dsqrt(r(1)*r(1) + r(2)*r(2) + r(3)*r(3))
b(:) = r(:) * ar_inv b(:) = r(:) * ar_inv
end subroutine drift end subroutine drift
#+END_SRC #+END_SRC
**** TODO Exercise 2
Now we can write the Monte Carlo sampling: #+begin_exercise
Sample $\Psi^2$ approximately using the drifted diffusion scheme,
with a diffusion constant $D=1/2$. You can use a time step of
0.001 a.u.
#+end_exercise
*Python* *Python*
#+BEGIN_SRC python #+BEGIN_SRC python :results output :tangle vmc.py
from hydrogen import *
from qmc_stats import *
def MonteCarlo(a,tau,nmax): def MonteCarlo(a,tau,nmax):
E = 0. E = 0.
N = 0. N = 0.
sq_tau = sqrt(tau) sq_tau = np.sqrt(tau)
r_old = np.random.normal(loc=0., scale=1.0, size=(3)) r_old = np.random.normal(loc=0., scale=1.0, size=(3))
d_old = drift(a,r_old) d_old = drift(a,r_old)
d2_old = np.dot(d_old,d_old) d2_old = np.dot(d_old,d_old)
psi_old = psi(a,r_old) psi_old = psi(a,r_old)
for istep in range(nmax): for istep in range(nmax):
eta = np.random.normal(loc=0., scale=1.0, size=(3)) chi = np.random.normal(loc=0., scale=1.0, size=(3))
r_new = r_old + tau * d_old + sq_tau * eta r_new = r_old + tau * d_old + sq_tau * chi
d_new = drift(a,r_new) N += 1.
d2_new = np.dot(d_new,d_new) E += e_loc(a,r_new)
psi_new = psi(a,r_new) r_old = r_new
# Metropolis
prod = np.dot((d_new + d_old), (r_new - r_old))
argexpo = 0.5 * (d2_new - d2_old)*tau + prod
q = psi_new / psi_old
q = np.exp(-argexpo) * q*q
if np.random.uniform() < q:
r_old = r_new
d_old = d_new
d2_old = d2_new
psi_old = psi_new
N += 1.
E += e_loc(a,r_old)
return E/N return E/N
a = 0.9
nmax = 100000 nmax = 100000
tau = 0.1 tau = 0.001
X = [MonteCarlo(a,tau,nmax) for i in range(30)] X = [MonteCarlo(a,tau,nmax) for i in range(30)]
E, deltaE = ave_error(X) E, deltaE = ave_error(X)
print(f"E = {E} +/- {deltaE}") print(f"E = {E} +/- {deltaE}")
#+END_SRC #+END_SRC
#+RESULTS: #+RESULTS:
: E = -0.4951783346213532 +/- 0.00022067316984271938 : E = -0.4112049153828464 +/- 0.00027934927432953063
*Fortran* *Fortran*
#+BEGIN_SRC f90 #+BEGIN_SRC f90
subroutine variational_montecarlo(a,nmax,energy) subroutine variational_montecarlo(a,nmax,energy)
implicit none implicit none
double precision, intent(in) :: a double precision, intent(in) :: a
@ -1111,15 +1180,153 @@ program qmc
call ave_error(X,nruns,ave,err) call ave_error(X,nruns,ave,err)
print *, 'E = ', ave, '+/-', err print *, 'E = ', ave, '+/-', err
end program qmc end program qmc
#+END_SRC #+END_SRC
#+begin_src sh :results output :exports both #+begin_src sh :results output :exports both
gfortran hydrogen.f90 qmc_stats.f90 vmc.f90 -o vmc gfortran hydrogen.f90 qmc_stats.f90 vmc.f90 -o vmc
./vmc ./vmc
#+end_src #+end_src
*** Metropolis algorithm
Discretizing the differential equation to generate the desired
probability density will suffer from a discretization error
leading to biases in the averages. The [[https://en.wikipedia.org/wiki/Metropolis%E2%80%93Hastings_algorithm][Metropolis-Hastings
sampling algorithm]] removes exactly the discretization errors, so
large time steps can be employed.
After the new position $\mathbf{r}_{n+1}$ has been computed, an
additional accept/reject step is performed. The acceptance
probability $A$ is chosen so that it is consistent with the
probability of leaving $\mathbf{r}_n$ and the probability of
entering $\mathbf{r}_{n+1}$:
\[ A(\mathbf{r}_{n} \rightarrow \mathbf{r}_{n+1}) = \min \left( 1,
\frac{g(\mathbf{r}_{n+1} \rightarrow \mathbf{r}_{n}) P(\mathbf{r}_{n+1})}
{g(\mathbf{r}_{n} \rightarrow \mathbf{r}_{n+1}) P(\mathbf{r}_{n})}
\right)
\]
In our Hydrogen atom example, $P$ is $\Psi^2$ and $g$ is a
solution of the discretized Fokker-Planck equation:
\begin{eqnarray*}
P(r_{n}) &=& \Psi^2(\mathbf{r}_n) \\
g(\mathbf{r}_{n} \rightarrow \mathbf{r}_{n+1}) & = &
\frac{1}{(4\pi\,D\,\tau)^{3/2}} \exp \left[ - \frac{\left(
\mathbf{r}_{n+1} - \mathbf{r}_{n} - 2D \frac{\nabla
\Psi(\mathbf{r}_n)}{\Psi(\mathbf{r}_n)} \right)^2}{4D\,\tau} \right]
\end{eqnarray*}
The accept/reject step is the following:
- Compute $A(\mathbf{r}_{n} \rightarrow \mathbf{r}_{n+1})$.
- Draw a uniform random number $u$
- if $u \le A(\mathbf{r}_{n} \rightarrow \mathbf{r}_{n+1})$, accept
the move
- if $u>A(\mathbf{r}_{n} \rightarrow \mathbf{r}_{n+1})$, reject
the move: set $\mathbf{r}_{n+1} = \mathbf{r}_{n}$, but *don't
remove the sample from the average!*
* Diffusion Monte Carlo **** TODO Exercise
#+begin_exercise
Modify the previous program to introduce the accept/reject step.
You should recover the unbiased result.
#+end_exercise
*Python*
#+BEGIN_SRC python
def MonteCarlo(a,tau,nmax):
E = 0.
N = 0.
sq_tau = np.sqrt(tau)
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):
eta = np.random.normal(loc=0., scale=1.0, size=(3))
r_new = r_old + tau * d_old + sq_tau * eta
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)*tau + prod
q = psi_new / psi_old
q = np.exp(-argexpo) * q*q
if np.random.uniform() < q:
r_old = r_new
d_old = d_new
d2_old = d2_new
psi_old = psi_new
N += 1.
E += e_loc(a,r_old)
return E/N
nmax = 100000
tau = 0.1
X = [MonteCarlo(a,tau,nmax) for i in range(30)]
E, deltaE = ave_error(X)
print(f"E = {E} +/- {deltaE}")
#+END_SRC
#+RESULTS:
: E = -0.4951783346213532 +/- 0.00022067316984271938
*Fortran*
#+BEGIN_SRC f90
subroutine variational_montecarlo(a,nmax,energy)
implicit none
double precision, intent(in) :: a
integer , 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 variational_montecarlo
program qmc
implicit none
double precision, parameter :: a = 0.9
integer , 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
#+END_SRC
#+begin_src sh :results output :exports both
gfortran hydrogen.f90 qmc_stats.f90 vmc.f90 -o vmc
./vmc
#+end_src
* TODO Diffusion Monte Carlo
We will now consider the H_2 molecule in a minimal basis composed of the We will now consider the H_2 molecule in a minimal basis composed of the
$1s$ orbitals of the hydrogen atoms: $1s$ orbitals of the hydrogen atoms:
@ -1133,3 +1340,4 @@ gfortran hydrogen.f90 qmc_stats.f90 vmc.f90 -o vmc
the nuclei. the nuclei.

1
energy_hydrogen.f90
View File

@ -26,7 +26,6 @@ program energy_hydrogen
r(3) = x(l) r(3) = x(l)
w = psi(a(j),r) w = psi(a(j),r)
w = w * w * delta w = w * w * delta
energy = energy + w * e_loc(a(j), r) energy = energy + w * e_loc(a(j), r)
norm = norm + w norm = norm + w
end do end do

36
energy_hydrogen.py
View File

@ -1,23 +1,23 @@
import numpy as np import numpy as np
from hydrogen import e_loc, psi from hydrogen import e_loc, psi
interval = np.linspace(-5,5,num=50) interval = np.linspace(-5,5,num=50)
delta = (interval[1]-interval[0])**3 delta = (interval[1]-interval[0])**3
r = np.array([0.,0.,0.]) r = np.array([0.,0.,0.])
for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]: for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]:
E = 0. E = 0.
norm = 0. norm = 0.
for x in interval: for x in interval:
r[0] = x r[0] = x
for y in interval: for y in interval:
r[1] = y r[1] = y
for z in interval: for z in interval:
r[2] = z r[2] = z
w = psi(a,r) w = psi(a,r)
w = w * w * delta w = w * w * delta
E += w * e_loc(a,r) E += w * e_loc(a,r)
norm += w norm += w
E = E / norm E = E / norm
print(f"a = {a} \t E = {E}") print(f"a = {a} \t E = {E}")

9
hydrogen.f90
View File

@ -23,3 +23,12 @@ double precision function e_loc(a,r)
double precision, external :: kinetic, potential double precision, external :: kinetic, potential
e_loc = kinetic(a,r) + potential(r) e_loc = kinetic(a,r) + potential(r)
end function e_loc end function e_loc
subroutine drift(a,r,b)
implicit none
double precision, intent(in) :: a, r(3)
double precision, intent(out) :: b(3)
double precision :: ar_inv
ar_inv = -a / dsqrt(r(1)*r(1) + r(2)*r(2) + r(3)*r(3))
b(:) = r(:) * ar_inv
end subroutine drift

4
hydrogen.py
View File

@ -11,3 +11,7 @@ def kinetic(a,r):
def e_loc(a,r): def e_loc(a,r):
return kinetic(a,r) + potential(r) return kinetic(a,r) + potential(r)
def drift(a,r):
ar_inv = -a/np.sqrt(np.dot(r,r))
return r * ar_inv

30
qmc_uniform.py
View File

@ -1,19 +1,19 @@
from hydrogen import * from hydrogen import *
from qmc_stats import * from qmc_stats import *
def MonteCarlo(a, nmax): def MonteCarlo(a, nmax):
E = 0. E = 0.
N = 0. N = 0.
for istep in range(nmax): for istep in range(nmax):
r = np.random.uniform(-5., 5., (3)) r = np.random.uniform(-5., 5., (3))
w = psi(a,r) w = psi(a,r)
w = w*w w = w*w
N += w N += w
E += w * e_loc(a,r) E += w * e_loc(a,r)
return E/N return E/N
a = 0.9 a = 0.9
nmax = 100000 nmax = 100000
X = [MonteCarlo(a,nmax) for i in range(30)] X = [MonteCarlo(a,nmax) for i in range(30)]
E, deltaE = ave_error(X) E, deltaE = ave_error(X)
print(f"E = {E} +/- {deltaE}") print(f"E = {E} +/- {deltaE}")

2
variance_hydrogen.f90
View File

@ -26,7 +26,6 @@ program variance_hydrogen
r(3) = x(l) r(3) = x(l)
w = psi(a(j),r) w = psi(a(j),r)
w = w * w * delta w = w * w * delta
energy = energy + w * e_loc(a(j), r) energy = energy + w * e_loc(a(j), r)
norm = norm + w norm = norm + w
end do end do
@ -44,7 +43,6 @@ program variance_hydrogen
r(3) = x(l) r(3) = x(l)
w = psi(a(j),r) w = psi(a(j),r)
w = w * w * delta w = w * w * delta
s2 = s2 + w * ( e_loc(a(j), r) - energy )**2 s2 = s2 + w * ( e_loc(a(j), r) - energy )**2
norm = norm + w norm = norm + w
end do end do

8
variance_hydrogen.py
View File

@ -20,8 +20,8 @@ for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]:
El = e_loc(a, r) El = e_loc(a, r)
E += w * El E += w * El
norm += w norm += w
E = E / norm E = E / norm
s2 = 0. s2 = 0.
for x in interval: for x in interval:
r[0] = x r[0] = x
for y in interval: for y in interval:
@ -32,5 +32,5 @@ for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]:
w = w * w * delta w = w * w * delta
El = e_loc(a, r) El = e_loc(a, r)
s2 += w * (El - E)**2 s2 += w * (El - E)**2
s2 = s2 / norm s2 = s2 / norm
print(f"a = {a} \t E = {E:10.8f} \t \sigma^2 = {s2:10.8f}") print(f"a = {a} \t E = {E:10.8f} \t \sigma^2 = {s2:10.8f}")

4
worg.css
View File

@ -306,6 +306,10 @@
/* font-lock-warning-face */ /* font-lock-warning-face */
background-color: #e3e3f7; background-color: #e3e3f7;
} }
.exercise {
/* font-lock-warning-face */
background-color: #e3f7e3;
}
.note { .note {
/* font-lock-warning-face */ /* font-lock-warning-face */
background-color: #f7f7d9; background-color: #f7f7d9;