qmc-lttc/QMC.org

#+TITLE: Quantum Monte Carlo
#+AUTHOR: Anthony Scemama, Claudia Filippi
#+LANGUAGE:  en
#+INFOJS_OPT: toc:t mouse:underline path:org-info.js
#+STARTUP: latexpreview
#+LATEX_CLASS: report
#+LATEX_HEADER_EXTRA: \usepackage{minted}
#+HTML_HEAD: <link rel="stylesheet" title="Standard" href="worg.css" type="text/css" />

#+OPTIONS: H:4 num:t toc:t \n:nil @:t ::t |:t ^:t -:t f:t *:t <:t
#+OPTIONS: TeX:t LaTeX:t skip:nil d:nil todo:t pri:nil tags:not-in-toc
# EXCLUDE_TAGS: solution solution2 noexport
#+EXCLUDE_TAGS: solution noexport

  #+BEGIN_SRC elisp :output none :exports none
(setq org-latex-listings 'minted
      org-latex-packages-alist '(("" "minted"))
      org-latex-pdf-process
      '("pdflatex -shell-escape -interaction nonstopmode -output-directory %o %f"
        "pdflatex -shell-escape -interaction nonstopmode -output-directory %o %f"
        "pdflatex -shell-escape -interaction nonstopmode -output-directory %o %f"))
(setq org-latex-minted-options '(("breaklines" "true")
                                 ("breakanywhere" "true")))
(setq org-latex-minted-options
      '(("frame" "lines")
        ("fontsize" "\\scriptsize")
        ("linenos" "")))
(org-beamer-export-to-pdf)
                            
  #+END_SRC   

  #+RESULTS:
  : /home/scemama/TREX/qmc-lttc/QMC.pdf

* Introduction

  This website contains the QMC tutorial of the 2021 LTTC winter school
  [[https://www.irsamc.ups-tlse.fr/lttc/Luchon][Tutorials in Theoretical Chemistry]].

  We propose different exercises to understand quantum Monte Carlo (QMC)
  methods. In the first section, we start with the computation of the energy of a
  hydrogen atom using numerical integration. The goal of this section is
  to familarize yourself with the concept of /local energy/.
  Then, we introduce the variational Monte Carlo (VMC) method which
  computes a statistical estimate of the expectation value of the energy
  associated with a given wave function, and apply this approach to the
  hydrogen atom.
  Finally, we present the diffusion Monte Carlo (DMC) method which
  we use here to estimate the exact energy of the hydrogen atom and of the H_2 molecule, 
  starting from an approximate wave function. 

  Code examples will be given in Python3 and Fortran. You can use
  whatever language you prefer to write the programs.

  We consider the stationary solution of the Schrödinger equation, so
  the wave functions considered here are real: for an $N$ electron
  system where the electrons move in the 3-dimensional space,
  $\Psi : \mathbb{R}^{3N} \rightarrow \mathbb{R}$. In addition, $\Psi$
  is defined everywhere, continuous, and infinitely differentiable.

  All the quantities are expressed in /atomic units/ (energies,
  coordinates, etc).

** Energy and local energy

  For a given system with Hamiltonian $\hat{H}$ and wave function $\Psi$, we define the local energy as
  
  $$
  E_L(\mathbf{r}) = \frac{\hat{H} \Psi(\mathbf{r})}{\Psi(\mathbf{r})},
  $$

  where $\mathbf{r}$ denotes the 3N-dimensional electronic coordinates.
  
  The electronic energy of a system, $E$, can be rewritten in terms of the 
  local energy $E_L(\mathbf{r})$ as

  \begin{eqnarray*}
  E & = & \frac{\langle \Psi| \hat{H} | \Psi\rangle}{\langle \Psi |\Psi \rangle} 
      =   \frac{\int \Psi(\mathbf{r})\, \hat{H} \Psi(\mathbf{r})\, d\mathbf{r}}{\int |\Psi(\mathbf{r}) |^2 d\mathbf{r}} \\
    & = & \frac{\int |\Psi(\mathbf{r})|^2\, \frac{\hat{H} \Psi(\mathbf{r})}{\Psi(\mathbf{r})}\,d\mathbf{r}}{\int |\Psi(\mathbf{r}) |^2 d\mathbf{r}} 
      =   \frac{\int |\Psi(\mathbf{r})|^2\, E_L(\mathbf{r})\,d\mathbf{r}}{\int |\Psi(\mathbf{r}) |^2 d\mathbf{r}}  
  \end{eqnarray*}
   
  For few dimensions, one can easily compute $E$ by evaluating the
  integrals on a grid but, for a high number of dimensions, one can
  resort to Monte Carlo techniques to compute $E$.
  
  To this aim, recall that the probabilistic /expected value/ of an
  arbitrary function $f(x)$ with respect to a probability density
  function $P(x)$ is given by

  $$ \langle f \rangle_P = \int_{-\infty}^\infty P(x)\, f(x)\,dx, $$

  where a probability density function $P(x)$ is non-negative
  and integrates to one:

  $$ \int_{-\infty}^\infty P(x)\,dx = 1. $$

  Similarly, we can view the the energy of a system, $E$, as the expected value of the local energy with respect to 
  a probability density $P(\mathbf{r})$ defined in 3$N$ dimensions:
  
  $$ E =  \int E_L(\mathbf{r}) P(\mathbf{r})\,d\mathbf{r} \equiv  \langle E_L \rangle_{P}\,, $$
  
  where the probability density is given by the square of the wave function:
  
  $$ P(\mathbf{r}) = \frac{|\Psi(\mathbf{r})|^2}{\int |\Psi(\mathbf{r})|^2 d\mathbf{r}}\,. $$
  
  If we can sample $N_{\rm MC}$ configurations $\{\mathbf{r}\}$
  distributed as $P$, we can estimate $E$ as the average of the local
  energy computed over these configurations:
  
  $$ E \approx \frac{1}{N_{\rm MC}} \sum_{i=1}^{N_{\rm MC}} E_L(\mathbf{r}_i) \,. $$
  
* Numerical evaluation of the energy of the hydrogen atom

  In this section, we consider the hydrogen atom with the following
  wave function:

  $$
  \Psi(\mathbf{r}) = \exp(-a |\mathbf{r}|)
  $$

  We will first verify that, for a particular value of $a$, $\Psi$ is an
  eigenfunction of the Hamiltonian

  $$
  \hat{H} = \hat{T} + \hat{V} = - \frac{1}{2} \Delta - \frac{1}{|\mathbf{r}|}
  $$

  To do that, we will compute the local energy and check whether it is constant.

** Local energy
   :PROPERTIES:
   :header-args:python: :tangle hydrogen.py
   :header-args:f90: :tangle hydrogen.f90
   :END:

   You will now program all quantities needed to compute the local energy of the H atom for the given wave function.
   
   Write all the functions of this section in a single file :
 ~hydrogen.py~ if you use Python, or ~hydrogen.f90~ is you use
   Fortran.
   
   #+begin_note
   - 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.
   #+end_note
   
*** Exercise 1

    #+begin_exercise
    Write a function which computes the potential at $\mathbf{r}$.
    The function accepts a 3-dimensional vector =r= as input argument
    and returns the potential.
    #+end_exercise

    $\mathbf{r}=\left( \begin{array}{c} x \\ y\\ z\end{array} \right)$, so
    $$
    V(\mathbf{r}) = -\frac{1}{\sqrt{x^2 + y^2 + z^2}}
    $$

 *Python*
     #+BEGIN_SRC python :results none :tangle none
#!/usr/bin/env python3
import numpy as np

def potential(r):
    # TODO
     #+END_SRC

 *Fortran*
     #+BEGIN_SRC f90 :tangle none
double precision function potential(r)
  implicit none
  double precision, intent(in) :: r(3)

  ! TODO

end function potential
     #+END_SRC

**** Solution                                                      :solution2:
 *Python*
     #+BEGIN_SRC python :results none
#!/usr/bin/env python3
import numpy as np

def potential(r):
    distance = np.sqrt(np.dot(r,r))
    assert (distance > 0)
    return -1. / distance
     #+END_SRC

 *Fortran*
     #+BEGIN_SRC f90 
double precision function potential(r)
  implicit none
  double precision, intent(in) :: 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 / distance
  else
     stop 'potential at r=0.d0 diverges'
  end if

end function potential
     #+END_SRC

*** 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
    input arguments, and returns a scalar.
    #+end_exercise
    
 *Python*
     #+BEGIN_SRC python :results none  :tangle none
def psi(a, r):
    # TODO
     #+END_SRC

 *Fortran*
     #+BEGIN_SRC f90  :tangle none
double precision function psi(a, r)
  implicit none
  double precision, intent(in) :: a, r(3)

  ! TODO

end function psi
     #+END_SRC
     
**** Solution                                                      :solution2:
 *Python*
     #+BEGIN_SRC python :results none
def psi(a, r):
    return np.exp(-a*np.sqrt(np.dot(r,r)))
     #+END_SRC

 *Fortran*
     #+BEGIN_SRC f90 
double precision function psi(a, r)
  implicit none
  double precision, intent(in) :: a, r(3)

  psi = dexp(-a * dsqrt( r(1)*r(1) + r(2)*r(2) + r(3)*r(3) ))
end function psi
     #+END_SRC
     
*** 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
    local kinetic energy.
    #+end_exercise

    The local kinetic energy is defined as $$T_L(\mathbf{r}) = -\frac{1}{2}\frac{\Delta \Psi(\mathbf{r})}{\Psi(\mathbf{r})}.$$
     
    We differentiate $\Psi$ with respect to $x$:
     
    \[ \Psi(\mathbf{r})  =  \exp(-a\,|\mathbf{r}|) \]
    \[\frac{\partial \Psi}{\partial x}
      = \frac{\partial \Psi}{\partial |\mathbf{r}|} \frac{\partial |\mathbf{r}|}{\partial x}   
      =  - \frac{a\,x}{|\mathbf{r}|} \Psi(\mathbf{r}) \]

    and we differentiate a second time:

    $$
    \frac{\partial^2 \Psi}{\partial x^2} =
    \left( \frac{a^2\,x^2}{|\mathbf{r}|^2}  -
    \frac{a(y^2+z^2)}{|\mathbf{r}|^{3}} \right) \Psi(\mathbf{r}).
    $$

    The Laplacian operator $\Delta = \frac{\partial^2}{\partial x^2} +
    \frac{\partial^2}{\partial y^2} + \frac{\partial^2}{\partial z^2}$
    applied to the wave function gives:

    $$
    \Delta \Psi (\mathbf{r}) = \left(a^2 - \frac{2a}{\mathbf{|r|}} \right) \Psi(\mathbf{r})\,.
    $$

    Therefore, the local kinetic energy is
    $$
    T_L (\mathbf{r}) = -\frac{1}{2}\left(a^2 - \frac{2a}{\mathbf{|r|}} \right) 
    $$
     
 *Python*
     #+BEGIN_SRC python :results none :tangle none
def kinetic(a,r):
    # TODO
     #+END_SRC

 *Fortran*
     #+BEGIN_SRC f90  :tangle none
double precision function kinetic(a,r)
  implicit none
  double precision, intent(in) :: a, r(3)

  ! TODO

end function kinetic
     #+END_SRC

**** Solution                                                      :solution2:
 *Python*
     #+BEGIN_SRC python :results none
def kinetic(a,r):
    distance = np.sqrt(np.dot(r,r))
    assert (distance > 0.)

    return a * (1./distance - 0.5 * a)
     #+END_SRC

 *Fortran*
     #+BEGIN_SRC f90 
double precision function kinetic(a,r)
  implicit none
  double precision, intent(in) :: a, 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 =  a * (1.d0 / distance - 0.5d0 * a)

  else
     stop 'kinetic energy diverges at r=0'
  end if

end function kinetic
     #+END_SRC

*** Exercise 4
    #+begin_exercise
    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.
    #+end_exercise
   
    $$
    E_L(\mathbf{r}) = -\frac{1}{2} \frac{\Delta \Psi}{\Psi} (\mathbf{r}) + V(\mathbf{r})
    $$

    
 *Python*
     #+BEGIN_SRC python :results none :tangle none
def e_loc(a,r):
    #TODO
     #+END_SRC

 *Fortran*

    #+begin_note
    When you call a function in Fortran, you need to declare its
    return type.
    You might by accident choose a function name which is the
    same as an internal function of Fortran. So it is recommended to
 *always* use the keyword ~external~ to make sure the function you
    are calling is yours.
    #+end_note

    #+BEGIN_SRC f90 :tangle none
double precision function e_loc(a,r)
  implicit none
  double precision, intent(in) :: a, r(3)

  double precision, external :: kinetic
  double precision, external :: potential

  ! TODO

end function e_loc
    #+END_SRC
   
**** Solution                                                      :solution2:
 *Python*
     #+BEGIN_SRC python :results none
def e_loc(a,r):
    return kinetic(a,r) + potential(r)
     #+END_SRC

 *Fortran*
     #+BEGIN_SRC f90
double precision function e_loc(a,r)
  implicit none
  double precision, intent(in) :: a, r(3)

  double precision, external :: kinetic
  double precision, external :: potential

  e_loc = kinetic(a,r) + potential(r)

end function e_loc
     #+END_SRC
   
*** Exercise 5

    #+begin_exercise
    Find the theoretical value of $a$ for which $\Psi$ is an eigenfunction of $\hat{H}$.
    #+end_exercise

**** Solution                                                      :solution2:

  \begin{eqnarray*}
  E &=& \frac{\hat{H} \Psi}{\Psi} = - \frac{1}{2} \frac{\Delta \Psi}{\Psi} -
  \frac{1}{|\mathbf{r}|}  \\
   &=& -\frac{1}{2}\left(a^2 - \frac{2a}{\mathbf{|r|}} \right) -
  \frac{1}{|\mathbf{r}|} \\
   &=&
  -\frac{1}{2} a^2 + \frac{a-1}{\mathbf{|r|}} 
  \end{eqnarray*}

  $a=1$ cancels the $1/|r|$ term, and makes the energy constant and
  equal to -0.5 atomic units.

** Plot of the local energy along the $x$ axis
   :PROPERTIES:
   :header-args:python: :tangle plot_hydrogen.py
   :header-args:f90: :tangle plot_hydrogen.f90
   :END:
   
   The program you will write in this section will be written in
   another file (~plot_hydrogen.py~ or ~plot_hydrogen.f90~ for
   example).
   It will use the functions previously defined.

   In Python, you should put at the beginning of the file
   #+BEGIN_SRC python :results none :tangle none
#!/usr/bin/env python3

from hydrogen import e_loc
   #+END_SRC
   to be able to use the ~e_loc~ function of the ~hydrogen.py~ file.
   
   In Fortran, you will need to compile all the source files together:
   #+begin_src sh :exports both
gfortran hydrogen.f90 plot_hydrogen.f90 -o plot_hydrogen
   #+end_src

*** Exercise

    #+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. 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. With Gnuplot, you will need 2 blank
    lines to separate the data corresponding to different values of $a$.
    #+end_exercise

   #+begin_note
   The potential and the kinetic energy both diverge at $r=0$, so we
   choose a grid which does not contain the origin to avoid numerical issues.
   #+end_note

 *Python*
     #+BEGIN_SRC python :results none :tangle none
#!/usr/bin/env python3

import numpy as np
import matplotlib.pyplot as plt

from hydrogen import e_loc

x=np.linspace(-5,5)
plt.figure(figsize=(10,5))

# TODO

plt.tight_layout()
plt.legend()
plt.savefig("plot_py.png")
     #+end_src

 *Fortran*
     #+begin_src f90  :tangle none
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
     #+end_src

     To compile and run:

     #+begin_src sh :exports both
gfortran hydrogen.f90 plot_hydrogen.f90 -o plot_hydrogen
./plot_hydrogen > data
     #+end_src

     To plot the data using Gnuplot:

     #+begin_src gnuplot :file plot.png :exports code
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'
     #+end_src

**** Solution                                                      :solution2:
 *Python*
     #+BEGIN_SRC python :results none
#!/usr/bin/env python3

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=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")
     #+end_src

     #+RESULTS:

     [[./plot_py.png]]

 *Fortran*
     #+begin_src f90 
program plot
  implicit none
  double precision, external :: e_loc

  double precision :: x(50), energy, dx, r(3), a(6)
  integer :: i, j

  a = (/ 0.1d0, 0.2d0, 0.5d0, 1.d0, 1.5d0, 2.d0 /)

  dx = 10.d0/(size(x)-1)
  do i=1,size(x)
     x(i) = -5.d0 + (i-1)*dx
  end do

  r(:) = 0.d0

  do j=1,size(a)
     print *, '# a=', a(j)
     do i=1,size(x)
        r(1) = x(i)
        energy = e_loc( a(j), r )
        print *, x(i), energy
     end do
     print *, ''
     print *, ''
  end do

end program plot
     #+end_src

     #+begin_src sh :exports none
gfortran hydrogen.f90 plot_hydrogen.f90 -o plot_hydrogen
./plot_hydrogen > data
     #+end_src

     #+begin_src gnuplot :file plot.png :exports results
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'
     #+end_src
     #+RESULTS:
     [[file:plot.png]]

** Numerical estimation of the energy
   :PROPERTIES:
   :header-args:python: :tangle energy_hydrogen.py
   :header-args:f90: :tangle energy_hydrogen.f90
   :END:

   If the space is discretized in small volume elements $\mathbf{r}_i$
   of size $\delta \mathbf{r}$, the expression of $\langle E_L \rangle_{\Psi^2}$
   becomes a weighted average of the local energy, where the weights
   are the values of the wave function square at $\mathbf{r}_i$
   multiplied by the volume element:
     
   $$
   \langle E \rangle_{\Psi^2} \approx \frac{\sum_i w_i E_L(\mathbf{r}_i)}{\sum_i w_i}, \;\;
   w_i = \left|\Psi(\mathbf{r}_i)\right|^2 \delta \mathbf{r}
   $$
     
   #+begin_note
   The energy is biased because:
   - The volume elements are not infinitely small (discretization error)
   - The energy is evaluated only inside the box (incompleteness of the space)
   #+end_note

   
*** Exercise
     #+begin_exercise
    Compute a numerical estimate of the energy using 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 :tangle none
#!/usr/bin/env python3

import numpy as np
from hydrogen import e_loc, psi

interval = np.linspace(-5,5,num=50)
delta = (interval[1]-interval[0])**3

r = np.array([0.,0.,0.])

for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]:
    # TODO
    print(f"a = {a} \t E = {E}")                

     #+end_src

 *Fortran*
     #+begin_src f90 
program energy_hydrogen
  implicit none
  double precision, external :: e_loc, psi
  double precision :: x(50), w, delta, energy, dx, r(3), a(6), norm
  integer :: i, k, l, j

  a = (/ 0.1d0, 0.2d0, 0.5d0, 1.d0, 1.5d0, 2.d0 /)

  dx = 10.d0/(size(x)-1)
  do i=1,size(x)
     x(i) = -5.d0 + (i-1)*dx
  end do

  do j=1,size(a)

     ! TODO

     print *, 'a = ', a(j), '    E = ', energy
  end do

end program energy_hydrogen
     #+end_src

     To compile the Fortran and run it:

     #+begin_src sh :results output :exports code
gfortran hydrogen.f90 energy_hydrogen.f90 -o energy_hydrogen
./energy_hydrogen 
     #+end_src

**** Solution                                                      :solution:
 *Python*
     #+BEGIN_SRC python :results none :exports both
#!/usr/bin/env python3

import numpy as np
from hydrogen import e_loc, psi

interval = np.linspace(-5,5,num=50)
delta = (interval[1]-interval[0])**3

r = np.array([0.,0.,0.])

for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]:
    E    = 0.
    norm = 0.

    for x in interval:
        r[0] = x
        for y in interval:
            r[1] = y
            for z in interval:
                r[2] = z

                w = psi(a,r)
                w = w * w * delta

                E    += w * e_loc(a,r)
                norm += w 

    E = E / norm
    print(f"a = {a} \t E = {E}")                

     #+end_src

     #+RESULTS:
     : a = 0.1 	 E = -0.24518438948809218
     : a = 0.2 	 E = -0.26966057967803525
     : a = 0.5 	 E = -0.3856357612517407
     : a = 0.9 	 E = -0.49435709786716214
     : a = 1.0 	 E = -0.5
     : a = 1.5 	 E = -0.39242967082602226
     : a = 2.0 	 E = -0.08086980667844901

 *Fortran*
     #+begin_src f90 
program energy_hydrogen
  implicit none